R/retrospective.R

In multiSA: Multi-Stock Assessment

#' @name retrospective
#' @title Retrospective analysis
#'
#' @description Perform a retrospective analysis, successive removals of most recent years of data to evaluate
#' consistency in model estimates of biomass, recruitment, etc.
#'
#' The `summary` method returns Mohn's rho and the plot method generates a markdown report.
#' @return A `MSAretro` object containing a named lists of arrays generated by the retrospective analysis:
#' - `S_yst` Spawning output array `[y, s, t]` where `t` indexes the retrospective peel
#' - `R_yst` Recruitment array `[y, s, t]`
#' - `F_yst` Apical fishing mortality `[y, s, t]`
#' - `VB_ymft` Vulnerable biomass available to each fishery `[y, m, f, t]`
#' @param MSAassess [MSAassess-class] object
#' @param yret Vector specifying the years (positive integers and include zero) to remove for the retrospective analysis
#' @param cores Integer for the number of cores to use for parallel processing (snowfall package)
#' @export
#' @import snowfall
retrospective <- function(MSAassess, yret = 0:5, cores = 1) {

  if (cores > 1 && !snowfall::sfIsRunning()) {
    snowfall::sfInit(parallel = TRUE, cpus = cores)
    on.exit(snowfall::sfStop())
  }

  .lapply <- pbapply::pblapply
  if (snowfall::sfIsRunning()) formals(.lapply)$cl <- substitute(snowfall::sfGetCluster())

  ret <- .lapply(yret, .ret, MSAassess)

  MSAdata <- get_MSAdata(MSAassess)
  ny <- MSAdata@Dmodel@ny
  nm <- MSAdata@Dmodel@nm
  nf <- MSAdata@Dfishery@nf
  ns <- MSAdata@Dmodel@ns

  nt <- length(yret)

  F_yst <-
    R_yst <-
    S_yst <-
    log_rdev_yst <- array(NA_real_, c(ny, ns, nt))

  VB_ymft <- array(NA_real_, c(ny, nm, nf, nt))

  F_yst[] <- sapply2(ret, function(x) apply(x@report$F_yas, c(1, 3), max))
  R_yst[] <- sapply2(ret, function(x) x@report$R_ys)
  S_yst[] <- sapply2(ret, function(x) apply(x@report$S_yrs, c(1, 3), sum))
  log_rdev_yst[] <- sapply2(ret, function(x) x@obj$env$parList(x@obj$env$last.par.best)$log_rdev_ys)
  VB_ymft[] <- sapply2(ret, function(x) apply(x@report$VB_ymfrs, 1:3, sum))

  ret_out <- list(S_yst = S_yst, R_yst = R_yst, F_yst = F_yst, VB_ymft = VB_ymft,
                  log_rdev_yst = log_rdev_yst) %>%
    structure(class = "MSAretro")
  attr(ret_out, "Dlabel") <- MSAdata@Dlabel
  attr(ret_out, "yret") <- yret
  return(ret_out)
}


.ret <- function(y, MSAassess) {

  MSAdata <- get_MSAdata(MSAassess)
  if (y == MSAdata@Dmodel@nyret) return(MSAassess)

  MSAdata@Dmodel@nyret <- y

  parameters <- MSAassess@obj$env$parList(par = MSAassess@obj$par)
  map <- MSAassess@obj$env$map
  if (length(map)) {
    map_dim <- lapply(names(map), function(i) {
      dim_i <- dim(parameters[[i]])
      if (length(dim_i)) {
        array(map[[i]], dim_i)
      } else {
        map[[i]]
      }
    }) %>%
      structure(names = names(map))
  } else {
    map_dim <- list()
  }
  random <- MSAassess@obj$env$random

  data_new <- check_data(MSAdata, silent = TRUE)
  parameters_new <- make_parameters(data_new, start = parameters, map = map_dim, silent = TRUE)
  fit <- fit_MSA(data_new, parameters_new$p, parameters_new$map, random, do_sd = FALSE, silent = TRUE)

  return(fit)
}

#' @rdname retrospective
#' @param x,object Output of `retrospective` function
#' @param var Character to indicate the metric, the item in the `MSAretro` list to be plotted. See details below.
#' @param s Integer for the stock index to plot
#' @param f Integer for the fleet index to plot
#' @param ... Not used
#' @importFrom gplots rich.colors
#' @importFrom graphics legend
#' @returns
#' `plot` generic for MSAretro returns individual figures using base graphics.
#' @export
#' @method plot MSAretro
plot.MSAretro <- function(x, var = c("S_yst", "R_yst", "F_yst", "log_rdev_yst", "VB_ymft"), s = 1, f = 1, ...) {
  var <- match.arg(var)
  Dlabel <- attr(x, "Dlabel")
  peel <- attr(x, "yret")
  ny <- length(Dlabel@year)
  nm <- max(length(Dlabel@season), 1)

  year <- Dlabel@year
  ylab <- retro_label(var)

  if (grepl("yst", var)) {
    y <- yannual <- x[[var]][, s, ]
    denom <- 1
  } else if (grepl("ymft", var)) {
    yannual <- x[[var]][, , f, ]
    if (length(dim(yannual)) > 2) { # nm = 1
      y <- collapse_yearseason(yannual)
      year <- make_yearseason(year, nm)
      ylab <- paste("Seasonal", tolower(ylab))
    } else {
      y <- yannual
    }
    denom <- nm
  }

  if (all(y >= 0, na.rm = TRUE)) {
    ylim <- c(0, 1.1) * range(y, na.rm = TRUE)
  } else {
    ylim <- 1.1 * range(y, na.rm = TRUE)
  }

  rcolor <- rich.colors(ncol(y))
  for (i in 1:length(peel)) {
    if (peel[i] > 0) y[seq(ny - peel[i] + 1, ny), i] <- NA_real_
  }

  make_tinyplot(year, y, ylab, name = peel/denom, rcolor, type = "l", leg = substitute(legend(title = "Years\nremoved:")))
  points(year[ny - peel], y[cbind(ny - peel, 1:ncol(y))], pch = 16, col = rcolor)

  invisible()
}

#' @rdname retrospective
#' @param by Character indicating whether to calculate to Mohn's rho on stock or fleet-based time series
#' @return
#' `summary` generic for MSAretro returns a matrix of Mohn's rho.
#' @export
#' @method summary MSAretro
summary.MSAretro <- function(object, by = c("stock", "fleet"), ...) {
  by <- match.arg(by)
  peel <- attr(object, "yret")

  if (by == "stock") {
    ns <- dim(object$S_yst)[2]
    sname <- attr(object, "Dlabel")@stock
    if (length(sname) <= 1) sname <- "rho"

    x <- object[grepl("yst", names(object))]

    rho <- sapply(x, function(i) apply(i, 2, Mohn_rho, peel)) %>%
      t() %>%
      matrix(length(x), ns) %>%
      structure(dimnames = list(sapply(names(x), retro_label), sname))
  } else {
    fname <- attr(object, "Dlabel")@fleet
    if (length(fname) <= 1) fname <- "rho"

    x <- object[grepl("ymft", names(object))]
    ny <- dim(x[[1]])[1]
    nf <- dim(x[[1]])[3]
    nt <- dim(x[[1]])[4]

    # First season
    rho <- sapply(x, function(i) {
      array(i[, 1, , ], c(ny, nf, nt)) %>% apply(2, Mohn_rho, peel)
    }) %>%
      t() %>%
      matrix(length(x), nf) %>%
      structure(dimnames = list(sapply(names(x), retro_label), fname))
  }
  return(rho)
}

Mohn_rho <- function(x, peel) {
  ny <- dim(x)[1]
  for (i in 1:length(peel)) { # Not necessary but good practice
    if (peel[i] > 0) x[seq(ny - peel[i] + 1, ny), i] <- NA_real_
  }
  num <- x[cbind(ny - peel, 1:ncol(x))]
  denom <- x[ny - peel, peel == 0]
  rho <- num/denom - 1
  mean(rho[peel > 0])
}



#' @inheritParams report
#' @return
#' `report` generic for MSAretro invisibly returns the output of [rmarkdown::render()]: character of the path of the rendered HTML markdown report.
#' @rdname retrospective
#' @importFrom rmarkdown render
#' @importFrom utils browseURL
#' @export
#' @method report MSAretro
report.MSAretro <- function(object, filename = "retro", dir = tempdir(), open_file = TRUE, render_args = list(), ...) {

  sname <- attr(object, "Dlabel")@stock
  fname <- attr(object, "Dlabel")@fleet

  rmd <- system.file("include", "retrospective.Rmd", package = "multiSA") %>% readLines()
  rmd_split <- split(rmd, 1:length(rmd))

  fleet_ind <- grep("*ADD FLEET RMD*", rmd)
  stock_ind <- grep("*ADD STOCK RMD*", rmd)
  rmd_split[[fleet_ind]] <- make_rmd_ret_fleet(fname)
  rmd_split[[stock_ind]] <- make_rmd_ret_stock(sname)

  ####### Function arguments for rmarkdown::render
  filename_rmd <- paste0(filename, ".Rmd")

  render_args$input <- file.path(dir, filename_rmd)
  if (is.null(render_args$quiet)) render_args$quiet <- TRUE

  # Generate markdown report
  if (!dir.exists(dir)) {
    message_info("Creating directory: ", dir)
    dir.create(dir)
  }
  write(do.call(c, rmd_split), file = file.path(dir, filename_rmd))

  # Rendering markdown file
  message_info("Rendering markdown file: ", file.path(dir, filename_rmd))
  output_filename <- do.call(rmarkdown::render, render_args)
  message("Rendered file: ", output_filename)

  if (open_file) browseURL(output_filename)
  invisible(output_filename)
}

make_rmd_ret_stock <- function(sname) {
  ns <- max(length(sname), 1)

  rmd <- lapply(1:ns, function(s) {
    header <- if (!length(sname) || length(sname) == 1) NULL else paste("###", sname[s], "\n")
    out <- c("```{r}",
             paste0("plot(object, \"S_yst\", s = ", s, ")"),
             "```\n",
             "```{r}",
             paste0("plot(object, \"R_yst\", s = ", s, ")"),
             "```\n",
             "```{r}",
             paste0("plot(object, \"F_yst\", s = ", s, ")"),
             "```\n",
             "```{r}",
             paste0("plot(object, \"log_rdev_yst\", s = ", s, ")"),
             "```\n")
    c(header, out)
  })
  do.call(c, rmd)
}

make_rmd_ret_fleet <- function(fname) {
  nf <- max(length(fname), 1)

  rmd <- lapply(1:nf, function(f) {
    out <- c("```{r}",
             paste0("plot(object, \"VB_ymft\", f = ", f, ")"),
             "```\n")
  })
  do.call(c, rmd)
}

retro_label <- function(var) {
  switch(
    var,
    "F_yst" = "Fishing mortality",
    "R_yst" = "Recruitment",
    "S_yst" = "Spawning output",
    "VB_ymft" = "Vulnerable biomass",
    "log_rdev_yst" = "log Recruitment deviations"
  )
}