A critical issue in model-based inference for studying trait-based community assembly and a solution

Statistical models for trait-environment interaction

The data required for the statistical analyses in this paper are abundance data of m species across n sites, together with trait data on the m species, and environmental data on the n sites. In its simplest form, there is one quantitative trait t (with m values {t_j}, j = 1, …, m) and one quantitative environmental variable e (with n values {e_i}, i =1, …, n). The abundance of species j in site i, denoted by y_ij, is assumed to be (similar to) a count that follows a negative-binomial distribution with mean μ_ij and unknown overdispersion. The statistical model that we use for detecting the trait-environment interaction is the GLM model (1)${\displaystyle log\left({\mu }_{ij}\right)=R_i+C_j+b_{te}{t}_j{e}_i,}$ where b_te is the coefficient measuring direction and strength of the trait-environment (t-e) interaction and R_i and C_j are the main effects for site i and species j, respectively. The main effects are thus formed by factors for site and species and can thus approximate any non-linear function of e and t, respectively. In total, the model has n + m +2 unknown parameters (including the unknown overdispersion).

An alternative for the GLM model is the GLMM model, which we present here for later reference, (2)${\displaystyle log\left({\mu }_{ij}\right)=R_i+C_j+{\beta }_j{e}_i\text{with}{\beta }_j=b_0+b_{te}{t}_j+{\epsilon }_{\beta j}}$ with β_j a species-specific slope with respect to e, modelled as a linear model of trait t, with intercept b₀ and slope b_te, and with ε_βj a normally distributed error term with mean 0 and variance ${\sigma }_{\beta }^2$. By inserting the model for β_j in the log-linear model, we see that b_te is indeed the coefficient of the interaction t_je_i. The term b₀e_i can be absorbed in the row main effect R_i, and ε_βje_i represents additional species-specific random variation that interacts with the observed environment.

Statistical tests on trait-environment interaction

The trait-environment interaction in Eq. (1) can be tested by fitting the model with and without the interaction term, the latter by setting b_te = 0, calculating the likelihood ratio (LR) statistic of the two models for the data and assessing its significance via resampling (Warton, Shipley & Hastie, 2015). We consider four ways of carrying out the test. The first test uses the LR based on a negative-binomial response distribution and is therefore rather slow. To investigate whether we could improve on speed without sacrificing the type I error rate and loosing (too much) power, we set the response distribution to Poisson in the other three tests, as any issue due to overdispersion is accounted for by the resampling procedure. Moreover, theory tells that Poisson likelihood is the only likelihood for non-negative data and models that gives consistent estimates under misspecification of the distribution; the normal likelihood/least-squares has this feature for unbounded data (Gourieroux, Monfort & Trognon, 1984a; Gourieroux, Monfort & Trognon, 1984b; Wooldridge, 1999). In detail, the four tests are:

• 1.

  anova.traitglm. The first test uses the ‘anova.traitglm’ function of version 3.11.5 of the R package mvabund (Wang et al., 2012) with site-based resampling and calculation of the LR assuming a negative-binomial distribution (as in the data-generating models of ‘Statistical models for trait-environment interaction’ and ‘Simulation models’). This function resamples sites by bootstrapping probability integral transform (PIT) residuals with all residuals across species of the same site kept together. The code for obtaining the significance of the interaction is simply:

  model1 <- traitglm(Y,E,T, composition=TRUE)

  anova.traitglm(model1, nBoot=nBoot)

  where Y, E and T contain the abundance data y_ij, the environmental values e = e_i and the trait values t =t_j, respectively. In the many-traits case, T contains the q simulated traits. The number of bootstrap samples is set to 39. With the observed sample this gives 40 samples and a minimum p-value of 1/40 = 0.025. This number is sufficient at the nominal level of 5% of the test that we used, as Monte Carlo significance tests are unbiased for any number of re-samples (Hope, 1968). We could increase the number of bootstraps to obtain a small increase in power.

• 2.

  sites. The second test also resamples sites, but differs from the first test in that it permutes sites (instead of bootstrapping them) and calculates the Poisson LR (instead of the negative-binomial LR). For statistical testing (rather than estimation), we prefer permutation to bootstrapping. Therefore, we wrote special purpose functions in R using glm on vectorised abundance data Y for calculating the LR statistic assuming Poisson distributed abundance data. When permuting sites, the values in e are permuted instead of the rows of the abundance table Y (Appendix S2). We used 39 permutations to make the test most similar to the anova.traitglm test.

• 3.

  species. This test is similar to the previous one but permutes species (instead of sites) and calculates the Poisson LR.

• 4.

  max r/c. The fourth test applies both tests 2 and 3 and takes the maximum of their p-values. It thus a GLM-based p_max test. The ‘r/c’ is a mnemonic for rows/columns (corresponding to sites/species).

Simulation models

Abundance data on m = 30 species in n = 30 sites were generated by two simulation models, a log-linear simulation model and a one-dimensional Gaussian response model, detailed in Appendix S1. In summary, two traits t and z (both m values) and two environmental variables e and x (both n values) were drawn independently from the standard normal distribution; t and e are taken as the observed trait and the observed environment, respectively, and z and x as unobserved. As such, variables z and x are latent variables that are unrelated to (i.e., independent of) the observed ones; alternatively, z and x represent simply noise, more specifically, variability among species in their response to the environment, as in GLMM models (Jamil et al., 2013; Pollock, Morris & Vesk, 2012), and unobserved variability among sites, respectively. Either way, unobserved variation is likely the case in most ecological data and needs to be considered more often in simulation studies.

The statistical test procedure seeks to detect whether the observed trait t and the observed environment e are associated (i.e., interact), without knowledge about the two latent variables z and x, as these are unobserved/unmeasured. In the null models, abundance data are generated without any interaction between the observed t and e, but with an interaction between, for example, the unobserved trait z and the observed environmental variable e.

In the Gaussian response model, this is achieved by generating the expected abundance μ_ij of species j at site i as a Gaussian response function of e and z: (3)${\displaystyle {\mu }_{ij}=h_{j}exp\left[-\frac{{\left(e_i-z_j\right)}^2}{2{\sigma }_j^2}\right],}$ where h_j is the maximum value and σ_j is the tolerance or niche breadth of species j that are both constants or random, with σ_j independent of t. As z and σ_j are independent of the observed trait t, this model by definition contains no association between t and e. This is the ‘environmentally structured but trait random’ case, in short the ‘trait random’ case. Similarly, we can define a ‘trait structured but environment random’ case, in short the ‘environment random’ case, by replacing e_i by x_i and z_j by t_j, and a ‘both random’ case by replacing e_i by x_i in Eq. (3).

The corresponding log-linear simulation model has free main effects (R_i and C_j) and one interaction, namely either z_je_i or t_jx_i respectively, e.g., for the ‘trait random’ case (4)${\displaystyle log\left({\mu }_{ij}\right)=R_i+C_j+b_{ze}{z}_j{e}_i}$ with b_ze the coefficient measuring the direction and strength of the z-e interaction, and similarly for other interactions later on. The interaction t_je_i is missing, so there is no association between t and e in the model. Equation (4) is like Eq. (3), the model of a GLMM with random species-specific slopes with respect to the environment e (Jamil et al., 2013), but without t–e interaction (b_te = 0).

The Gaussian model of Eq. (3) and log-linear model of Eq. (4) are closely related, at least when the tolerance is constant (σ_j = σ). First, the maximum h_j does not depend on the trait and thus log(h_j) can be absorbed in the free coefficient C_j in Eq. (4). Working out the square in Eq. (3) gives squared terms in e_i and z_j that can be absorbed in the free main effects of Eq. (4) and a product z_je_i with associated coefficient1∕σ², which is then equal to b_ze in Eq. (4). With variable niche breadths, the Gaussian response model is not an easy GLM.

The log-linear model allows one case that is not available in the one-dimensional Gaussian model. This case has two interaction terms: t interacts with x and e interacts with z, but t does not interact with e, (5)${\displaystyle log\left({\mu }_{ij}\right)=R_i+C_j+b_{ze}{z}_j{e}_i+b_{tx}{t}_j{x}_i.}$ These are simple models. In the reported simulations, we also included structured and unstructured noise in the expected abundances, respectively (Appendix S1); these are important to make the data more realistic, but are not essential for our main results.

Abundance data y_ij (i = 1, …, n; j = 1, …, m) were drawn from the negative-binomial distribution with mean μ_ij and overdispersion parameter 1, giving variance function ${\mu }_{ij}+{\mu }_{ij}^2$.

To mimic a situation with many traits and a single environmental variable, t_j in the above equations was taken as the sum of q independent traits, normalized to a theoretical variance of 1.

To study the power of different methods to assess parameter significance, once the type I error had been controlled, the interaction between the observed trait and environment, i.e., the term b_tet_je_i, was added to the log-linear model with non-zero regression coefficient b_te, and z_j was replaced by t_j in the Gaussian response model. Each model was simulated 1,000 times.

Predictive modelling by lasso

Brown et al. (2014) proposed selecting trait-environment interactions on the basis of their predictive power using GLM with the least absolute shrinkage and selection operator (lasso). This approach uses penalized regression and a key point is the selection of the penalty parameter; if it is set very high, no predictor variable enters the model, if it is set to 0 all variables enter. The penalty is usually selected by cross-validation. Brown et al. (2014) used a site-based cross-validation approach, treating species as fixed factors. To evaluate this approach, we simulated 1000 data sets as in the ‘trait random’ case of the Gaussian model with n = m = 30 and 10 traits and a single environmental variable. We applied the function traitglm of the R package mvabund with arguments method = “cv.glm1path” and composition = TRUE, then counted the number of simulated data sets in which the best model, i.e., that gave the “best predictive performance” (Wang et al., 2012), contained any non-zero trait-environment coefficients. As a control experiment, we analyzed data from a model without any interaction of a latent trait with the observed variables. We also applied species-based cross-validation to the same data sets to determine whether this led to similar results.

Case study

As an illustration, we re-analysed the data of Choler (2005) who investigated the shift in Alpine plant traits along a snow-melt gradient. This is the aravo data set in the R package ade4 (Dray & Dufour, 2007). It has 82 species and 75 sites. Here we considered only one trait and one environmental variable, namely Spread (maximum lateral spread of clonal plants) and Snow (mean snowmelt date in Julian day averaged over 1997–1999). Because the distribution of Spread is right-skewed, it was logarithmically transformed before analysis.

First, we performed statistical tests based on the squared fourth-corner correlation, using the R package ade4 (Dray & Dufour, 2007). Second, we performed regression-based analyses as described in Brown et al. (2014) and Warton, Shipley & Hastie (2015). In such analyses, the abundance of the species is the response variable for which we must specify a distribution. We used the default distribution of the R package mvabund (Wang et al., 2012), the negative-binomial distribution, and checked whether this distribution is tenable for these data. Following Warton, Shipley & Hastie (2015), we did this by plotting the Dunn-Smith residuals of the model ‘site+species+species:Snow’ against the fitted values. The plot (Fig. S1) does not show any clear pattern that would suggest an ill-specified mean–variance relation. We therefore continued analysis with this default distribution. We then applied the anova.traitglm test and the other three tests in ‘Statistical tests on trait-environment interaction.’ We used 999 resamples in each test. Third, as an alternative to statistical testing, we applied predictive modelling as in ‘Predictive modelling by lasso.’ With a single trait-environment interaction, the choice is between a model with and a model without the interaction.

Statistical testing results

In the Gaussian model simulations, the ‘trait random’ case is the only one that generates a significantly inflated type I error rate (>0.40) in the anova.traitglm test, whereas the max r/c test provides appropriate rates (Table 1). With six species and five or ten random traits, the type I error rate of anova.traitglm was even more inflated (>0.75), and for twenty species with 10–30 traits the type I error rate was even greater than 0.90.

In the log-linear simulation model, similar inflation of type I error rates is found. The type I error rate of the anova.traitglm test increases with the size of the coefficient b_ze in Eqs. (4) and (5) (Fig. 1, grey lines) up to values around 0.50. The sites test results in similar p-values. The species test results in much lower p-values than the sites test when b_tx = 0, but similar p-values if b_tx = b_ze (solid and dotted blue lines in Fig. 1, respectively). The latter is expected as, when b_tx = b_ze, the model is perfectly symmetric in species and sites. The max r/c test, which takes the maximum of the p-values per simulated dataset, results in p-values around or below 0.05 if b_tx = 0 but has a slightly inflated type I error rate if b_tx = b_ze > 0.4 (solid and dotted black lines in Fig. 1, respectively). The highest type I error rate of the p_max test is 0.08 (Fig. 1).

Figure 2 shows the power of max r/c test. Recall that power here is estimated at a higher than 0.05 nominal level because the test does not fully control the type I error. The actual size of the test is 0.08 (the maximum observed type I error rate for this test in Fig. 1). As expected, power increases with the trait-environment coefficient b_te, but decreases with increasing b_ze, which represents either the interaction of e with an unmeasured trait or noise in the model for the species-specific slopes in Eq. (2). The effect of an additional noise component (modelled by b_tx) on power is not very large, except near b_te = 0, where power is in fact type I error rate, which was inflated when b_ze > 0.4, as in Fig. 1.

Predictive modelling results

In the ‘trait random’ case, site-based cross-validation to select the penalty parameter of the lasso often resulted in models with trait-environment terms that have, in fact, no predictive value; the number of data sets with such terms was 851 (in 1,000 simulated data sets) with 463 data sets giving one interaction coefficient that was at least 0.1 in absolute value. In the control experiment using a model without any interaction of a latent trait with the observed variables, we found low numbers of false positive (<20 out of 1,000). With m = 10 (instead of 30), there were 855 false positives numbers with 710 and 403 data sets with a coefficient larger than 0.1 and 0.3, respectively. Applying species-based cross-validation resulted for m = 30 in two data sets with nonzero trait-environment coefficients and for m = 10 in ten such data sets.

Case study results

Statistical tests on the interaction between Spread and Snow using the fourth-corner correlation resulted in p-values of 0.002 and 0.361 for the site- and species-based permutation tests, respectively (999 permutations). In the p_max test, these p-values are combined by taking their maximum and the final p-value was thus 0.361. The conclusion from the fourth-corner analysis is that there is no evidence that Spread and Snow are associated.

Using GLM-models, the site-based approaches found evidence for interaction between Spread and Snow (p-value of 0.001 for both the anova.traitglm and sites tests). By contrast, the species and max r/c tests did not detect evidence for this interaction (p-value of 0.376 for both).

Site-based cross-validation in lasso-based predictive modelling selected the model with interaction, with a fourth-corner regression coefficient of 0.066 with respect to the standardized trait and the standardized environment variable. Species-based cross-validation selected the model without interaction.