Calculation of the so-called “functional stability metrics” (in the context of the estar package) relies on the interaction strengths between species composing a community. This information is classically summarized in the community matrix (sensu Novak et al. 2016, detailed in the Glossary in the main text). In community models, this is equivalent to a time-discrete version (Downing et al. 2020) of the Jacobian matrix composed of the partial derivatives of the dynamical equations describing each species population growth (Domínguez-García, Dakos, and Kéfi 2019). In real-world communities, interaction strengths can be estimated with multivariate autoregressive state-space (MARSS) models.

In this vignette, we focus on explaining how a MARSS model works, how these models can be used with the functions provided by the estar package, and the issues that may arise when empirical data is limited. To facilitate understanding, we illustrate the definition of a MARSS model in the context of this application to an aquatic community. For a more general and formal definition, refer to the MARSS package documentation.

Overview of demonstration data

Similar to what we did for the univariate stability metrics, we illustrate the use of functional stability metrics on a dataset compiled from an ecotoxicological study about the effects of insecticide on a community of aquatic macroinvertebrates (Brink et al. 1996; Wijngaarden et al. 1996). The insecticide chlorpyrifos negatively affects freshwater macroinvertebrates, and, to a lesser degree, zooplankton species. This insecticide was applied at four different concentrations (0.1, 0.9, 6, and 44 \(\mu\)g/ L), with two replicates per concentration level. Additionally, four replicates were used as control, i.e. an undisturbed systems (‘baseline’ in estar terminology). The community is composed of 128 species, classified into 5 functional groups: herbivores, detritiherbivores, carnivores, omnivores, and detritivores.

Organism mean count (log axis) of the five functional groups over the course of ecotoxicological experiment. The vertical dashed line identifies the time when chloryfiros insecticide was applied (at week 0), and facets represent the different concentrations (micro g/L).

MARSS models

Introduction to MARSS models

A complete MARSS model that describes a state process (\(\mathbf{X}\)) assessed via an observation process (\(\mathbf{Y}\)) is defined as: \[\begin{equation} \begin{gathered} \mathbf{x}_t = \mathbf{B}\mathbf{x}_{t-1}+\mathbf{u}+\mathbf{C}_{t}\mathbf{c}_{t}+\mathbf{w}_t \text{ where } \mathbf{w}_t \sim \ \text{MVN}(0,\mathbf{Q}) \\ \mathbf{y}_t = \mathbf{Z}\mathbf{x}_t+\mathbf{a}+\mathbf{D}_{t}\mathbf{d}_{t}+\mathbf{v}_t \text{ where } \mathbf{v}_t \sim \ \text{MVN}(0,\mathbf{R}) \\ \mathbf{x}_1 \sim\ \text{MVN}(\pi,\Lambda) \text{ or } \mathbf{x}_0 \sim\ \text{MVN}(\pi,\Lambda) \end{gathered} \end{equation}\]

where

\(\mathbf{X}\) is a \(m \times T\) matrix, composed of the unknown abundance values of \(m\) species or functional groups (\(fg\), rows), over time \(T\) (\(t\), columns)

\[\begin{equation} X= \begin{bmatrix} x_{fg1, t1} & x_{fg1, t2} & ... & x_{fg1, tT}\\ x_{fg2, t1} & x_{fg2, t2} & ... & x_{fg2, tT}\\ ... & ... & ... & ...\\ x_{fgm, t1} & x_{fgm, t2} & ... & x_{fgm, tT} \end{bmatrix} \end{equation}\]

\(\mathbf{Y}\) is a \(n \times T\) matrix, composed of \(n\) observations (rows) over time \(T\) (\(t\), columns)

\[\begin{equation} Y= \begin{bmatrix} y_{1, t1} & y_{1, t2} & ... & y_{1, tT}\\ y_{2, t1} & y_{2, t2} & ... & y_{2, tT}\\ ... & ... & ... & ...\\ y_{n, t1} & y_{n, t2} & ... & y_{n, tT} \end{bmatrix} \end{equation}\]

Note that:

considering that the observations in \(\mathbf{Y}\) refer to processes \(\mathbf{X}\), \(n \geq m\). Missing values can be passed as NA.
from \(\mathbf{X}\) and \(\mathbf{Y}\) alone, it is impossible to trace how observations refer to processes. We will address this in the description of \(\mathbf{Y}\).

The state process (\(\mathbf{X}\)) is described by:

\(\mathbf{x}_t\) is a \(m \times 1\) vector of unknown abundances of the \(m\) species at time \(t\).
\(\mathbf{B}\) is a \(m \times m\) matrix of states covariances. This is the matrix of interactions strengths (\(\beta\)) we want to estimate to calculate community-based stability metrics.
\(\mathbf{u}\) is a \(m \times 1\) vector of states fixed values. For example, expected abundance values for each functional group.
\(\mathbf{C}_t\) is a \(m \times p\) matrix of coefficients (\(\gamma\)) of each of the \(p\) covariates (enviromental conditions) effects on each of the \(m\) states (abundances).
\(\mathbf{c}_t\) is a \(p \times 1\) vector of covariate values. For example, environmental conditions which are expected to affect abundances.
\(\mathbf{w}_t\) is a \(m \times 1\) vector of the state (\(X\)) process errors. These account for the stochasticity of the system.

Thus, in matrix formulation:

\[\begin{equation} \begin{bmatrix} x_{1}\\ ... \\ x_{m} \end{bmatrix}_t = \begin{bmatrix} \beta_{\mathrm{1,1}} & ... & \beta_{\mathrm{1,m}} \\ ... & ... & ... \\ \beta_{\mathrm{m,1}} & ... & \beta_{\mathrm{m,m}} \end{bmatrix} \begin{bmatrix} x_{1}\\ ... \\ x_{m} \end{bmatrix}_{t-1} + \begin{bmatrix} u_{1}\\ ...\\ u_{m} \end{bmatrix} + \begin{bmatrix} \gamma_{\mathrm{1,1}} & ... & \gamma_{\mathrm{1,p}} \\ ... & ... & ... \\ \gamma_{\mathrm{m,1}} & ... & \gamma_{\mathrm{m,p}} \end{bmatrix} \begin{bmatrix} c_{\mathrm{1}} \\ ... \\ c_{\mathrm{m}} \end{bmatrix}_{t} + \begin{bmatrix} w_{1}\\ ...\\ w_{m} \end{bmatrix}_t \\ \end{equation}\]

The \(w_{t}\) errors follow a multivariate normal distribution of mean 0, and covariance \(\mathbf{Q}_t\) (a \(m \times m\) matrix). \[\begin{equation} \mathbf{w}_t \sim \,\text{MVN} \begin{pmatrix} 0, \begin{bmatrix} \sigma_{\mathrm{1}} & ... & \sigma_{\mathrm{1,m}} \\ ... & ... & ... \\ \sigma_{\mathrm{m,1}} & ... & \sigma_{\mathrm{m}} \\ \end{bmatrix} \end{pmatrix}\\ \end{equation}\]

where (\(\sigma\)) refer to variances of abundances (diagonal) and covariances among abundances of different species (off-diagonal).

The model initial conditions can be specified at \(t=0\) (\(x_0\)) or \(t=1\) (\(x_1\)).

\[\begin{equation} \begin{bmatrix} x_{fg1}\\ ...\\ x_{fgm} \end{bmatrix}_0 = \begin{bmatrix} \mu_1\\ ...\\ \mu_m \end{bmatrix}_t \end{equation}\]

The observation process (\(\mathbf{Y}\)) is described by:

\(\mathbf{y}_t\) is a \(m \times 1\) vector with the number of \(n\) observations taken at time \(t\).
\(\mathbf{Z}_t\) is a \(n \times m\) matrix which describes the relationships (\(\zeta\)) between hidden states \(\mathbf{x}\) (real abundances of functional groups) and the observations \(\mathbf{y}\) (samples).
\(\mathbf{a}\) is a \(n \times 1\) vector of fixed values.
\(\mathbf{D}_t\) is a \(n \times q\) matrix of \(\delta\) coefficients of each of the \(q\) covariates on each of the \(n\) observations.
\(\mathbf{d}_t\) s a \(q \times 1\) vector of covariate values affecting the observation process.
\(\mathbf{v}_t\) is a \(n \times 1\) vector of observation errors. They follow a multivariate normal distribution with mean 0, and covariance \(\mathbf{R}_t\) (a \(n \times n\) matrix) .

In matrix formulation:

\[\begin{equation} \begin{bmatrix} y_{1}\\ ... \\ y_{n} \end{bmatrix}_t = \begin{bmatrix} \zeta_{\mathrm{1,1}} & ... & \zeta_{\mathrm{1,m}} \\ ... & ... & ... \\ \zeta_{\mathrm{m,1}} & ... & \zeta_{\mathrm{m,m}} \end{bmatrix} \begin{bmatrix} x_{1}\\ ... \\ x_{m} \end{bmatrix}_{t-1} + \begin{bmatrix} a_{obs1}\\ ...\\ a_{obsn} \end{bmatrix} + \begin{bmatrix} \delta_{\mathrm{1,1}} & ... & \delta_{\mathrm{1,q}} \\ ... & ... & ... \\ \delta_{\mathrm{n,1}} & ... & \delta_{\mathrm{q,q}} \end{bmatrix} \begin{bmatrix} d_{\mathrm{1}} \\ ... \\ d_{\mathrm{q}} \end{bmatrix}_{t} + \begin{bmatrix} v_{obs1}\\ ...\\ v_{obsn} \end{bmatrix}_t \\ \end{equation}\]

The \(v_{t}\) errors follow a multivariate normal distribution of mean 0, and covariance \(\mathbf{R}_t\) (a \(n \times n\) matrix). \[\begin{equation} \mathbf{v}_t \sim \,\text{MVN} \begin{pmatrix} 0, \begin{bmatrix} \sigma_{\mathrm{1}} & ... & \sigma_{\mathrm{1,n}} \\ ... & ... & ... \\ \sigma_{\mathrm{n,1}} & ... & \sigma_{\mathrm{n}} \\ \end{bmatrix} \end{pmatrix}\\ \end{equation}\]

The terms \(\mathbf{y}_t\), \(\mathbf{c}_t\) and \(\mathbf{d}_t\) are empirical data that are provided by the user. The remaining terms \(\mathbf{B}\), \(\mathbf{u}\), \(\mathbf{C}\), \(\mathbf{Q}\), \(\mathbf{Z}\), \(\mathbf{a}\), \(\mathbf{D}\), and \(\mathbf{R}\) can be estimated or some of them (e.g. \(\mathbf{u}\), \(\mathbf{Q}\)) can be fixed to certain values based on the assumptions about the system. The terms \(\mathbf{c}_t\) and \(\mathbf{d}_t\) must be given as an input, without any missing values.

MARSS in estar

When calculating the functional metrics the user can either provide the required community matrix (\(\mathbf{B}\)) or have it estimated from the empirical data using a set of predefined assumptions about the modeled community. Empirical data in this case constitutes time series of population abundances (or biomass) of all the species (or functional groups) in the community. None of the parameters detailed here have default values in functions that calculate the community-based stability metrics. The motivation for it is to avoid the indiscriminate use, and thus force the user to take into consideration the implications of the chosen assumptions.

When fitting MARSS models, the user can have up to eight matrices or vectors of parameters estimated (Holmes, Scheuerell, et al., 2024). This is true for all other parameters of the model, which refer to the (co ) variances and errors of abundance and observation processes. The species interaction matrix (B) is usually fully estimated, i.e. all interaction strengths are estimated. It is however possible to constrain values to specific values or in relation to each other. We showcase the use of MARSS in a system assumed to be at equilibrium, with zero observational error, zero covariates (e..g. environmental conditions affecting abundances or the observation). These assumptions reflect knowledge about the system, which, in our case we assume to be minimal. However, for that reason, we also leave it for the user to define their own assumptions. They will affect the estimation of B. Moreover, it is important to point out that the calculation of asymptotic resilience (Downing et al., 2020) implemented in the package requires that process error and constant covariates.

Pay attention to careful specification of the parameters for the options. For example, a care should be made when specifying the initial and the final time step of the time series. If the user knows the exact timing when the disturbance took place, then the time series should start right from that time point or one time step (whatever appropriate time step is useful) after.

We present a minimal case of one single observation per state process and other simplifying assumptions.

Minimal case

The calculation of asymptotic resilience (Downing et al. (2020)) that is used in the package requires zero process error and constant covariates. Therefore, our minimal case assumes that:

Each observation corresponds to one state process (\(y_i = x_i\)), thus \(\mathbf{Z}\) = “identity”. In our case we have 5 functional groups, therefore \(\mathbf{Z}\) corresponds to a \(5 \times 5\) identity matrix, since \(n = m\):

Code

Z_I5 <- matrix(list(0), 5, 5)
diag(Z_I5) <- 1

All functional groups have the same state processes variance: \(w_t = \sigma\), thus Q = "diagonal and equal" (in MARSS terminology), which equates to \[\begin{equation} \begin{bmatrix} \sigma & 0 & 0 & 0 & 0 \\ 0 & \sigma & 0 & 0 & 0 \\ 0 & 0 & \sigma & 0 & 0 \\ 0 & 0 & 0 & \sigma & 0 \\ 0 & 0 & 0 & 0 & \sigma\\ \end{bmatrix} \end{equation}\]
Observational error is the same for all functional groups and is assumed to be 0 (we assume complete knowledge of the system): \(v_t = v\), thus R = 0
Constant covariates: \(C_{t} = C\) and \(c_{t} = c = 0\)
We assume that the system is at equilibrium, therefore: \(U\) = \(A\) = 0
\(X_0 = X_{t=1}\), which is passed to the MARSS() as the argument tinitx = 1

Such model is formalized as:

\[\begin{equation} X = \mathbf{B} \times X_{t-1} + C.0 + 0 + 0 \\ Y = Z \times X_t + 0 + D_t.0 + v \\ \end{equation}\]

Thus, our aquatic commmunity data must be split by the concentration of the pesticide, into matrices \(Y_{control}, Y_{0.1}, Y_{0.9}, Y_6, Y_{44}\), each containing the time-series of functional group counts under one treatment. We have two replicates per treatment and four replicates for the control but we are modelling one observation as corresponding to one state processes, therefore a choice must be made how to deal with this. One possibility is to randomly pick one of the replicates per treatment. Alternatively, one can choose to average the data over replicates. We here use the mean of species abundances over all replicates:

Code

## calculate z-scores of abundances
aquacommz_allscen.ldf <- aquacomm_fgps %>%
    dplyr::filter(time >= 1 , time <= 28) %>%
    ungroup() %>%
    dplyr::mutate_at(vars(herb, detr_herb, carn, omni, detr), 
                     ~MARSS::zscore(.)) %>%
    dplyr::mutate(across(c(herb, detr_herb, carn, omni, detr),
                         ~dplyr::na_if(., 0))) %>%
    tidyr::pivot_longer(cols = c(herb, detr_herb, carn, omni, detr),
                        names_to = "fgp",
                        values_to = "abund_z")

## summarize abunds over replicates
aquacommz_allscen.summldf <- aquacommz_allscen.ldf %>%
    dplyr::group_by(time, treat, fgp) %>%
    dplyr::summarize(abundz_mu = mean (abund_z),
                     abundz_sd = sd(abund_z)) %>%
    dplyr::ungroup()
#> `summarise()` has grouped output by 'time', 'treat'. You can override using the
#> `.groups` argument.

## convert into time-series matrix
aquacommz_allscen.summmxls <- aquacommz_allscen.summldf %>%
    dplyr::select(time, treat, fgp, abundz_mu) %>%
    split(.$treat) %>%
    purrr::map(~ dplyr::select(., time, fgp, abundz_mu) %>%
        unique() %>%
        tidyr::pivot_wider(id_cols = fgp,
                           names_from = time, values_from = "abundz_mu", 
                           names_prefix = "time ") %>%
        tibble::column_to_rownames(var = "fgp") %>%
        as.matrix())

Next, we prepare the arguments to be used for specifying the MARSS model according to the minimal requirements case outlined above. As we want to estimate all interaction strengths, we set B = "unconstrained".

Code

## 5 groups
### no observation error variance: 0 reps (summarized data), 5 gps
R_05 <- matrix(list(0), 5, 5)

aquacommz_allscen.marssls <- aquacommz_allscen.summmxls %>%
    purrr::map(~ MARSS(.,
                list(B = "unconstrained",
                     U = "zero", A = "zero",
                     Z = "identity", ## Z_I5 could also have been provided here 
                     Q = "diagonal and equal", 
                     R = R_05,
                     tinitx = 1),
                method = "BFGS"))
#> Success! Converged in 102 iterations.
#> Function MARSSkfas used for likelihood calculation.
#> 
#> MARSS fit is
#> Estimation method: BFGS 
#> Estimation converged in 102 iterations. 
#> Log-likelihood: 12.89424 
#> AIC: 36.21151   AICc: 188.8269   
#>  
#>                Estimate
#> B.(1,1)        -0.00384
#> B.(2,1)        -0.48776
#> B.(3,1)         0.49413
#> B.(4,1)        -0.63563
#> B.(5,1)        -0.88599
#> B.(1,2)        -0.79441
#> B.(2,2)        -0.71431
#> B.(3,2)         0.34129
#> B.(4,2)        -1.06959
#> B.(5,2)        -0.23848
#> B.(1,3)         1.06766
#> B.(2,3)         0.88814
#> B.(3,3)        -0.59104
#> B.(4,3)         1.19712
#> B.(5,3)         1.03068
#> B.(1,4)         1.06455
#> B.(2,4)         2.28849
#> B.(3,4)        -0.48379
#> B.(4,4)         1.78597
#> B.(5,4)         1.15372
#> B.(1,5)         0.64510
#> B.(2,5)         0.26483
#> B.(3,5)         0.23857
#> B.(4,5)         0.59291
#> B.(5,5)         0.76098
#> Q.diag          0.03073
#> x0.X.carn      -1.23766
#> x0.X.detr      -1.76648
#> x0.X.detr_herb -1.84701
#> x0.X.herb       0.25242
#> x0.X.omni      -0.26112
#> Initial states (x0) defined at t=1
#> 
#> Standard errors have not been calculated. 
#> Use MARSSparamCIs to compute CIs and bias estimates.
#> 
#> Success! Converged in 64 iterations.
#> Function MARSSkfas used for likelihood calculation.
#> 
#> MARSS fit is
#> Estimation method: BFGS 
#> Estimation converged in 64 iterations. 
#> Log-likelihood: 22.73044 
#> AIC: 16.53912   AICc: 169.1545   
#>  
#>                Estimate
#> B.(1,1)         -0.4165
#> B.(2,1)         -0.0461
#> B.(3,1)          0.7853
#> B.(4,1)         -0.6549
#> B.(5,1)          1.0068
#> B.(1,2)         -0.3959
#> B.(2,2)         -0.0530
#> B.(3,2)         -0.2897
#> B.(4,2)         -1.3143
#> B.(5,2)          0.4224
#> B.(1,3)         -0.1699
#> B.(2,3)          0.8898
#> B.(3,3)         -0.1784
#> B.(4,3)         -0.2465
#> B.(5,3)          0.1291
#> B.(1,4)          0.7078
#> B.(2,4)         -0.7165
#> B.(3,4)          0.2834
#> B.(4,4)          0.6668
#> B.(5,4)         -0.4993
#> B.(1,5)          0.2784
#> B.(2,5)         -0.1199
#> B.(3,5)         -0.4755
#> B.(4,5)          0.0676
#> B.(5,5)          0.4803
#> Q.diag           0.0188
#> x0.X.carn       -0.9269
#> x0.X.detr        0.7523
#> x0.X.detr_herb  -0.5069
#> x0.X.herb       -0.8082
#> x0.X.omni       -1.2319
#> Initial states (x0) defined at t=1
#> 
#> Standard errors have not been calculated. 
#> Use MARSSparamCIs to compute CIs and bias estimates.
#> 
#> Success! Converged in 457 iterations.
#> Function MARSSkfas used for likelihood calculation.
#> 
#> MARSS fit is
#> Estimation method: BFGS 
#> Estimation converged in 457 iterations. 
#> Log-likelihood: -17.78594 
#> AIC: 97.57189   AICc: 250.1873   
#>  
#>                Estimate
#> B.(1,1)         -0.6455
#> B.(2,1)         -1.6432
#> B.(3,1)         -3.2392
#> B.(4,1)          0.8320
#> B.(5,1)          2.3106
#> B.(1,2)          0.9553
#> B.(2,2)          1.1013
#> B.(3,2)          1.2936
#> B.(4,2)         -0.8920
#> B.(5,2)         -2.5817
#> B.(1,3)         -0.1544
#> B.(2,3)         -0.2534
#> B.(3,3)         -1.0449
#> B.(4,3)          0.3990
#> B.(5,3)          0.3992
#> B.(1,4)         -0.1744
#> B.(2,4)         -0.4754
#> B.(3,4)          0.0987
#> B.(4,4)          0.6302
#> B.(5,4)          1.2261
#> B.(1,5)          0.0178
#> B.(2,5)          0.2035
#> B.(3,5)          0.8882
#> B.(4,5)         -0.4878
#> B.(5,5)          0.5761
#> Q.diag           0.1425
#> x0.X.carn        0.3796
#> x0.X.detr       -1.4360
#> x0.X.detr_herb  -5.8122
#> x0.X.herb       -1.1427
#> x0.X.omni       -2.7909
#> Initial states (x0) defined at t=1
#> 
#> Standard errors have not been calculated. 
#> Use MARSSparamCIs to compute CIs and bias estimates.
#> 
#> Success! Converged in 220 iterations.
#> Function MARSSkfas used for likelihood calculation.
#> 
#> MARSS fit is
#> Estimation method: BFGS 
#> Estimation converged in 220 iterations. 
#> Log-likelihood: -23.34317 
#> AIC: 108.6863   AICc: 261.3017   
#>  
#>                Estimate
#> B.(1,1)         -2.3823
#> B.(2,1)          1.0896
#> B.(3,1)         -0.5606
#> B.(4,1)         -0.0179
#> B.(5,1)          0.7345
#> B.(1,2)         -2.6838
#> B.(2,2)          1.1122
#> B.(3,2)         -0.6071
#> B.(4,2)          0.5630
#> B.(5,2)         -0.0574
#> B.(1,3)          3.5143
#> B.(2,3)         -0.4325
#> B.(3,3)          1.2171
#> B.(4,3)         -1.1742
#> B.(5,3)         -0.7598
#> B.(1,4)         -0.0724
#> B.(2,4)          0.4573
#> B.(3,4)          0.2121
#> B.(4,4)          0.5342
#> B.(5,4)          1.4439
#> B.(1,5)          1.0358
#> B.(2,5)         -1.1006
#> B.(3,5)          0.2629
#> B.(4,5)          0.2850
#> B.(5,5)         -0.7703
#> Q.diag           0.1881
#> x0.X.carn        1.7363
#> x0.X.detr       -1.4806
#> x0.X.detr_herb  -0.4085
#> x0.X.herb       -1.2921
#> x0.X.omni        0.8345
#> Initial states (x0) defined at t=1
#> 
#> Standard errors have not been calculated. 
#> Use MARSSparamCIs to compute CIs and bias estimates.
#> 
#> Success! Converged in 187 iterations.
#> Function MARSSkfas used for likelihood calculation.
#> 
#> MARSS fit is
#> Estimation method: BFGS 
#> Estimation converged in 187 iterations. 
#> Log-likelihood: -8.653387 
#> AIC: 79.30677   AICc: 231.9222   
#>  
#>                Estimate
#> B.(1,1)         -0.0208
#> B.(2,1)          0.6099
#> B.(3,1)          0.1791
#> B.(4,1)         -0.0165
#> B.(5,1)          0.5966
#> B.(1,2)          1.1379
#> B.(2,2)          0.6065
#> B.(3,2)          2.7719
#> B.(4,2)         -0.5077
#> B.(5,2)          1.7027
#> B.(1,3)         -0.4974
#> B.(2,3)         -0.2568
#> B.(3,3)         -1.7395
#> B.(4,3)          0.8265
#> B.(5,3)         -0.5155
#> B.(1,4)          0.3383
#> B.(2,4)         -0.2229
#> B.(3,4)         -0.3574
#> B.(4,4)          0.4981
#> B.(5,4)         -0.6097
#> B.(1,5)         -1.6698
#> B.(2,5)          0.3865
#> B.(3,5)         -2.8546
#> B.(4,5)          0.5559
#> B.(5,5)         -0.9264
#> Q.diag           0.0902
#> x0.X.carn       -1.6448
#> x0.X.detr       -0.6463
#> x0.X.detr_herb   0.2271
#> x0.X.herb       -2.2936
#> x0.X.omni       -0.3627
#> Initial states (x0) defined at t=1
#> 
#> Standard errors have not been calculated. 
#> Use MARSSparamCIs to compute CIs and bias estimates.
names(aquacommz_allscen.marssls) <- paste0("Conc. = ", c("0", "0.1", "0.9", "6", "44" ), " micro g/L")

We provide the function estar::extractB() which extracts \(\mathbf{B}\) in matrix form from the fitted MARSS object and names the interacting groups

Code

(aquacomm.Bls <- aquacommz_allscen.marssls %>%
  purrr::map(~estar::extractB(., 
                states_names = c("Herbivores", "DetHerbivores", "Carnivores", "Omnivores", "Detrivores"))))

#> $`Conc. = 0 micro g/L`
#>                Herbivores DetHerbivores Carnivores  Omnivores Detrivores
#> Herbivores    -0.00384054    -0.7944114  1.0676640  1.0645532  0.6450999
#> DetHerbivores -0.48775802    -0.7143126  0.8881352  2.2884880  0.2648297
#> Carnivores     0.49412621     0.3412916 -0.5910375 -0.4837852  0.2385650
#> Omnivores     -0.63562999    -1.0695943  1.1971242  1.7859672  0.5929082
#> Detrivores    -0.88599113    -0.2384849  1.0306840  1.1537171  0.7609775
#> 
#> $`Conc. = 0.1 micro g/L`
#>                Herbivores DetHerbivores Carnivores  Omnivores  Detrivores
#> Herbivores    -0.41652765   -0.39590416 -0.1699447  0.7078464  0.27842473
#> DetHerbivores -0.04607136   -0.05295295  0.8898280 -0.7164686 -0.11986370
#> Carnivores     0.78529501   -0.28965584 -0.1784333  0.2834296 -0.47550796
#> Omnivores     -0.65492182   -1.31427312 -0.2464757  0.6667780  0.06763712
#> Detrivores     1.00676461    0.42235586  0.1290599 -0.4993199  0.48032025
#> 
#> $`Conc. = 0.9 micro g/L`
#>               Herbivores DetHerbivores Carnivores   Omnivores  Detrivores
#> Herbivores    -0.6455068     0.9553172 -0.1544395 -0.17435612  0.01778414
#> DetHerbivores -1.6432166     1.1013087 -0.2533819 -0.47536805  0.20354936
#> Carnivores    -3.2391640     1.2935895 -1.0449297  0.09869478  0.88816165
#> Omnivores      0.8320019    -0.8920486  0.3989834  0.63018015 -0.48775593
#> Detrivores     2.3105803    -2.5817354  0.3991956  1.22614773  0.57609858
#> 
#> $`Conc. = 6 micro g/L`
#>                Herbivores DetHerbivores Carnivores  Omnivores Detrivores
#> Herbivores    -2.38234463   -2.68379362  3.5142992 -0.0723508  1.0358009
#> DetHerbivores  1.08957278    1.11219717 -0.4324573  0.4572830 -1.1006100
#> Carnivores    -0.56060383   -0.60711379  1.2170955  0.2120916  0.2629435
#> Omnivores     -0.01790964    0.56304823 -1.1742214  0.5341647  0.2849972
#> Detrivores     0.73450653   -0.05738783 -0.7598349  1.4439453 -0.7703449
#> 
#> $`Conc. = 44 micro g/L`
#>                Herbivores DetHerbivores Carnivores  Omnivores Detrivores
#> Herbivores    -0.02083992     1.1379372 -0.4973880  0.3383258 -1.6698239
#> DetHerbivores  0.60990448     0.6065011 -0.2567718 -0.2228877  0.3865387
#> Carnivores     0.17913798     2.7718747 -1.7394515 -0.3573693 -2.8545996
#> Omnivores     -0.01646089    -0.5076764  0.8264983  0.4981301  0.5559434
#> Detrivores     0.59659292     1.7027159 -0.5154681 -0.6096649 -0.9263735

Finally, having estimated \(\mathbf{B}\) we can calculate all the functional stability metrics:

Reactivity

A positive value indicates the system is reactive and the higher its value, the faster the perturbation can grow. The aquatic community in our example seems to be reactive even when no insecticide is applied.

Code

aquacomm.Bls %>%
  purrr::map(estar::reactivity)

#> $`Conc. = 0 micro g/L`
#> [1] 2.561053
#> 
#> $`Conc. = 0.1 micro g/L`
#> [1] 1.46357
#> 
#> $`Conc. = 0.9 micro g/L`
#> [1] 2.772984
#> 
#> $`Conc. = 6 micro g/L`
#> [1] 2.50594
#> 
#> $`Conc. = 44 micro g/L`
#> [1] 1.569554

Maximal amplification

The largest response of the community to the perturbation, it is a measure of the transient behavior of the system.

Code

aquacomm.Bls %>%
  purrr::map(estar::max_amp)

#> $`Conc. = 0 micro g/L`
#> [1] 7.235982
#> 
#> $`Conc. = 0.1 micro g/L`
#> [1] 4.568985
#> 
#> $`Conc. = 0.9 micro g/L`
#> [1] 7.173882
#> 
#> $`Conc. = 6 micro g/L`
#> [1] 5.567664
#> 
#> $`Conc. = 44 micro g/L`
#> [1] 3.938252

Initial resilience

Initial resilience is the initial rate of return to equilibrium. The larger its value, the more stable the system, as its “worst case” initial rate of return to equilibrium is faster (Downing et al. 2020).

Code

aquacomm.Bls %>%
  purrr::map(estar::init_resil)

#> $`Conc. = 0 micro g/L`
#> [1] -0.8034616
#> 
#> $`Conc. = 0.1 micro g/L`
#> [1] -0.1462433
#> 
#> $`Conc. = 0.9 micro g/L`
#> [1] -0.5265344
#> 
#> $`Conc. = 6 micro g/L`
#> [1] -1.025015
#> 
#> $`Conc. = 44 micro g/L`
#> [1] -0.686173

Asymptotic resilience

The slowest/long-term asymptotic rate of return to equilibrium after a pulse perturbation. Asymptotic resilience can take on positive real numbers. The larger its value, the more stable the system as such communities have a faster long-term return to stable state.

Code

aquacomm.Bls %>%
  purrr::map(estar::asympt_resil)

#> $`Conc. = 0 micro g/L`
#> [1] -1.249486
#> 
#> $`Conc. = 0.1 micro g/L`
#> [1] -0.8430753
#> 
#> $`Conc. = 0.9 micro g/L`
#> [1] -1.79057
#> 
#> $`Conc. = 6 micro g/L`
#> [1] -1.169717
#> 
#> $`Conc. = 44 micro g/L`
#> [1] -1.112916

Stochastic invariability

Inverse of the maximal response variance to white noise. The larger its value, the more stable the system as it indicates less variability due to random displacements of the system from the equilibrium.

Code

aquacomm.Bls %>%
  map(estar::stoch_var)

#> $`Conc. = 0 micro g/L`
#> [1] 23.9071
#> 
#> $`Conc. = 0.1 micro g/L`
#> [1] 5.286142
#> 
#> $`Conc. = 0.9 micro g/L`
#> [1] 13.90415
#> 
#> $`Conc. = 6 micro g/L`
#> [1] 11.75702
#> 
#> $`Conc. = 44 micro g/L`
#> [1] 13.94729

Function call	Min.	Lower quart.	Mean	Median	Upper quart.	Max.	Signif.
reactivity	0.0002651	0.0002924	0.0003797	0.0003213	0.0003614	0.0013451	a
max_amp	0.0013689	0.0014863	0.0082912	0.0016132	0.0019386	0.6394164	a
init_resil	0.0002659	0.0002910	0.0003660	0.0003150	0.0003753	0.0016038	a
asympt_resil	0.0002217	0.0002457	0.0003180	0.0002694	0.0003093	0.0016270	a
stoch_var	0.0011134	0.0012181	0.0017989	0.0013416	0.0016967	0.0186624	a

Jacobian stability metrics

Overview of demonstration data

MARSS models

Introduction to MARSS models

MARSS in estar

Minimal case

Reactivity

Maximal amplification

Initial resilience

Asymptotic resilience

Stochastic invariability

Performance of functional metrics

References