Ecological Archives E093-167-A2

Claudio de Sassi, Owen T. Lewis, and Jason M. Tylianakis. 2012. Plant-mediated and nonadditive effects of two global change drivers on an insect herbivore community. Ecology 93:1892–1901.

Appendix B. Multivariate analyses and mixed effects models.

Dissimilarity measures

We conducted two sets of analyses, each based on a different dissimilarity measure. The first used only species presence–absence data, with the Jaccard dissimilarity, to focus strictly on changes in community composition. The second test used the Modified-Gower distance with base 10 (Anderson et al. 2006). This distance measure considers an order-of-magnitude change in abundance (e.g., from 0.01 to 0.1) equal to a change in composition (i.e., from 0 to 1 species), and therefore accounts for the changes in relative abundance of species in addition to changes in the community composition per se. This approach allowed us to specify explicitly the relative importance given to changes in species relative abundance vs. changes in composition in the analysis (Anderson et al. 2006).

Plant composition in the herbivore community composition analyses

To include the effect of plant composition relative to the experimental drivers, we performed a Principal Component Analysis (PCA) and extracted the scores of the four PCA axes that each explained more than 5% of the variation in plant composition (see Table B1). Together, these 4 axes explained 75.9% of the total variation in plant community composition and were included in the model as fixed effects alongside temperature and nitrogen. We tested the maximal model first, and then removed all non-significant terms until the best-fitting model was obtained.

In both plant and herbivore analyses, the effect of temperature was significant even when it entered the model after elevation, indicating that temperature had an effect on plant community structure even after controlling for other effects correlated with elevation (e.g., radiation, partial gases concentration). In contrast, elevation was not significant, even when it entered the model before temperature (P > 0.05).

TABLE B1. Factor loadings of the Principal Component Analysis on plant composition, showing the relative (% Variation) and cumulative (Cum. % Variation) contribution of the individual axes (PC) to the explained variation.





















Test for biotic homogenization

This test computes a distance matrix between the species composition of groups (in our case, between the coldest, mid, and warmest plot in each transect), and the individual distances from each site to its group centroid are used in a one-way permutational ANOVA to test for differences in multivariate dispersion between groups (Anderson et al. 2006). In other words, it tests whether, for example, the warmest sites in each transect were on average more similar to each other in their composition than were the coldest sites, or vice versa. Note that this required us to treat temperature as a categorical factor for this analysis (to have groups within which to assess similarity), rather than a variate as in all our other analyses.

Mixed effects models

To test the effect of vegetation composition on total herbivore abundance and biomass (using a Poisson error), we included the first four axes of the plant composition PCA in the initial model, and subsequently removed all non-significant scores.

When testing species richness, we included the total sample size as an additional covariate, to determine whether changes in richness were simply driven by changes in sample size.

For models using a Poisson error (abundance data), we directly tested the coefficients of our fixed effects (as recommended by Bolker et al. 2009). Due to issues associated with calculating P values from mixed effects models with a Gaussian error structure (Bolker et al. 2009), we used Markov Chain Monte Carlo (MCMC) resampling to estimate P values from Gaussian models. The MCMC procedure was carried out using the pvals.fnc function in the languageR package (Baayen 2010) for R.

TABLE B2. Comparison of the maximal models in our phenology analyses, using either temperature or elevation as predictor. Asterisks (*) between two or more predictors signify that the model included main effects of each predictor and all possible interactions between them. A colon (:) between two predictors indicates an interaction effect between them. In all cases, temperature provided a better fit than elevation. The maximal is simplified to a reduced best model if necessary. We provide a P value from a Likelihood Ratio test between the elevation and temperature model, and between the initial and the best model (where applicable).

Response variable




L.R Test

Herbivore abundance through time


elevation × nitrogen × time





temperature × nitrogen × time


P < 0.001


Best final model

temperature × nitrogen × time








Herbivore bodyweight through time


elevation × nitrogen × time





temperature × nitrogen × time


P = 0.004


Best final model

temperature × nitrogen × time








Herbivore biomass through time


elevation × nitrogen × time





temperature × nitrogen × time


P = 0.001


Best final model



P = 0.030


Anderson, M. J., K. E. Ellingsen, and B. H. McArdle. 2006. Multivariate dispersion as a measure of beta diversity. Ecology Letters 9:683–693.

Baayen, R. H. 2010. languageR: Data sets and functions with "Analyzing Linguistic Data: A practical introduction to statistics". R package version 1.0.

Bolker, B. M., M. E. Brooks, C. J. Clark, S. W. Geange, J. R. Poulsen, M. H. H. Stevens, and J. S. S. White. 2009. Generalized linear mixed models: a practical guide for ecology and evolution. Trends in Ecology & Evolution, 24:127–135.

Crawley, M. J. 2007. The R book. John Wiley and Sons, Chichester, UK.

[Back to E093-167]