Merge pull request #73 from nicholasjclark/stability

add VAR stability metrics ala Ives et al

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: mvgam
 Title: Multivariate (Dynamic) Generalized Additive Models
 Version: 1.1.3
-Date: 2024-09-03
+Date: 2024-09-05
 Authors@R: person("Nicholas J", "Clark", , "[email protected]", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-7131-3301"))
 Description: Fit Bayesian Dynamic Generalized Additive Models to sets of time series. Users can build dynamic nonlinear State-Space models that can incorporate semiparametric effects in observation and process components, using a wide range of observation families. Estimation is performed using Markov Chain Monte Carlo with Hamiltonian Monte Carlo in the software 'Stan'. References: Clark & Wells (2022) <doi:10.1111/2041-210X.13974>.
 URL: https://github.com/nicholasjclark/mvgam, https://nicholasjclark.github.io/mvgam/

NAMESPACE

@@ -60,6 +60,7 @@ S3method(score,mvgam_forecast)
 S3method(set_coef,mvgam)
 S3method(smooth.construct,mod.smooth.spec)
 S3method(smooth.construct,moi.smooth.spec)
+S3method(stability,mvgam)
 S3method(stancode,mvgam)
 S3method(stancode,mvgam_prefit)
 S3method(standata,mvgam_prefit)
@@ -136,6 +137,7 @@ export(series_to_mvgam)
 export(set_prior)
 export(sim_mvgam)
 export(slopes)
+export(stability)
 export(stancode)
 export(standata)
 export(student)

R/stability.R

@@ -0,0 +1,227 @@
+#' Calculate measures of latent VAR community stability
+#'
+#' Compute reactivity, return rates and contributions of interactions to
+#' stationary forecast variance from
+#' \code{mvgam} models with Vector Autoregressive dynamics
+#'
+#'@name stability.mvgam
+#'@param object \code{list} object of class \code{mvgam} resulting from a call to [mvgam()]
+#'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:
+#'\deqn{
+#'\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 \eqn{\mu_{\infty}} and variance \eqn{\Sigma_{\infty}}, to then calculate the following quantities:
+#'\itemize{
+#'    \item `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 \eqn{A} for explaining the shape of the stationary forecast distribution?):
+#'     \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}:
+#'     \deqn{
+#'     det(A)^{2/p} \quad
+#'     }
+#'    \item `prop_int_offdiag`: Sensitivity of `prop_int` to intra-series
+#'    interactions (i.e. how important are the off-diagonals of the autoregressive coefficient
+#'    matrix \eqn{A} for shaping `prop_int`?), calculated as the relative magnitude of the *off-diagonals* in
+#'    the partial derivative matrix:
+#'     \deqn{
+#'     [2~det(A) (A^{-1})^T] \quad
+#'     }
+#'    \item `prop_int_diag`: Sensitivity of `prop_int` to inter-series
+#'    interactions (i.e. how important are the diagonals of the autoregressive coefficient matrix \eqn{A}
+#'    for shaping `prop_int`?), calculated as the relative magnitude of the *diagonals* in the partial derivative
+#'    matrix:
+#'     \deqn{
+#'     [2~det(A) (A^{-1})^T] \quad
+#'     }
+#'   \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 `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
+#'    \eqn{\sigma_{max}(A)} is the largest singular value of \eqn{A}, then reactivity is defined as:
+#'     \deqn{
+#'     log\sigma_{max}(A) \quad
+#'     }
+#'    \item `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 \eqn{A}:
+#'    \deqn{
+#'    max(\lambda_{A}) \quad
+#'    }
+#'    Lower values suggest greater stability
+#'    \item `var_return_rate`: Asymptotic (long-term) return rate of the variance of the transition distribution
+#'    to the stationary variance:
+#'    \deqn{
+#'    max(\lambda_{A \otimes{A}}) \quad
+#'    }
+#'    Again, lower values suggest greater stability
+#'   }
+#'@return An \code{data.frame} containing posterior draws for each of the above stability metrics.
+#'@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}}
+#' @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)
+#'
+#' # 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')
+#'
+#' # 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
+#' 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)
+#' }
+#'@export
+stability <- function(object, ...){
+  UseMethod("stability", object)
+}
+
+#'@rdname stability.mvgam
+#'@method stability mvgam
+#'@export
+stability.mvgam = function(object){
+
+  # Check trend_model
+  trend_model <- attr(object$model_data, 'trend_model')
+  if(!trend_model %in% c('VAR', 'VARcor', 'VAR1', 'VAR1cor')){
+    stop('Only VAR(1) models currently supported for calculating stability metrics',
+         call. = FALSE)
+  }
+
+  # Take posterior draws of the interaction matrix
+  B_post <- as.matrix(object, variable = 'A', regex = TRUE)
+
+  # Take posterior draws of Sigma
+  Sigma_post <- as.matrix(object, variable = 'Sigma', regex = TRUE)
+
+  # 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
+      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 = mean(abs(int_env[lower.tri(int_env)])) /
+                          (mean(abs(diag(int_env))) + mean(abs(int_env[lower.tri(int_env)]))))
+
+      # 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)
+}