Calculate measures of latent VAR community stability

Compute reactivity, return rates and contributions of interactions to stationary forecast variance from mvgam models with Vector Autoregressive dynamics

Usage

stability(object, ...)

# S3 method for mvgam
stability(object, ...)

Arguments

object: list object of class mvgam resulting from a call to mvgam() that used a Vector Autoregressive latent process model (either as VAR(cor = FALSE) or VAR(cor = TRUE))
...: ignored

Value

A data.frame containing posterior draws for each stability metric.

Details

These measures of stability can be used to assess how important inter-series dependencies are to the variability of a multivariate system and to ask how systems are expected to respond to environmental perturbations. Using the formula for a latent VAR(1) as: $$ \mu_t \sim \text{MVNormal}(A(\mu_{t - 1}), \Sigma) \quad $$ this function will calculate the long-term stationary forecast distribution of the system, which has mean $\mu_{\infty}$ and variance $\Sigma_{\infty}$, to then calculate the following quantities:

prop_int: Proportion of the volume of the stationary forecast distribution that is attributable to lagged interactions (i.e. how important are the autoregressive interaction coefficients in $A$ for explaining the shape of the stationary forecast distribution?): $$ det(A)^2 \quad $$
prop_int_adj: Same as prop_int but scaled by the number of series $p$ to facilitate direct comparisons among systems with different numbers of interacting variables: $$ det(A)^{2/p} \quad $$
prop_int_offdiag: Sensitivity of prop_int to inter-series interactions (i.e. how important are the off-diagonals of the autoregressive coefficient matrix $A$ for shaping prop_int?), calculated as the relative magnitude of the off-diagonals in the partial derivative matrix: $$ [2~det(A) (A^{-1})^T] \quad $$
prop_int_diag: Sensitivity of prop_int to intra-series interactions (i.e. how important are the diagonals of the autoregressive coefficient matrix $A$ for shaping prop_int?), calculated as the relative magnitude of the diagonals in the partial derivative matrix: $$ [2~det(A) (A^{-1})^T] \quad $$
prop_cov_offdiag: Sensitivity of $\Sigma_{\infty}$ to inter-series error correlations (i.e. how important are off-diagonal covariances in $\Sigma$ for shaping $\Sigma_{\infty}$?), calculated as the relative magnitude of the off-diagonals in the partial derivative matrix: $$ [2~det(\Sigma_{\infty}) (\Sigma_{\infty}^{-1})^T] \quad $$
prop_cov_diag: Sensitivity of $\Sigma_{\infty}$ to error variances (i.e. how important are diagonal variances in $\Sigma$ for shaping $\Sigma_{\infty}$?), calculated as the relative magnitude of the diagonals in the partial derivative matrix: $$ [2~det(\Sigma_{\infty}) (\Sigma_{\infty}^{-1})^T] \quad $$
reactivity: A measure of the degree to which the system moves away from a stable equilibrium following a perturbation. Values > 0 suggest the system is reactive, whereby a perturbation of the system in one period can be amplified in the next period. If $\sigma_{max}(A)$ is the largest singular value of $A$, then reactivity is defined as: $$ log\sigma_{max}(A) \quad $$
mean_return_rate: Asymptotic (long-term) return rate of the mean of the transition distribution to the stationary mean, calculated using the largest eigenvalue of the matrix $A$: $$ max(\lambda_{A}) \quad $$ Lower values suggest greater stability
var_return_rate: Asymptotic (long-term) return rate of the variance of the transition distribution to the stationary variance: $$ max(\lambda_{A \otimes{A}}) \quad $$ Again, lower values suggest greater stability

Major advantages of using mvgam to compute these metrics are that well-calibrated uncertainties are available and that VAR processes are forced to be stationary. These properties make it simple and insightful to calculate and inspect aspects of both long-term and short-term stability. But it is also possible to more directly inspect possible interactions among the time series in a latent VAR process. To do so, you can calculate and plot Generalized or Orthogonalized Impulse Response Functions using the irf function.

References

AR Ives, B Dennis, KL Cottingham & SR Carpenter (2003). Estimating community stability and ecological interactions from time-series data. Ecological Monographs. 73, 301-330.

Author

Nicholas J Clark

Examples

# \donttest{
# Simulate some time series that follow a latent VAR(1) process
simdat <- sim_mvgam(family = gaussian(),
                    n_series = 4,
                    trend_model = VAR(cor = TRUE),
                    prop_trend = 1)
plot_mvgam_series(data = simdat$data_train, series = 'all')


# Fit a model that uses a latent VAR(1)
mod <- mvgam(y ~ -1,
             trend_formula = ~ 1,
             trend_model = VAR(cor = TRUE),
             family = gaussian(),
             data = simdat$data_train,
             silent = 2)
#> In file included from stan/lib/stan_math/stan/math/prim/prob/von_mises_lccdf.hpp:5,
#>                  from stan/lib/stan_math/stan/math/prim/prob/von_mises_ccdf_log.hpp:4,
#>                  from stan/lib/stan_math/stan/math/prim/prob.hpp:359,
#>                  from stan/lib/stan_math/stan/math/prim.hpp:16,
#>                  from stan/lib/stan_math/stan/math/rev.hpp:16,
#>                  from stan/lib/stan_math/stan/math.hpp:19,
#>                  from stan/src/stan/model/model_header.hpp:4,
#>                  from C:/Users/uqnclar2/AppData/Local/Temp/Rtmp2bnpq5/model-3cd07f564707.hpp:2:
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp: In function 'stan::return_type_t<T_x, T_sigma, T_l> stan::math::von_mises_cdf(const T_x&, const T_mu&, const T_k&)':
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: '-Wmisleading-indentation' is disabled from this point onwards, since column-tracking was disabled due to the size of the code/headers
#>   194 |       if (cdf_n < 0.0)
#>       | 
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: adding '-flarge-source-files' will allow for more column-tracking support, at the expense of compilation time and memory

# Calulate stability metrics for this system
metrics <- stability(mod)

# Proportion of stationary forecast distribution
# attributable to lagged interactions
hist(metrics$prop_int,
     xlim = c(0, 1),
     xlab = 'Prop_int',
     main = '',
     col = '#B97C7C',
     border = 'white')


# Within this contribution of interactions, how important
# are inter-series interactions (offdiagonals of the A matrix) vs
# intra-series density dependence (diagonals of the A matrix)?
layout(matrix(1:2, nrow = 2))
hist(metrics$prop_int_offdiag,
     xlim = c(0, 1),
     xlab = '',
     main = 'Inter-series interactions',
     col = '#B97C7C',
     border = 'white')

hist(metrics$prop_int_diag,
     xlim = c(0, 1),
     xlab = 'Contribution to interaction effect',
     main = 'Intra-series interactions (density dependence)',
     col = 'darkblue',
     border = 'white')

layout(1)

# How important are inter-series error covariances
# (offdiagonals of the Sigma matrix) vs
# intra-series variances (diagonals of the Sigma matrix) for explaining
# the variance of the stationary forecast distribution?
layout(matrix(1:2, nrow = 2))
hist(metrics$prop_cov_offdiag,
     xlim = c(0, 1),
     xlab = '',
     main = 'Inter-series covariances',
     col = '#B97C7C',
     border = 'white')

hist(metrics$prop_cov_diag,
     xlim = c(0, 1),
     xlab = 'Contribution to forecast variance',
     main = 'Intra-series variances',
     col = 'darkblue',
     border = 'white')

layout(1)

# Reactivity, i.e. degree to which the system moves
# away from a stable equilibrium following a perturbation
# (values > 1 suggest a more reactive, less stable system)
hist(metrics$reactivity,
     main = '',
     xlab = 'Reactivity',
     col = '#B97C7C',
     border = 'white',
     xlim = c(-1*max(abs(metrics$reactivity)),
              max(abs(metrics$reactivity))))
abline(v = 0, lwd = 2.5)

# }