Model specifications
Joint responses of species groups to environmental predictors, interactions among species groups, and the combination of these processes were estimated using the gjamTime model as described in Clarket al (2020) via the package gjam (Clark et al. 2017) with the gjamTime supplemental functions https://github.com/jimclarkatduke/gjam/blob/master/gjamTimeFunctions.R?raw=True.
Relative abundances of species groups were modeled as raw counts (‘hits’) of the a priori dominant, subdominant, moderate and rare species groups in each plot-year combination censored by the total number of vegetative hits within the same plot-year with a maximum of 100 hits per plot (see Plant community surveys). This censoring value reflects the observation effort term as described in Clark et al (2020) and we used the ‘DA’ (discrete abundance) data type specification for count data. Censored response data are then stored as a latent vector (ws) with a joint multivariate normal distribution with a mean of μs, which is a length s mean vector, and an error Σ, which is an s x s covariance matrix. In other words:
\(w_{s}\sim MVN(u\)s, Σ) (eq. 1).
Changes in population density of each species group over time is modeled using a Lotka-Volterra (LV) model specification from which the gjamTime model is derived:
\(\frac{dw_{s}}{\text{dt}}\) =\({(w}_{s}\ x\ X)\rho_{s}+\ {(w}_{s}\text{\ x\ }w_{s^{\prime}})\alpha_{s}+\ \varepsilon_{s}\)(eq. 2).
The first term defines the density-independent growth rate of a species group (ρs.) multiplied by the density of species groups and the environmental impact (ws × X). The second term defines the species-group’s density-dependent growth rate αs, which is modified by the density of two interacting species-groups s and s ’ (ws × ws’). Finally, the last term encompasses residual species group error (εs) (eq. 2).
Because this is a community of functionally similar herbaceous plant species competing for limited resources during a short growing season, we set model priors for α parameters to allow for negative (-1, 0) species group interactions (i.e. competition) only. For ρ intercepts, we set wide model priors from (-1, 1) to allow for species groups to increase or decrease by a maximum of 100% of their cover in a given time step (1 year). We set priors on ρ coefficients as (-0.5, 0.5) to allow a 50% change (positive or negative) in ρ in response to a given environmental driver at each time step.
Equation 2 is then reorganized as the discrete-time version of the LV model (eq 2.1) for model fitting:
\(\Delta w_{\text{st}}=P^{\prime}v_{\text{st}}+A\ u_{\text{st}}\ +\Sigma^{1/2}\varepsilon_{\text{st}}\)(eq. 2.1)
where \(\Delta w_{\text{st}}\) is the growth increment for population abundances of s species group, P and A are sparse matrices which reorganize \(\rho\) and α coefficients respectively to optimize posterior simulation and allow for direct sampling (see Clark et al 2020 SI Appendix S2.9, 2.10), \(v_{\text{st}}\) is a length-V vector where V is a block matrix of all possible combinations of species abundances\(w_{\text{st}}\) and q environmental variables (\(w_{\text{st}}\text{x\ }X_{\text{qt}}\)), \(u_{\text{st}}\) is a length-U vector where U is a block matrix of all possible combinations of species group interactions (\(w_{\text{st}}\text{\ x\ }w_{s^{\prime}t}),\ \text{\ ε}_{\text{st}}\) is a random vector and \(\Sigma^{1/2}\) is a square root matrix for thes x s process error covariance.
Finally, steady-state abundance distributions, i.e. probabilistic predicted equilibrium abundances of species groups, were estimated by numerical integration of the modeled parameter estimates of environmental effects on growth rates and interactions among species groups allowing for interactive and non-linear responses to emerge across environmental gradients (i.e. Environment x Species interactions- ESIs, Clark et al 2020). For each model output, we simulated 100 equilibrium abundance values for each species group at 10^3 discrete steps (10 steps for each covariate x covariate combination) across observed gradients of snow depth, N deposition and temperature, calculating a mean and standard deviation of the \(w_{s}^{*}\) estimates for each set of 100 simulations (See supplementary methods).
We chose to run models at the species group level (dominant, subdominant, moderate and rare) rather than at the individual species level to improve model predictions and to address broad questions about shifts in community structure under global change. This species-group approach captures rare species’ joint contribution to community dynamics in a biologically meaningful way, and allows for conceptual comparisons across multiple global change scenarios that would be too complex using the entire species set (n=20). The species-grouping model also better predicts the data (Fig S1), in particular for species with moderate and rare coverage, thus increasing parameter estimate confidence.