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.