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).