Model outputs
Models were run in a (state-space) hierarchical Bayesian framework, with model fitting by Markov chain Monte Carlo for 10000 iterations with a burn-in period of 2000 using the function ‘gjam’ in the package gjam (Clark et al. 2017). Model convergence was confirmed by visual assessment of the mixing of chains as well as model-fit diagnostic plots generated in the gjamPlot function of the gjam package (Fig S6 a-e). We ran models separately for each of the four treatment types (W, NW, SW, SNW), as well as the control (CTL), to compare the influence of global change drivers on growth rates and biotic interactions over time. Because the current gjamTime model does not test the influence of environmental covariates on density-dependent interactions, we ran a separate gjamTime model for each treatment type and then compared estimated species interaction matrices between models of each global change treatment versus control plots.
The effects of environmental drivers on density-independent growth rates (\(\rho\)s) were assessed via the mean and 95% Bayesian credible intervals of parameter estimates. For density-dependent interactions of species groups (αs), we calculated the difference in the mean estimates (αµ) between control plots and each global change treatment type for all species group pairs (i.e. Δαs). We then summed all changes in interspecific competition on a given species group and combined the interspecific and intraspecific Δαs to estimate the net change in competition on each species group within each treatment type. We discuss predicted steady-state distributions when one or more species groups showed non-linear patterns in equilibrium abundances over a given environmental gradient (see Clark et al. 2020).
Finally, to assess community stability, we used eigenvalue analysis from modeled interaction matrices; communities are considered stable if all real eigenvalues are negative (Allesina & Tang 2012). We also compared the rightmost (highest) real eigenvalues to compare stability across communities whereby lower (more negative) rightmost real eigenvalues denoted higher stability (Carpentier et al. 2021).