Skip to content

Commit

Permalink
update documentation and tests for stability
Browse files Browse the repository at this point in the history
  • Loading branch information
Nicholas Clark committed Sep 9, 2024
1 parent aaa423d commit b4b5af2
Show file tree
Hide file tree
Showing 10 changed files with 158 additions and 215 deletions.
256 changes: 112 additions & 144 deletions R/stability.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,9 @@
#'that used a Vector Autoregressive latent process model (either as `VAR(cor = FALSE)` or
#'`VAR(cor = TRUE)`)
#'@param ... ignored
#'@details These measures of stability can be used to assess how systems respond to
#'environmental perturbations. Using the formula for a latent VAR(1) as:
#'@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:
#'\deqn{
#'\mu_t \sim \text{MVNormal}(A(\mu_{t - 1}), \Sigma) \quad
#'}
Expand All @@ -23,8 +24,8 @@
#' \deqn{
#' det(A)^2 \quad
#' }
#' \item `prop_int_adj`: Same as `prop_int` but scaled by the number of series to facilitate
#' direct comparisons among systems with different numbers of interacting variables \eqn{p}:
#' \item `prop_int_adj`: Same as `prop_int` but scaled by the number of series \eqn{p} to facilitate
#' direct comparisons among systems with different numbers of interacting variables:
#' \deqn{
#' det(A)^{2/p} \quad
#' }
Expand All @@ -42,20 +43,13 @@
#' \deqn{
#' [2~det(A) (A^{-1})^T] \quad
#' }
#' \item `prop_cov_offdiag`: Sensitivity of \eqn{\Sigma_{\infty}} to inter-series error correlations
#' \item `prop_cov`: Sensitivity of \eqn{\Sigma_{\infty}} to intra-series error correlations
#' (i.e. how important are off-diagonal covariances in \eqn{\Sigma} for shaping
#' \eqn{\Sigma_{\infty}}?), calculated as the relative magnitude of the *off-diagonals* in
#' the partial derivative matrix:
#' \deqn{
#' [2~det(\Sigma_{\infty}) (\Sigma_{\infty}^{-1})^T] \quad
#' }
#' \item `prop_cov_diag`: Sensitivity of \eqn{\Sigma_{\infty}} to intra-series error variances
#' (i.e. how important are diagonal variances in \eqn{\Sigma} for shaping
#' \eqn{\Sigma_{\infty}}?), calculated as the relative magnitude of the *diagonals* in
#' the partial derivative matrix:
#' \deqn{
#' [2~det(\Sigma_{\infty}) (\Sigma_{\infty}^{-1})^T] \quad
#' }
#' \item `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
Expand All @@ -77,12 +71,14 @@
#' }
#' Again, lower values suggest greater stability
#' }
#'@return A \code{data.frame} containing posterior draws for each of the above stability metrics.
#' To more directly inspect possible interactions among the time series in a latent VAR process,
#' you can inspect Generalized Orthogonalized Impulse Response Functions using the \code{\link{irf}} function.
#'@return An \code{data.frame} containing posterior draws for each stability metric.
#'@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
#'@seealso \code{\link{VAR}},\code{\link{irf}}
#'@seealso \code{\link{VAR}}, \code{\link{irf}}
#' @examples
#' \donttest{
#' # Simulate some time series that follow a latent VAR(1) process
Expand All @@ -92,12 +88,11 @@
#' prop_trend = 1)
#' plot_mvgam_series(data = simdat$data_train, series = 'all')
#'
#' # Fit a State-Space VAR(1)
#' mod <- mvgam(y ~ series,
#' trend_formula = ~ -1,
#' # Fit a model that uses a latent VAR(1)
#' mod <- mvgam(y ~ -1,
#' trend_formula = ~ 1,
#' trend_model = VAR(cor = TRUE),
#' family = gaussian(),
#' share_obs_params = TRUE,
#' data = simdat$data_train,
#' silent = 2)
#'
Expand All @@ -107,61 +102,30 @@
#' # 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)
#' xlim = c(0, 1),
#' xlab = 'Prop_int',
#' main = '',
#' col = '#B97C7C',
#' border = 'white')
#'
#' # 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)
#' # Proportion of stationary forecast distribution
#' # attributable to correlated process errors
#' hist(metrics$prop_env,
#' xlim = c(0, 1),
#' xlab = 'Prop_env',
#' main = '',
#' col = '#B97C7C',
#' border = 'white')
#'
#' # 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))))
#' main = '',
#' xlab = 'Reactivity',
#' col = '#B97C7C',
#' border = 'white',
#' xlim = c(-1*max(abs(metrics$reactivity)),
#' max(abs(metrics$reactivity))))
#' abline(v = 0, lwd = 2.5)
#' }
#'@export
Expand All @@ -172,7 +136,7 @@ stability <- function(object, ...){
#'@rdname stability.mvgam
#'@method stability mvgam
#'@export
stability.mvgam = function(object, ...){
stability.mvgam = function(object){

# Check trend_model
trend_model <- attr(object$model_data, 'trend_model')
Expand All @@ -190,82 +154,86 @@ stability.mvgam = function(object, ...){
# Number of series in the VAR process
n_series <- object$n_lv

metrics <- do.call(rbind, lapply(seq_len(NROW(B_post)), function(i){

B <- matrix(B_post[i,],
nrow = n_series,
ncol = n_series)
p <- dim(B)[1]

# If we want to get the variance of the stationary distribution (Sigma_inf)
Sigma <- matrix(Sigma_post[i,],
nrow = n_series,
ncol = n_series)
vecS_inf <- solve(diag(p * p) - kronecker(B, B)) %*% as.vector(Sigma)
Sigma_inf <- matrix(vecS_inf, nrow = p)

# The difference in volume between Sigma_inf and Sigma is:
# det(Sigma_inf - Sigma) = det(Sigma_inf) * det(B) ^ 2
# according to Ives et al 2003 (eqn 24)

# We can take partial derivatives to determine which elements of
# Sigma_inf contribute most to rates of change in the
# proportion of Sigma_inf that is due to process error
# Thanks to Mark Scheuerell for providing inspirational code
# https://github.com/mdscheuerell/safs-quant-sem-2022/blob/main/lwa_analysis.R
int_env <- det(Sigma_inf) * t(solve(Sigma_inf))

# Proportion of inter-series covariance to
# to overall environmental variation contribution (i.e. how important are
# correlated errors for controlling the shape of the stationary forecast
# distribution?)
dat <- data.frame(prop_cov_offdiag = mean(abs(int_env[lower.tri(int_env)])) /
(mean(abs(diag(int_env))) + mean(abs(int_env[lower.tri(int_env)]))))

# Proportion of error variances to stationary forecast distribution
dat$prop_cov_diag <- 1 - dat$prop_cov_offdiag

# Proportion of volume of Sigma_inf attributable to series interactions,
# measuring the degree to which interactions increase
# the variance of the stationary distribution (Sigma_inf) relative
# to the variance of the process error (Sigma)
# lower values = more stability
dat$prop_int = abs(det(B)) ^ 2

# Ives et al 2003 suggest to scale this by the number of series for more direct
# comparisons among different studies
dat$prop_int_adj <- abs(det(B)) ^ (2 / p)

# Sensitivity of the species interaction proportion to particular
# interactions is also calculated using partial derivatives
# (note the use of 2 here because we squared det(B) in the above eqn)
int_sens <- 2 * det(B) * t(solve(B))

# Proportion of interspecific contributions to
# to overall interaction contribution
dat$prop_int_offdiag <- mean(abs(int_sens[lower.tri(int_sens)])) /
(mean(abs(diag(int_sens))) + mean(abs(int_sens[lower.tri(int_sens)])))

# Proportion of density dependent contributions to
# to overall interaction contribution
dat$prop_int_diag <- 1 - dat$prop_int_offdiag

# Reactivity, measuring 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
# Following Neubert et al 2009 Ecology (Detecting reactivity)
dat$reactivity <- log(max(svd(B)$d))
if(is.null(n_series)){
n_series <- nlevels(object$obs_data$series)
}

# Return rate of transition distribution to the stationary distribution
# Asymptotic return rate of the mean
# lower values = more stability
dat$mean_return_rate <- max(abs(eigen(B)$values))
metrics <- do.call(rbind, lapply(seq_len(NROW(B_post)), function(i){

# Asymptotic return rate of the variance
# lower values = more stability
dat$var_return_rate <- max(abs(eigen(B %x% B)$values))
dat
}))
B <- matrix(B_post[i,],
nrow = n_series,
ncol = n_series)
p <- dim(B)[1]

# If we want to get the variance of the stationary distribution (Sigma_inf)
Sigma <- matrix(Sigma_post[i,],
nrow = n_series,
ncol = n_series)
vecS_inf <- solve(diag(p * p) - kronecker(B, B)) %*% as.vector(Sigma)
Sigma_inf <- matrix(vecS_inf, nrow = p)

# The difference in volume between Sigma_inf and Sigma is:
# det(Sigma_inf - Sigma) = det(Sigma_inf) * det(B) ^ 2
# according to Ives et al 2003 (eqn 24)

# We can take partial derivatives to determine which elements of
# Sigma_inf contribute most to rates of change in the
# proportion of Sigma_inf that is due to process error
# Thanks to Mark Scheuerell for providing inspirational code
# https://github.com/mdscheuerell/safs-quant-sem-2022/blob/main/lwa_analysis.R
int_env <- det(Sigma_inf) * t(solve(Sigma_inf))

# Proportion of inter-series covariance to
# to overall environmental variation contribution (i.e. how important are
# correlated errors for controlling the shape of the stationary forecast
# distribution?)
dat <- data.frame(prop_cov_offdiag = mean(abs(int_env[lower.tri(int_env)])) /
(mean(abs(diag(int_env))) + mean(abs(int_env[lower.tri(int_env)]))))

# Proportion of error variances to stationary forecast distribution
dat$prop_cov_diag <- 1 - dat$prop_cov_offdiag

# Proportion of volume of Sigma_inf attributable to series interactions,
# measuring the degree to which interactions increase
# the variance of the stationary distribution (Sigma_inf) relative
# to the variance of the process error (Sigma)
# lower values = more stability
dat$prop_int = abs(det(B)) ^ 2

# Ives et al 2003 suggest to scale this by the number of series for more direct
# comparisons among different studies
dat$prop_int_adj <- abs(det(B)) ^ (2 / p)

# Sensitivity of the species interaction proportion to particular
# interactions is also calculated using partial derivatives
# (note the use of 2 here because we squared det(B) in the above eqn)
int_sens <- 2 * det(B) * t(solve(B))

# Proportion of interspecific contributions to
# to overall interaction contribution
dat$prop_int_offdiag <- mean(abs(int_sens[lower.tri(int_sens)])) /
(mean(abs(diag(int_sens))) + mean(abs(int_sens[lower.tri(int_sens)])))

# Proportion of density dependent contributions to
# to overall interaction contribution
dat$prop_int_diag <- 1 - dat$prop_int_offdiag

# Reactivity, measuring 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
# Following Neubert et al 2009 Ecology (Detecting reactivity)
dat$reactivity <- log(max(svd(B)$d))

# Return rate of transition distribution to the stationary distribution
# Asymptotic return rate of the mean
# lower values = more stability
dat$mean_return_rate <- max(abs(eigen(B)$values))

# Asymptotic return rate of the variance
# lower values = more stability
dat$var_return_rate <- max(abs(eigen(B %x% B)$values))
dat
}))
return(metrics)
}
Binary file modified R/sysdata.rda
View file
Binary file not shown.
2 changes: 0 additions & 2 deletions man/mvgam_marginaleffects.Rd
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

Loading

0 comments on commit b4b5af2

Please sign in to comment.