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.