R/getResiduals.R

In CamFreshwater/samSim: Closed-loop salmon dynamics simulator

#' Gather stock-recruit residuals
#'
#' This function fits simple stock-recruit models (Ricker or Larkin) to
#' generate model residuals that can be merged with simulated recruitment
#' deviations. Useful for posterior predictive comparisons and generating
#' output figures if model residuals are a variable of interest.
#'
#' @param recList A list of length nCU with each element a dataframe with nrow
#' = to length of the maximum observed timeseries. Each dataframe should
#' contain the following columns: an estimate of total recruitment (by brood
#' year), spawner abundance, and brood year. Lagged spawner abundances are
#' automatically generated by the function.
#' @param modelVec A character vector of model names with length nCU
#' (\code{ricker} or \code{larkin}).
#' @return A numeric matrix of model residuals with nrow = nrow(input data)
#' and ncol = nCU.
#'
#' @examples
#' EXAMPLE DATA NEEDS TO BE GENERATED
#' @export

#This function pulls residuals from observed recruitment data
getResiduals <- function(recList, modelVec) {
  residMat <- sapply(seq_along(recList), function(h) {
    d <- recList[[h]] %>%
      dplyr::mutate(model = modelVec[h], logProd = log(totalRec/ets),
                    spwnRetYr = yr + 4, ets1 = NA, ets2 = NA, ets3 = NA,
                    resid = NA)
    for (i in 1:nrow(d)) { #add top-fitting SR model
      d$ets1[i] <- ifelse(i < 1, NA, d$ets[i - 1])
      d$ets2[i] <- ifelse(i < 2, NA, d$ets[i - 2])
      d$ets3[i] <- ifelse(i < 3, NA, d$ets[i - 3])
    }
    if (unique(d$model) == "ricker") {
      d2 <- d %>%
        dplyr::filter(!is.na(logProd))
      d2$resid = lm(logProd ~ ets, data = d2)$residuals
    } else {
      d2 <- d %>%
        dplyr::filter(!is.na(logProd), !is.na(ets3))
      d2$resid = lm(logProd ~ ets + ets1 + ets2 + ets3, data = d2)$residuals
    }
    d[d$yr %in% d2$yr, "resid"] <- d2$resid
    return(d$resid)
  })
  return(residMat)
}