Phylogenetic Generalized Linear Mixed Model for Community Data

This function performs Generalized Linear Mixed Models for binary, count, and continuous data, estimating regression coefficients with approximate standard errors. It is specifically designed for community data in which species occur within multiple sites (locations). A Bayesian version of PGLMM uses the package INLA, which is not available on CRAN yet. If you wish to use this option, you must first install INLA from https://www.r-inla.org/ by running install.packages('INLA', repos='https://www.math.ntnu.no/inla/R/stable') in R.

pglmm(
  formula,
  data = NULL,
  family = "gaussian",
  cov_ranef = NULL,
  random.effects = NULL,
  REML = TRUE,
  optimizer = c("nelder-mead-nlopt", "bobyqa", "Nelder-Mead", "subplex"),
  repulsion = FALSE,
  add.obs.re = TRUE,
  verbose = FALSE,
  cpp = TRUE,
  bayes = FALSE,
  s2.init = NULL,
  B.init = NULL,
  reltol = 10^-6,
  maxit = 500,
  tol.pql = 10^-6,
  maxit.pql = 200,
  marginal.summ = "mean",
  calc.DIC = TRUE,
  calc.WAIC = TRUE,
  prior = "inla.default",
  prior_alpha = 0.1,
  prior_mu = 1,
  ML.init = FALSE,
  tree = NULL,
  tree_site = NULL,
  sp = NULL,
  site = NULL,
  bayes_options = NULL,
  bayes_nested_matrix_as_list = FALSE
)

communityPGLMM(
  formula,
  data = NULL,
  family = "gaussian",
  cov_ranef = NULL,
  random.effects = NULL,
  REML = TRUE,
  optimizer = c("nelder-mead-nlopt", "bobyqa", "Nelder-Mead", "subplex"),
  repulsion = FALSE,
  add.obs.re = TRUE,
  verbose = FALSE,
  cpp = TRUE,
  bayes = FALSE,
  s2.init = NULL,
  B.init = NULL,
  reltol = 10^-6,
  maxit = 500,
  tol.pql = 10^-6,
  maxit.pql = 200,
  marginal.summ = "mean",
  calc.DIC = TRUE,
  calc.WAIC = TRUE,
  prior = "inla.default",
  prior_alpha = 0.1,
  prior_mu = 1,
  ML.init = FALSE,
  tree = NULL,
  tree_site = NULL,
  sp = NULL,
  site = NULL,
  bayes_options = NULL,
  bayes_nested_matrix_as_list = FALSE
)

Arguments

formula

A two-sided linear formula object describing the mixed effects of the model.

To specify that a random term should have phylogenetic covariance matrix along with non-phylogenetic one, add __ (two underscores) at the end of the group variable; e.g., + (1 | sp__) will construct two random terms, one with phylogenetic covariance matrix and another with non-phylogenetic (identity) matrix. In contrast, __ in the nested terms (below) will only create a phylogenetic covariance matrix. Nested random terms have the general form (1|sp__@site__) which represents phylogenetically related species nested within correlated sites. This form can be used for bipartite questions. For example, species could be phylogenetically related pollinators and sites could be phylogenetically related plants, leading to the random effect (1|insects__@plants__). If more than one phylogeny is used, remember to add all to the argument cov_ranef = list(insects = insect_phylo, plants = plant_phylo). Phylogenetic correlations can be dropped by removing the __ underscores. Thus, the form (1|sp@site__) excludes the phylogenetic correlations among species, while the form (1|sp__@site) excludes the correlations among sites.

data

A data.frame containing the variables named in formula.

family

Either "gaussian" for a Linear Mixed Model, or "binomial" or "poisson" for Generalized Linear Mixed Models. "family" should be specified as a character string (i.e., quoted). For binomial and Poisson data, we use the canonical logit and log link functions, respectively. Binomial data can be either presence/absence, or a two-column array of 'successes' and 'failures'. For both binomial and Poisson data, we add an observation-level random term by default via add.obs.re = TRUE. If bayes = TRUE there are two additional families available: "zeroinflated.binomial", and "zeroinflated.poisson", which add a zero inflation parameter; this parameter gives the probability that the response is a zero. The rest of the parameters of the model then reflect the "non-zero" part part of the model. Note that "zeroinflated.binomial" only makes sense for success/failure response data.

cov_ranef

A named list of covariance matrices of random terms. The names should be the group variables that are used as random terms with specified covariance matrices (without the two underscores, e.g. list(sp = tree1, site = tree2)). The actual object can be either a phylogeny with class "phylo" or a prepared covariance matrix. If it is a phylogeny, pglmm will prune it and then convert it to a covariance matrix assuming Brownian motion evolution. pglmm will also standardize all covariance matrices to have determinant of one. Group variables will be converted to factors and all covariance matrices will be rearranged so that rows and columns are in the same order as the levels of their corresponding group variables.

random.effects

Optional pre-build list of random effects. If NULL (the default), the function prep_dat_pglmm will prepare the random effects for you from the information in formula, data, and cov_ranef. random.effect allows a list of pre-generated random effects terms to increase flexibility; for example, this makes it possible to construct models with both phylogenetic correlation and spatio-temporal autocorrelation. In preparing random.effect, make sure that the orders of rows and columns of covariance matrices in the list are the same as their corresponding group variables in the data. Also, this should be a list of lists, e.g. random.effects = list(re1 = list(matrix_a), re2 = list(1, sp = sp, covar = Vsp)).

REML

Whether REML or ML is used for model fitting the random effects. Ignored if bayes = TRUE.

optimizer

nelder-mead-nlopt (default), bobyqa, Nelder-Mead, or subplex. Nelder-Mead is from the stats package and the other optimizers are from the nloptr package. Ignored if bayes = TRUE.

repulsion

When there are nested random terms specified, repulsion = FALSE tests for phylogenetic underdispersion while repulsion = FALSE tests for overdispersion. This argument is a logical vector of length either 1 or >1. If its length is 1, then all covariance matrices in nested terms will be either inverted (overdispersion) or not. If its length is >1, then you can select which covariance matrix in the nested terms to be inverted. Make sure to get the length right: for all the terms with @, count the number of "__" to determine the length of repulsion. For example, sp__@site and sp@site__ will each require one element of repulsion, while sp__@site__ will take two elements (repulsion for sp and repulsion for site). Therefore, if your nested terms are (1|sp__@site) + (1|sp@site__) + (1|sp__@site__), then you should set the repulsion to be something like c(TRUE, FALSE, TRUE, TRUE) (length of 4).

add.obs.re

Whether to add an observation-level random term for binomial or Poisson distributions. Normally it would be a good idea to add this to account for overdispersion, so add.obs.re = TRUE by default.

verbose

If TRUE, the model deviance and running estimates of s2 and B are plotted each iteration during optimization.

cpp

Whether to use C++ function for optim. Default is TRUE. Ignored if bayes = TRUE.

bayes

Whether to fit a Bayesian version of the PGLMM using r-inla.

s2.init

An array of initial estimates of s2 for each random effect that scales the variance. If s2.init is not provided for family="gaussian", these are estimated using lm assuming no phylogenetic signal. A better approach might be to run link[lme4:lmer]{lmer} and use the output random effects for s2.init. If s2.init is not provided for family = "binomial", these are set to 0.25.

B.init

Initial estimates of $B$, a matrix containing regression coefficients in the model for the fixed effects. This matrix must have dim(B.init) = c(p + 1, 1), where p is the number of predictor (independent) variables; the first element of B corresponds to the intercept, and the remaining elements correspond in order to the predictor (independent) variables in the formula. If B.init is not provided, these are estimated using lm or glm assuming no phylogenetic signal. A better approach might be to run lmer and use the output fixed effects for B.init. When bayes = TRUE, initial values are estimated using the maximum likelihood fit unless ML.init = FALSE, in which case the default INLA initial values will be used.

reltol

A control parameter dictating the relative tolerance for convergence in the optimization; see optim.

maxit

A control parameter dictating the maximum number of iterations in the optimization; see optim.

tol.pql

A control parameter dictating the tolerance for convergence in the PQL estimates of the mean components of the GLMM. Ignored if family = "gaussian" or bayes = TRUE.

maxit.pql

A control parameter dictating the maximum number of iterations in the PQL estimates of the mean components of the GLMM. Ignored if family = "gaussian" or bayes = TRUE.

marginal.summ

Summary statistic to use for the estimate of coefficients when doing a Bayesian PGLMM (when bayes = TRUE). Options are: "mean", "median", or "mode", referring to different characterizations of the central tendency of the Bayesian posterior marginal distributions. Ignored if bayes = FALSE.

calc.DIC

Should the Deviance Information Criterion be calculated and returned when doing a Bayesian PGLMM? Ignored if bayes = FALSE.

calc.WAIC

Should the WAIC be calculated and returned when doing a Bayesian PGLMM? Ignored if bayes = FALSE.

prior

Which type of default prior should be used by pglmm? Only used if bayes = TRUE. There are currently four options: "inla.default", which uses the default INLA priors; "pc.prior.auto", which uses a complexity penalizing prior (as described in Simpson et al. (2017)) designed to automatically choose good parameters (only available for gaussian and binomial responses); "pc.prior", which allows the user to set custom parameters on the "pc.prior" prior, using the prior_alpha and prior_mu parameters (Run INLA::inla.doc("pc.prec") for details on these parameters); and "uninformative", which sets a very uninformative prior (nearly uniform) by using a very flat exponential distribution. The last option is generally not recommended but may in some cases give estimates closer to the maximum likelihood estimates. "pc.prior.auto" is only implemented for family = "gaussian" and family = "binomial" currently.

prior_alpha

Only used if bayes = TRUE and prior = "pc.prior", in which case it sets the alpha parameter of INLA's complexity penalizing prior for the random effects. The prior is an exponential distribution where prob(sd > mu) = alpha, where sd is the standard deviation of the random effect.

prior_mu

Only used if bayes = TRUE and prior = "pc.prior", in which case it sets the mu parameter of INLA's complexity penalizing prior for the random effects. The prior is an exponential distribution where prob(sd > mu) = alpha, where sd is the standard deviation of the random effect.

ML.init

Only relevant if bayes = TRUE. Should maximum likelihood estimates be calculated and used as initial values for the Bayesian model fit? Sometimes this can be helpful, but it may not help; thus, we set the default to FALSE. Also, it does not work with the zero-inflated families.

tree

A phylogeny for column sp, with "phylo" class, or a covariance matrix for sp. Make sure to have all species in the matrix; if the matrix is not standardized, (i.e., det(tree) != 1), pglmm will try to standardize it for you. No longer used: keep here for compatibility.

tree_site

A second phylogeny for "site". This is required only if the site column contains species instead of sites. This can be used for bipartitie questions; tree_site can also be a covariance matrix. Make sure to have all sites in the matrix; if the matrix is not standardized (i.e., det(tree_site) != 1), pglmm` will try to standardize it for you. No longer used: keep here for compatibility.

sp

No longer used: keep here for compatibility.

site

No longer used: keep here for compatibility.

bayes_options

Additional options to pass to INLA for if bayes = TRUE. A named list where the names correspond to parameters in the inla function. One special option is diagonal: if an element in the options list is names diagonal this tells INLA to add its value to the diagonal of the random effects precision matrices. This can help with numerical stability if the model is ill-conditioned (if you get a lot of warnings, try setting this to list(diagonal = 1e-4)).

bayes_nested_matrix_as_list

For bayes = TRUE, prepare the nested terms as a list of length of 4 as the old way?

Value

An object (list) of class communityPGLMM with the following elements:

formula: the formula for fixed effects
formula_original: the formula for both fixed effects and random effects
data: the dataset
family: gaussian, binomial, or poisson depending on the model fit
random.effects: the list of random effects
B: estimates of the regression coefficients
B.se: approximate standard errors of the fixed effects regression coefficients. This is set to NULL if bayes = TRUE.
B.ci: approximate Bayesian credible interval of the fixed effects regression coefficients. This is set to NULL if bayes = FALSE
B.cov: approximate covariance matrix for the fixed effects regression coefficients
B.zscore: approximate Z scores for the fixed effects regression coefficients. This is set to NULL if bayes = TRUE
B.pvalue: approximate tests for the fixed effects regression coefficients being different from zero. This is set to NULL if bayes = TRUE
ss: standard deviations of the random effects for the covariance matrix $\sigma^2V$ for each random effect in order. For the linear mixed model, the residual variance is listed last.
s2r: random effects variances for non-nested random effects
s2n: random effects variances for nested random effects
s2resid: for linear mixed models, the residual variance
s2r.ci: Bayesian credible interval for random effects variances for non-nested random effects. This is set to NULL if bayes = FALSE
s2n.ci: Bayesian credible interval for random effects variances for nested random effects. This is set to NULL if bayes = FALSE
s2resid.ci: Bayesian credible interval for linear mixed models, the residual variance. This is set to NULL if bayes = FALSE
logLik: for linear mixed models, the log-likelihood for either the restricted likelihood (REML=TRUE) or the overall likelihood (REML=FALSE). This is set to NULL for generalized linear mixed models. If bayes = TRUE, this is the marginal log-likelihood
AIC: for linear mixed models, the AIC for either the restricted likelihood (REML = TRUE) or the overall likelihood (REML = FALSE). This is set to NULL for generalised linear mixed models
BIC: for linear mixed models, the BIC for either the restricted likelihood (REML = TRUE) or the overall likelihood (REML = FALSE). This is set to NULL for generalised linear mixed models
DIC: for Bayesian PGLMM, this is the Deviance Information Criterion metric of model fit. This is set to NULL if bayes = FALSE.
REML: whether or not REML is used (TRUE or FALSE).
bayes: whether or not a Bayesian model was fit.
marginal.summ: The specified summary statistic used to summarize the Bayesian marginal distributions. Only present if bayes = TRUE
s2.init: the user-provided initial estimates of s2
B.init: the user-provided initial estimates of B
Y: the response (dependent) variable returned in matrix form
X: the predictor (independent) variables returned in matrix form (including 1s in the first column)
H: the residuals. For linear mixed models, this does not account for random terms, To get residuals after accounting for both fixed and random terms, use residuals(). For the generalized linear mixed model, these are the predicted residuals in the logit -1 space.
iV: the inverse of the covariance matrix for the entire system (of dimension (nsp * nsite) by (nsp * nsite)). This is NULL if bayes = TRUE.
mu: predicted mean values for the generalized linear mixed model (i.e., similar to fitted(merMod)). Set to NULL for linear mixed models, for which we can use fitted().
nested: matrices used to construct the nested design matrix. This is set to NULL if bayes = TRUE
Zt: the design matrix for random effects. This is set to NULL if bayes = TRUE
St: diagonal matrix that maps the random effects variances onto the design matrix
convcode: the convergence code provided by optim. This is set to NULL if bayes = TRUE
niter: number of iterations performed by optim. This is set to NULL if bayes = TRUE
inla.model: Model object fit by underlying inla function. Only returned if bayes = TRUE

Details

For Gaussian data, pglmm analyzes the phylogenetic linear mixed model

$$Y = \beta_0 + \beta_1x + b_0 + b_1x$$ $$b_0 ~ Gaussian(0, \sigma_0^2I_{sp})$$ $$b_1 ~ Gaussian(0, \sigma_0^2V_{sp})$$ $$\eta ~ Gaussian(0,\sigma^2)$$

where $\beta_0$ and $\beta_1$ are fixed effects, and $V_{sp}$ is a variance-covariance matrix derived from a phylogeny (typically under the assumption of Brownian motion evolution). Here, the variation in the mean (intercept) for each species is given by the random effect $b_0$ that is assumed to be independent among species. Variation in species' responses to predictor variable $x$ is given by a random effect $b_0$ that is assumed to depend on the phylogenetic relatedness among species given by $V_{sp}$; if species are closely related, their specific responses to $x$ will be similar. This particular model would be specified as

z <- pglmm(Y ~ X + (1|sp__), data = data, family = "gaussian", cov_ranef = list(sp = phy))

Or you can prepare the random terms manually (not recommended for simple models but may be necessary for complex models):

re.1 <- list(1, sp = dat$sp, covar = diag(nspp))

re.2 <- list(dat$X, sp = dat$sp, covar = Vsp)

z <- pglmm(Y ~ X, data = data, family = "gaussian", random.effects = list(re.1, re.2))

The covariance matrix covar is standardized to have its determinant equal to 1. This in effect standardizes the interpretation of the scalar $\sigma^2$. Although mathematically this is not required, it is a very good idea to standardize the predictor (independent) variables to have mean 0 and variance 1. This will make the function more robust and improve the interpretation of the regression coefficients. For categorical (factor) predictor variables, you will need to construct 0-1 dummy variables, and these should not be standardized (for obvious reasons).

For binary generalized linear mixed models (family = 'binomial'), the function estimates parameters for the model of the form, for example,

$$y = \beta_0 + \beta_1x + b_0 + b_1x$$ $$Y = logit^{-1}(y)$$ $$b_0 ~ Gaussian(0, \sigma_0^2I_{sp})$$ $$b_1 ~ Gaussian(0, \sigma_0^2V_{sp})$$

where $\beta_0$ and $\beta_1$ are fixed effects, and $V_{sp}$ is a variance-covariance matrix derived from a phylogeny (typically under the assumption of Brownian motion evolution).

z <- pglmm(Y ~ X + (1|sp__), data = data, family = "binomial", cov_ranef = list(sp = phy))

As with the linear mixed model, it is a very good idea to standardize the predictor (independent) variables to have mean 0 and variance 1. This will make the function more robust and improve the interpretation of the regression coefficients.

References

Ives, A. R. and M. R. Helmus. 2011. Generalized linear mixed models for phylogenetic analyses of community structure. Ecological Monographs 81:511-525.

Ives A. R. 2018. Mixed and phylogenetic models: a conceptual introduction to correlated data. https://leanpub.com/correlateddata.

Rafferty, N. E., and A. R. Ives. 2013. Phylogenetic trait-based analyses of ecological networks. Ecology 94:2321-2333.

Simpson, Daniel, et al. 2017. Penalising model component complexity: A principled, practical approach to constructing priors. Statistical science 32(1): 1-28.

Li, D., Ives, A. R., & Waller, D. M. 2017. Can functional traits account for phylogenetic signal in community composition? New Phytologist, 214(2), 607-618.

Author

Anthony R. Ives, Daijiang Li, Russell Dinnage

Examples

## Structure of examples:
# First, a (brief) description of model types, and how they are specified
# - these are *not* to be run 'as-is'; they show how models should be organised
# Second, a run-through of how to simulate, and then analyse, data
# - these *are* to be run 'as-is'; they show how to format and work with data

# \donttest{
#############################################
### Brief summary of models and their use ###
#############################################
## Model structures from Ives & Helmus (2011)
if(FALSE){
  # dat = data set for regression (note: must have a column "sp" and a column "site")
  # phy = phylogeny of class "phylo"
  # repulsion = to test phylogenetic repulsion or not
  
  # Model 1 (Eq. 1)
  z <- pglmm(freq ~ sp + (1|site) + (1|sp__@site), data = dat, family = "binomial",
             cov_ranef = list(sp = phy), REML = TRUE, verbose = TRUE, s2.init = .1)
  
  # Model 2 (Eq. 2)
  z <- pglmm(freq ~ sp + X + (1|site) + (X|sp__), data = dat, family = "binomial",
             cov_ranef = list(sp = phy), REML = TRUE, verbose = TRUE, s2.init = .1)
  
  # Model 3 (Eq. 3)
  z <- pglmm(freq ~ sp*X + (1|site) + (1|sp__@site), data = dat, family = "binomial",
             cov_ranef = list(sp = phy), REML = TRUE, verbose = TRUE, s2.init = .1)
  
  ## Model structure from Rafferty & Ives (2013) (Eq. 3)
  # dat = data set
  # phyPol = phylogeny for pollinators (pol)
  # phyPlt = phylogeny for plants (plt)
  
  z <- pglmm(freq ~ pol * X + (1|pol__) + (1|plt__) + (1|pol__@plt) +
               (1|pol@plt__) + (1|pol__@plt__),
             data = dat, family = "binomial",
             cov_ranef = list(pol = phyPol, plt = phyPlt),
             REML = TRUE, verbose = TRUE, s2.init = .1)
}

#####################################################
### Detailed analysis showing covariance matrices ###
#####################################################

# This is the example from section 4.3 in Ives, A. R. (2018) Mixed 
# and phylogenetic models: a conceptual introduction to correlated data.

library(ape)
#> Warning: package 'ape' was built under R version 4.2.2
library(mvtnorm)

# Investigating covariance matrices for different types of model structure
nspp <- 6
nsite <- 4

# Simulate a phylogeny that has a lot of phylogenetic signal (power = 1.3)
phy <- compute.brlen(rtree(n = nspp), method = "Grafen", power = 1.3)

# Simulate species means
sd.sp <- 1
mean.sp <- rTraitCont(phy, model = "BM", sigma=sd.sp^2)

# Replicate values of mean.sp over sites
Y.sp <- rep(mean.sp, times=nsite)

# Simulate site means
sd.site <- 1
mean.site <- rnorm(nsite, sd=sd.site)

# Replicate values of mean.site over sp
Y.site <- rep(mean.site, each=nspp)

# Compute a covariance matrix for phylogenetic attraction
sd.attract <- 1
Vphy <- vcv(phy)

# Standardize the phylogenetic covariance matrix to have determinant = 1. 
# (For an explanation of this standardization, see subsection 4.3.1 in Ives (2018))
Vphy <- Vphy/(det(Vphy)^(1/nspp))

# Construct the overall covariance matrix for phylogenetic attraction. 
# (For an explanation of Kronecker products, see subsection 4.3.1 in the book)
V <- kronecker(diag(nrow = nsite, ncol = nsite), Vphy)
Y.attract <- array(t(rmvnorm(n = 1, sigma = sd.attract^2*V)))

# Simulate residual errors
sd.e <- 1
Y.e <- rnorm(nspp*nsite, sd = sd.e)

# Construct the dataset
d <- data.frame(sp = rep(phy$tip.label, times = nsite), 
                site = rep(1:nsite, each = nspp))

# Simulate abundance data
d$Y <- Y.sp + Y.site + Y.attract + Y.e

# Analyze the model
pglmm(Y ~ 1 + (1|sp__) + (1|site) + (1|sp__@site), data = d, cov_ranef = list(sp = phy))
#> as(<matrix>, "dgTMatrix") is deprecated since Matrix 1.5-0; do as(as(as(., "dMatrix"), "generalMatrix"), "TsparseMatrix") instead
#> Linear mixed model fit by restricted maximum likelihood
#> 
#> Call:Y ~ 1
#> <environment: 0x000001ff72f0ee60>
#> 
#> logLik    AIC    BIC 
#> -45.13 102.26 102.46 
#> 
#> Random effects:
#>              Variance   Std.Dev
#> 1|sp        2.105e-01 0.4587672
#> 1|sp__      1.450e-07 0.0003809
#> 1|site      2.580e-06 0.0016062
#> 1|sp__@site 1.946e+00 1.3950764
#> residual    5.906e-01 0.7685368
#> 
#> Fixed effects:
#>                 Value Std.Error  Zscore Pvalue
#> (Intercept) -0.086412  0.664073 -0.1301 0.8965
#> 

# Display random effects: the function `pglmm_plot_ranef()` does what 
# the name implies. You can set `show.image = TRUE` and `show.sim.image = TRUE` 
# to see the matrices and simulations.
re <- pglmm_plot_ranef(Y ~ 1 + (1|sp__) + (1|site) + (1|sp__@site), data = d, 
                    cov_ranef = list(sp = phy), show.image = FALSE, 
                    show.sim.image = FALSE)

#################################################
### Example of a bipartite phylogenetic model ###
#################################################

# Investigating covariance matrices for different types of model structure
nspp <- 20
nsite <- 15

# Simulate a phylogeny that has a lot of phylogenetic signal (power = 1.3)
phy.sp <- compute.brlen(rtree(n = nspp), method = "Grafen", power = 1.3)
phy.site <- compute.brlen(rtree(n = nsite), method = "Grafen", power = 1.3)

# Simulate species means
mean.sp <- rTraitCont(phy.sp, model = "BM", sigma = 1)

# Replicate values of mean.sp over sites
Y.sp <- rep(mean.sp, times = nsite)

# Simulate site means
mean.site <- rTraitCont(phy.site, model = "BM", sigma = 1)

# Replicate values of mean.site over sp
Y.site <- rep(mean.site, each = nspp)

# Generate covariance matrix for phylogenetic attraction among species
sd.sp.attract <- 1
Vphy.sp <- vcv(phy.sp)
Vphy.sp <- Vphy.sp/(det(Vphy.sp)^(1/nspp))
V.sp <- kronecker(diag(nrow = nsite, ncol = nsite), Vphy.sp)
Y.sp.attract <- array(t(rmvnorm(n = 1, sigma = sd.sp.attract^2*V.sp)))
# Generate covariance matrix for phylogenetic attraction among sites
sd.site.attract <- 1
Vphy.site <- vcv(phy.site)
Vphy.site <- Vphy.site/(det(Vphy.site)^(1/nsite))
V.site <- kronecker(Vphy.site, diag(nrow = nspp, ncol = nspp))
Y.site.attract <- array(t(rmvnorm(n = 1, sigma = sd.site.attract^2*V.site)))

# Generate covariance matrix for phylogenetic attraction of species:site interaction
sd.sp.site.attract <- 1
V.sp.site <- kronecker(Vphy.site, Vphy.sp)
Y.sp.site.attract <- array(t(rmvnorm(n = 1, sigma = sd.sp.site.attract^2*V.sp.site)))

# Simulate residual error
sd.e <- 0.5
Y.e <- rnorm(nspp*nsite, sd = sd.e)

# Construct the dataset
d <- data.frame(sp = rep(phy.sp$tip.label, times = nsite), 
                site = rep(phy.site$tip.label, each = nspp))

# Simulate abundance data
d$Y <- Y.sp + Y.site + Y.sp.attract + Y.site.attract + Y.sp.site.attract + Y.e

# Plot random effects covariance matrices and then add phylogenies
# Note that, if show.image and show.sim are not specified, pglmm_plot_ranef() shows
# the covariance matrices if nspp * nsite < 200 and shows simulations 
# if nspp * nsite > 100
re <- pglmm_plot_ranef(Y ~ 1 + (1|sp__) + (1|site__) + (1|sp__@site) + 
                    (1|sp@site__) + (1|sp__@site__),
                    data=d, cov_ranef = list(sp = phy.sp, site = phy.site))


# This flips the phylogeny to match to covariance matrices
rot.phy.site <- phy.site
for(i in (nsite+1):(nsite+Nnode(phy.site))) 
   rot.phy.site <- rotate(rot.phy.site, node = i)

plot(phy.sp, main = "Species", direction = "upward")

plot(rot.phy.site, main = "Site")


# Analyze the simulated data and compute a P-value for the (1|sp__@site__) 
# random effect using a LRT. It is often better to fit the reduced model before
# the full model, because it s numerically easier to fit the reduced model, 
# and then the parameter estimates from the reduced model can be given to the
# full model. In this case, I have used the estimates of the random effects 
# from the reduce model, mod.r$ss, as the initial estimates for the same 
# parameters in the full model in the statement s2.init=c(mod.r$ss, 0.01)^2. 
# The final 0.01 is for the last random effect in the full model, (1|sp__@site__). 
# Note also that the output of the random effects from communityPGLMM(), mod.r$ss, 
# are the standard deviations, so they have to be squared for use as initial 
# values of variances in mod.f.

mod.r <- pglmm(Y ~ 1 + (1|sp__) + (1|site__) + (1|sp__@site) + (1|sp@site__), 
                        data = d, cov_ranef = list(sp = phy.sp, site = phy.site))
mod.f <- pglmm(Y ~ 1 + (1|sp__) + (1|site__) + (1|sp__@site) + (1|sp@site__) + 
               (1|sp__@site__), data = d, 
               cov_ranef = list(sp = phy.sp, site = phy.site), 
               s2.init = c(mod.r$ss, 0.01)^2)
mod.f
#> Linear mixed model fit by restricted maximum likelihood
#> 
#> Call:Y ~ 1
#> <environment: 0x000001ff72f0ee60>
#> 
#> logLik    AIC    BIC 
#> -698.6 1415.1 1438.2 
#> 
#> Random effects:
#>                Variance   Std.Dev
#> 1|sp          1.176e-05 3.429e-03
#> 1|sp__        3.044e+00 1.745e+00
#> 1|site        2.102e-04 1.450e-02
#> 1|site__      1.179e+00 1.086e+00
#> 1|sp__@site   1.954e+00 1.398e+00
#> 1|sp@site__   1.373e+00 1.172e+00
#> 1|sp__@site__ 4.392e-10 2.096e-05
#> residual      1.048e-05 3.237e-03
#> 
#> Fixed effects:
#>               Value Std.Error Zscore Pvalue
#> (Intercept) 0.75921   2.77504 0.2736 0.7844
#> 
pvalue <- pchisq(2*(mod.f$logLik - mod.r$logLik), df = 1, lower.tail = FALSE)
pvalue
#> [1] 0.9998261
# }