R/stability.R

In mvgam: Multivariate (Dynamic) Generalized Additive Models

#' Calculate measures of latent VAR community stability
#'
#' Compute reactivity, return rates and contributions of interactions to
#' stationary forecast variance from \pkg{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 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)
#'   }
#'
#'   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:
#'     \deqn{ det(A)^2 }
#'
#'     \item `prop_int_adj`: Same as `prop_int` but scaled by the number of
#'     series \eqn{p}:
#'     \deqn{ det(A)^{2/p} }
#'
#'     \item `prop_int_offdiag`: Sensitivity of `prop_int` to inter-series
#'     interactions (off-diagonals of \eqn{A}):
#'     \deqn{ [2~det(A) (A^{-1})^T] }
#'
#'     \item `prop_int_diag`: Sensitivity of `prop_int` to intra-series
#'     interactions (diagonals of \eqn{A}):
#'     \deqn{ [2~det(A) (A^{-1})^T] }
#'
#'     \item `prop_cov_offdiag`: Sensitivity of \eqn{\Sigma_{\infty}} to
#'     inter-series error correlations:
#'     \deqn{ [2~det(\Sigma_{\infty}) (\Sigma_{\infty}^{-1})^T] }
#'
#'     \item `prop_cov_diag`: Sensitivity of \eqn{\Sigma_{\infty}} to error
#'     variances:
#'     \deqn{ [2~det(\Sigma_{\infty}) (\Sigma_{\infty}^{-1})^T] }
#'
#'     \item `reactivity`: Degree to which the system moves away from a stable
#'     equilibrium following a perturbation. If \eqn{\sigma_{max}(A)} is the
#'     largest singular value of \eqn{A}:
#'     \deqn{ \log\sigma_{max}(A) }
#'
#'     \item `mean_return_rate`: Asymptotic return rate of the mean of the
#'     transition distribution to the stationary mean:
#'     \deqn{ \max(\lambda_{A}) }
#'
#'     \item `var_return_rate`: Asymptotic return rate of the variance of the
#'     transition distribution to the stationary variance:
#'     \deqn{ \max(\lambda_{A \otimes A}) }
#'   }
#'
#'   Major advantages of using \pkg{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.
#'
#'   You can also inspect interactions among the time series in a latent VAR
#'   process using \code{\link{irf}} for impulse response functions or
#'   \code{\link{fevd}} for forecast error variance decompositions.
#'
#' @return A \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}},
#'   \code{\link{fevd}}
#'
#' @examples
#' \dontrun{
#' # 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,
#'   chains = 2,
#'   silent = 2
#' )
#'
#' # Calculate stability metrics for this system
#' metrics <- stability(mod)
#'
#' # Proportion of stationary forecast distribution attributable to interactions
#' hist(
#'   metrics$prop_int,
#'   xlim = c(0, 1),
#'   xlab = 'Prop_int',
#'   main = '',
#'   col = '#B97C7C',
#'   border = 'white'
#' )
#'
#' # Inter- vs intra-series interaction contributions
#' 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)
#'
#' # Inter- vs intra-series contributions to forecast variance
#' 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: system response to 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 <- mcmc_chains(object$model_output, 'A')

  # Take posterior draws of Sigma
  Sigma_post <- mcmc_chains(object$model_output, 'Sigma')

  # Number of series in the VAR process
  n_series <- object$n_lv

  if (is.null(n_series)) {
    n_series <- nlevels(object$obs_data$series)
  }

  metrics <- do.call(
    rbind,
    lapply(seq_len(NROW(B_post)), function(i) {
      B <- matrix(B_post[i, ], nrow = n_series, ncol = n_series, byrow = TRUE)
      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,
        byrow = TRUE
      )
      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)
}