Statistical Models
We examined how nutrient addition (NPK treatment) influenced the five response variables associated with the Price equation partition, plus plot-level species richness and aboveground biomass (g/m2) using multilevel regression models. We fitted seven separate univariate multilevel regression models, one to each metric to quantify the effect of NPK treatments on local communities across time compared to community change across time in control plots.
Each univariate model included treatment (NPK or Control) as a categorical fixed effect, time since experimental start as a continuous fixed effect (in years), and their interaction. These same covariates were also allowed to vary as random intercepts and slopes among sites, blocks (nested within sites), and plots (nested within blocks).
To quantify the joint response of these metrics to NPK treatments across time, we also fitted two multivariate multilevel regression models that included multiple response variables in the same model. The first multivariate model was fitted to examine the joint response of species richness and biomass to NPK treatments; the second examined the joint response of all five components of the partition (species loss, species gain, and biomass change associated with species loss, gain and persistent species) in control and NPK plots. This multivariate approach allows for correlations between responses to be quantified. For the multivariate models assessing the joint responses between variables, we could only allow treatment, year, and their interaction to vary among sites, as models did not converge when finer grouping variables were included. We report results from the univariate models for our main results, and report the strength of the correlation between different responses estimated with the multivariate models. We visually examined plots of residuals for all models to assess whether model assumptions (e.g., homogeneity of variance) were met. Posterior predictive plots were used to visually determine how well models reproduced the data (Supplementary Information Fig S5 a-n). Our results did not qualitatively change when only sites with experiments running for varying numbers of years (all years, ≥ 3, ≥ 6, or ≥ 10 years) were included (Fig S6), and we present results using a minimum of 3 years in the main text.
For Bayesian inferences and estimates of uncertainty, all models described above were fitted using the Hamiltonian Monte Carlo (HMC) sampler Stan (Carpenter et al. 2017), and coded using the ‘brms’ package (Bürkner 2018) in the R for Statistical Computing and Graphics environment (v.4.0.2; (R Core Development Team 2019). All models were fitted with 4 chains, and varying iterations (Supplementary Information). We report the 95% Credible Intervals (hereafter CI) around the mean overall slope for each metric in the main results (Table S2, S3, S4, S5). We used weakly regularizing priors and visual inspection of HMC chains showed excellent convergence.