R/calcSynchMetrics.R

In CamFreshwater/samSim: Closed-loop salmon dynamics simulator

#' Calculate synchrony metrics
#'
#' This function calculates rolling window estimates of several metrics related
#' to Thibaut and Connoly's (2013 Ecology Letters) analysis of aggregate
#' variability, including the synchrony index (phi), mean temporal coefficient
#' of variation among components (with weighting options), and mean pairwise
#' correlation coefficients. It is dependent on the \code{synchrony} package.
#' Contains multiple conditionals because computation can slows as the number
#' of time series, time series length, or the number of trials increase.
#'
#' @importFrom synchrony community.sync meancorr
#' @importFrom zoo rollapplyr
#'
#' @param synchList The output list \code{x_y_synchArrays.Rdata} generated by
#' \code{recoverySimulator.R} saved in the directory \code{outputs/simData/}.
#' @param windowSize A numeric representing the size of the moving window.
#' @param synch A logical (defaults to TRUE) specifying whether the synchrony
#' metric phi should be calculated.
#' @param compCV A logical (defaults to TRUE) specifying whether mean component
#' CV should be calculated.
#' @param corr A logical (defaults to FALSE) specifying whether mean pairwise
#' correlation coefficients should be calculated (approximately equivalent to
#' phi).
#' @param weight A logical specifying whether compCV should be weighted by mean
#' abundance of recruits.
#' @return Returns a list containing three metrics (synchrony, component CV,
#' and mean pairwise correlation coefficient) for four simulation output time
#' series (recruits by brood year, log(R/S), SR model residuals, and spawner
#' abundance). Each list element is a matrix with dimensions nYears x nTrials.
#'
#' @examples
#' #relatively large number of trials so will take some time to run
#' calcSynchMetrics(synchList, windowSize = 10, synch = TRUE,
#'                  compCV = TRUE, corr = FALSE, weight = TRUE)
#'
#' @export
calcSynchMetrics <- function(synchList, windowSize = 10, synch = TRUE,
                             compCV = TRUE, corr = FALSE, weight = TRUE) {
  spwn <- synchList$S
  recBY <- synchList$recBY
  nYears <- dim(recBY)[1]
  nTrials <- dim(recBY)[3]

  synchRecBY <- matrix(NA, nrow = nYears, ncol = nTrials)
  synchSpwn <- matrix(NA, nrow = nYears, ncol = nTrials)
  compCVRecBY <- matrix(NA, nrow = nYears, ncol = nTrials)
  compCVSpwn <- matrix(NA, nrow = nYears, ncol = nTrials)
  corrRecBY <- matrix(NA, nrow = nYears, ncol = nTrials)
  corrSpwn <- matrix(NA, nrow = nYears, ncol = nTrials)

  for (n in 1:nTrials) {
    if (synch == TRUE) {
      synchRecBY[ , n] <- rollapplyr(recBY[ , , n], width = windowSize,
                                     function(x) community.sync(x)$obs,
                                     fill = NA, by.column = FALSE)
      synchSpwn[ , n] <- rollapplyr(spwn[ , , n], width = windowSize,
                                    function(x) community.sync(x)$obs,
                                    fill = NA, by.column = FALSE)
    }
    if (compCV == TRUE) {
      compCVRecBY[ , n] <- rollapplyr(recBY[ , , n], width = windowSize,
                                      function(x) wtdCV(x, weight = weight),
                                      fill = NA, by.column = FALSE)
      compCVSpwn[ , n] <- rollapplyr(spwn[ , , n], width = windowSize,
                                     function(x) wtdCV(x, weight = weight),
                                     fill = NA, by.column = FALSE)
    }
    if (corr == TRUE) {
      corrRecBY[ , n] <- rollapplyr(recBY[ , , n], width = windowSize,
                                    function(x) meancorr(x)$obs,
                                    fill = NA, by.column = FALSE)
      corrSpwn[ , n] <- rollapplyr(spwn[ , , n], width = windowSize,
                                   function(x) meancorr(x)$obs,
                                   fill = NA, by.column = FALSE)
    }
  }
  outputList <- list(synchRecBY, synchSpwn, compCVRecBY, compCVSpwn, corrRecBY,
                     corrSpwn)
  names(outputList) <- c("synchRecBY", "synchSpwn", "compCVRecBY", "compCVSpwn",
                         "corrRecBY", "corrSpwn")
  return(outputList)
}