Generalized linear mixed models can detect unimodal species-environment relationships

Generalized linear mixed models and unimodal response

For ease of exposition we start with a logistic linear mixed model for presence–absence data as example GLMM. The same approach can be followed for count data and loglinear models, which would relate to the RC model (de Rooij, 2007). Consider the logistic linear mixed model that relates the probability of occurrence p_ij of species j in site i to a quantitative environmental variable x_i by the formula (1)${\displaystyle \text{logit}\left(p_{ij}\right)={\alpha }_j+{\beta }_j{x}_i+{\gamma }_i(i=1,\dots ,n;j=1,\dots ,m)}$ with α_j an intercept, β_j a slope and γ_i a site effect, which we all take as random parameters with normal distributions with zero mean and variances ${\sigma }_{\alpha }^2$, ${\sigma }_{\beta }^2$ and ${\sigma }_{\gamma }^2$. In this random intercept, random slope model (Gelman & Hill, 2007; Zuur et al., 2009) it is prudent to have an additional parameter ρ for the correlation between the intercepts {α_j} and slopes {β_j}; otherwise the model would change by just centering the environmental variable. Inclusion of the random site effects {γ_i} are a means to avoid pseudoreplication (Hurlbert, 1984) as they introduce correlation among species. These correlations were not modelled by Pollock, Morris & Vesk (2012) which makes their statistical tests liberal. The site effects may account for the size of the site, the fertility of the site or any other unknown factors that influence the probability of occurrence of all species in the site. The site effect γ_i will thus be expected to be related to the expected number of species in a site, that is to ∑_jp_ij and, in terms of the data, to the number of species that is observed in a site, for short the site total, defined as S_i = ∑_jy_ij. The site total and the site effect are expected to have a monotonic positive relationship.

We now turn to one of the simplest unimodal curves for presence–absence, the Gaussian logistic curve (Oksanen et al., 2001; ter Braak & Looman, 1986) (2)${\displaystyle \mathrm{logit}\left(p_{ij}\right)=a_j-\frac{{(x_i-u_j)}^2}{2t_j^2}}$ with a_j a coefficient related to maximum probability of occurrence, u_j the species optimum and t_j the tolerance of species j. This model thus has a logistic form but is nonlinear in this parameterization. By expanding the quadratic term, (3)${\displaystyle a_j-\frac{{(x_i-u_j)}^2}{2t_j^2}=a_j-\frac{1}{2t_j^2}{x}_i^2-\frac{1}{2t_j^2}{u}_j^2+\frac{1}{t_j^2}{x}_i{u}_j=\left(a_j-\frac{1}{2t_j^2}{u}_j^2\right)+\left(\frac{u_j}{t_j^2}\right)x_i-\frac{1}{2t_j^2}{x}_i^2}$ it can be fitted to the data of each individual species by a generalized linear model (GLM) by using x and x² as predictors (Jongman, ter Braak & van Tongeren, 1995; ter Braak & Looman, 1986). The fit is unimodal if the regression coefficient of x² is negative.

Might we be able to model unimodal response even without the squared term? On assuming t_j = t and setting (4)${\displaystyle {\alpha }_j=a_j-\frac{1}{2t^2}{u}_j^2,{\beta }_j=\frac{u_j}{t^2}\text{and}{\gamma }_i=\frac{-1}{2t^2}{x}_i^2,}$ we obtain Eq. (1) again. If t would vary among species then Eq. (1) does not exactly hold because $x_i^2/t_j^2$ then also depends on j. With equal tolerances, unimodal response can thus be represented by a simple linear model with site effects and, as we propose, be fitted by a GLMM based on Eq. (1). The GLMM has additional normality assumptions. In case of unimodal response, the assumption that the site effects {γ_i} are independent normal is false, as the site effects then depend on x_i through $x_i^2/t^2$. This will be the basis of our test on unimodal response in the next section. The site effects in Eq. (4) also have a nonzero mean, but that is not a problem, as the mean can be taken out and transferred to the intercepts (α_j).

The unimodal model with two or more environmental variables (ter Braak & Prentice, 1988) can similarly be rewritten as a simple linear model without squared terms if the tolerances are equal (Appendix 1). In conclusion, up to distributional assumptions, the GLMM (e.g., Eq. (1)) can be interpreted as a Gaussian logistic model with equal tolerances for the species.

A graphical tool and statistical test for unimodal response

Equation (4) suggests a graphical tool for detecting unimodal response and also a statistical test. The idea is to fit a GLMM to the binary data {y_ij} with respect to the environmental variable with values {x_i} (i = 1, …, n). In the R package lme4 (Bates, Maechler & Bolker, 2011) , the model can be fitted by ${\displaystyle \mathtt{lmer(y\; \sim \; 1\; +\; x\; +\; (1\; +\; x\mid sp)\; +\; (1\mid site),}\mathtt{family\; =\; binomial(link\; =\; "logit"),data),}}$ where y represents the vectorized response data while sp and site are factors indicating species and sites. The site effects {γ_i} obtained from the fit are then plotted against the environmental variable {x_i}. There is an indication of unimodal response in terms of the species response with respect to the environmental variable x if this graph shows a n-shaped (as opposed to u-shaped) quadratic relationship. If the shape is not quadratic but curved, a transformation of x may improve it.

In the statistical test for unimodal response, the null model is the GLMM of Eq. (1) with, specifically, independent and normally distributed site effects. The alternative model is that the site effects depend quadratically on the environmental variable. As the site effects typically also depend on the site total S, a sensitive test on unimodal response is obtained by regressing the site effects on x, x² and S, according to the model formula (5)${\displaystyle \gamma \sim x+x^2+S.}$ There is evidence of unimodal response if the squared term is significant as judged by a z-test or, equivalently, an ANOVA test on its regression coefficient, the null model being γ ∼ x + S. The R-code for making the graph and performing the test on unimodality is supplied in Supplemental Information 1.

Simulated data

In the first example series the procedure to simulate data is the following:

1. Generate n = 50 values of an environmental variable x as a random sample from the uniform distribution, x ∼ U(−2, 2).

2. Generate a vector u of length m from a uniform distribution U(−τ, τ), where τ = 2 + t, for a fixed value of t, to ensure that optima are also placed outside the sample range of x.

3. Generate a vector a of length m drawn at random from the standard normal distribution.

4. Generate binomial probabilities p_ij from the unimodal response curve (6)${\displaystyle p_{ij}={\mathrm{logit}}^{-1}\left(a_j-\frac{{(x_i-u_j)}^2}{2t_j^2}\right)}$

In the second example series, we vary the tolerance lognormally among species, t ∼ Log N(0, σ) with σ = 0.25, 0.5 and 1 so that the median tolerance is 1. For σ = 1, the coefficient of variation is larger than 100%. In the third example series we vary the number of species (m = 10, 50 and 100) with t = 1. For the rest, the simulation is set up as in the first series. We also simulated data according to Eq. (1) with the strict assumptions of the GLMM of independent normal site effects.

Each dataset was characterized by beta diversity and length of gradient. The index of beta diversity is β_w = T/S−1, where T is the total number of species, and S is the average number of species per site (Whittaker, 1960). Length of gradient, obtained by analyzing the data with detrended correspondence analysis (DCA), is expressed in standard deviation (SD) units (Hill & Gauch, 1980). Values greater than 4 are commonly taken to indicate unimodal response. Beta diversity was calculated using the asbio package (Aho, 2011) and DCA was performed in the vegan package (Oksanen et al., 2011), both in R.

Real data

We illustrate our method with three real data sets. The first is the Dune Meadow data (Jongman, ter Braak & van Tongeren, 1995). This is a small data set of 28 higher plants in 20 sites in a dune area in the Netherlands. Environmental variables, related to soil and management, were measured at each site; we use the variable Moisture.

The second data set includes the vegetation of the rising seashore on the island Skabbholmen in the Stockholm archipelago, eastern central Sweden (Cramer & Hytteborn, 1987) and is part of the Canoco package (ter Braak & Smilauer, 1998). The data set consists of 63 sites sampled in both 1978 and 1984 and contains 68 species. The environmental variable is Elevation.

The third data set involves phytoplankton communities of 203 lakes located within four climate zones and associated measurements on various environmental variables and morphological species traits of 60 species (Kruk et al., 2011). We consider the environmental variable Temperature.

For each data set we fit a GLMM according to Eq. (1) with x the noted environmental variable, plot the resulting site effects against x and test for unimodal response at the 5% significance level explained by and below Eq. (5). For the seashore data, we analyzed 1978 and 1984 separately. We also compare the regression coefficient β_j as estimated by GLMM with the optimum u_j as obtained by explicitly fitting Eq. (2) using GLM (ter Braak & Looman, 1986) for species with a well-defined optimum, that is, the squared term of which has a z-ratio smaller than −1 in the GLM model. In the small Dune Meadow data set we used z-ratio < 0.

Simulated data

Figure 1 shows the simulated response curves in example series 1; the sampled range of the environmental variable is the range of x shown. With increasing tolerance the part of the curves that is sampled shows less unimodal response. This is expressed quantitatively in the length of gradient SD units which varies between about 1 SD (not so unimodal) to 6 SD (very unimodal); the beta diversity varies correspondingly between 1 and 5.

In Fig. 2, the site effects estimated by the GLMM analysis of each of the simulated data sets are plotted against the environmental variable. In all three series, site effects shows a clear quadratic relationship with the environmental variable except for large tolerance (t = 4) in the first series. Note the decreasing range of site effects as the tolerance increases in the first row of Fig. 2; for large tolerance, the site effects are close to zero. Nevertheless, the squared term in Eq. (5) was significant in all examples (P < 0.001) so that the method detects unimodal response even if it is moderate (t = 4). When the data are simulated according to Eq. (1) with the strict normality assumptions of the GLMM, the squared term was not judged significant more often than expected on the basis of Type I error of the test.

Figure 3 shows the relationship between the random slopes (β_j) estimated by GLMM and the true optima (u_j) in the three simulation series. The relation is positive as predicted by Eq. (4). The relationship is weaker the larger the tolerance (Fig. 3, first row). With tolerance varying across species, the relationship continues to hold true surprisingly well (Fig. 3, second row), except perhaps when the coefficient of variation of the tolerance is large (>100%). The larger the number of species the clearer the predicted relationships (last row of Figs. 2 and 3).

Real data example

The site effects estimated by GLMM show a quadratic relationship with the noted environmental variable in each of the three data sets (Fig. 4 top row), with the least unimodal response in the Dune Meadow data. Unimodal response is significant in all cases, as judged by our proposed significance test (P < 0.001) and the relationship between the random slopes (β_j) estimated by GLMM and the optima (u_j) obtained from a fit of the unimodal model of Eq. (1) is close to linear (Fig. 4, bottom row). In the Dune Meadow data there is one outlier for a species with a z-ratio close to 0. In the phytoplankton data, the species with similar low values for the optimum received differential values for the slope, but otherwise there is a good agreement.