Statistical Analyses
To estimate the effect of temperature on size spectra slopes (bexponent) and standing stock community biomass, we fit separate
hierarchical models with varying intercepts across sites and year
(McElreath 2020). We chose a Bayesian approach because it easily
incorporates prior information and hierarchical model structures (Hobbs
& Hooten 2015; Dietze 2017). Because the exponent b is
continuous and can be positive or negative, we used a Gaussian
likelihood to model size spectra slopes. The predictor variable was mean
annual temperature. To account for unexplained variation among sites and
years, we included random intercepts of site and year .
Priors for the intercept were set as Normal (-1.5, 1). This puts
low prior probabilities on positive exponents and on exponents with
extreme negative values (e.g., < -4), reflecting a wide range
of possible values reported in the size spectra literature (White et al.
2007, Blanchard et al. 2009, Edwards et al. 2017). Prior distributions
for the effect of mean annual temperature were set asNormal (0,0.1), and Exponential (2) for sigma and the
standard deviation of the random intercepts. To ensure that the prior
distributions contained reasonable prior prediction but did not
overwhelm the posterior inference, we used prior predictive simulation
(see Supplemental Information, Fig. S1; Gabry et al. 2019).
For community biomass, we used a Gamma likelihood with a log link,
because biomass is a continuous and positive measure (Hobbs & Hooten
2015). The model structure for biomass was the same as the model
structure for size spectra slopes. The prior for the intercept was set
as Normal (7,1). This reflects a prior expectation that community
biomass values of ~150 to ~9000 mg dry
mass/m2 are reasonable, with a mean of
~1000. These point estimates are obtained by
exponentiating, due to the log-link, the prior mean (i.e. 7) plus and
minus two standard deviations (7-2 or 7+2). These values are compatible
with the range of values reported in the literature (Benke et al.1984; Grimm 1988; Warmbold & Wesner 2018). The remaining priors wereNormal (0,0.1) for the slope, Exponential (2) for the
standard deviation of the random intercepts, and Gamma (0.01,
0.01) for the shape parameter of the Gamma likelihood. As before, these
priors were specified using prior predictive simulation (Supplemental
Information, Fig. S1). Our prior and posterior distributions for
parameter coefficients are plotted in SI Figures S2 and S3. A
sensitivity analysis of our priors indicated our parameter coefficient
estimates were robust to halving and doubling the SD value (Fig. S6).
The full hierarchical model structure can be found in Supplemental
Information S1.
Models were specified using the brms package (Bürkner 2018) in R
(R Development Core Team 2017), with posterior distributions derived
using Hamiltonian Monte Carlo in rstan (Stan Development Team
2018). The models were run using 4 chains each with 6000 iterations, in
which the first half were discarded as warmup, resulting in 12000
posterior draws. Convergence was checked by ensuring that all r-hats
were <1.1, and by visually assessing trace plots (Gelman &
Rubin 1992). To assess model fit, we used posterior predictive checks in
which we simulated ten data sets from the posterior distribution and
graphically compared them to the original data set (Gelman et al. 2013;
Gabry et al. 2019; Figure S4 and S5). Strong visual discrepancies
between the original data and simulated data indicate poor model fit
(Gabry et al. 2019). In our simulations, both models produced data that
were qualitatively similar to the original data, indicating good model
fit (SI, Figures S4 & S5).