R/SSplotRecdevs.R

In r4ss/r4ss: R Code for Stock Synthesis

#' Plot recruitment deviations
#'
#' Plot recruitment deviations and associated quantities including derived
#' measures related to bias adjustment.
#'
#'
#' @template replist
#' @param subplots vector controlling which subplots to create
#' @template plot
#' @template print
#' @param add add to existing plot (not yet implemented)
#' @param uncertainty include plots showing uncertainty?
#' @param minyr optional input for minimum year to show in plots
#' @param maxyr optional input for maximum year to show in plots
#' @param forecastplot include points from forecast years?
#' @param col1 first color used
#' @param col2 second color used
#' @param col3 third color used
#' @param col4 fourth color used
#' @template legendloc
#' @template labels
#' @template pwidth
#' @template pheight
#' @template punits
#' @template res
#' @template ptsize
#' @template cex.main
#' @template plotdir
#' @template verbose
#' @author Ian Taylor, Ian Stewart
#' @export
#' @seealso [SS_plots()], [SS_fitbiasramp()]
SSplotRecdevs <-
  function(
    replist,
    subplots = 1:3,
    plot = TRUE,
    print = FALSE,
    add = FALSE,
    uncertainty = TRUE,
    minyr = -Inf,
    maxyr = Inf,
    forecastplot = FALSE,
    col1 = "black",
    col2 = "blue",
    col3 = "green3",
    col4 = "red",
    legendloc = "topleft",
    labels = c(
      "Year", # 1
      "Asymptotic standard error estimate", # 2
      "Log recruitment deviation", # 3
      "Bias adjustment fraction, 1 - stddev^2 / sigmaR^2"
    ), # 4
    pwidth = 6.5,
    pheight = 5.0,
    punits = "in",
    res = 300,
    ptsize = 10,
    cex.main = 1,
    plotdir = "default",
    verbose = TRUE
  ) {
    # Plot of recrecruitment deviations,  asymptotic error check, and bias adjustment

    # table to store information on each plot
    plotinfo <- NULL
    if (plotdir == "default") {
      plotdir <- replist[["inputs"]][["dir"]]
    }

    parameters <- replist[["parameters"]]
    recruit <- replist[["recruit"]]
    startyr <- replist[["startyr"]]
    endyr <- replist[["endyr"]]
    sigma_R_in <- replist[["sigma_R_in"]]

    recdevEarly <- parameters[
      substring(parameters[["Label"]], 1, 13) %in% c("Early_RecrDev"),
    ]
    early_initage <- parameters[
      substring(parameters[["Label"]], 1, 13) %in% c("Early_InitAge"),
    ]
    main_initage <- parameters[
      substring(parameters[["Label"]], 1, 12) %in% c("Main_InitAge"),
    ]
    recdev <- parameters[
      substring(parameters[["Label"]], 1, 12) %in% c("Main_RecrDev"),
    ]
    recdevFore <- parameters[
      substring(parameters[["Label"]], 1, 8) == "ForeRecr",
    ]
    recdevLate <- parameters[
      substring(parameters[["Label"]], 1, 12) == "Late_RecrDev",
    ]

    if (nrow(recdev) == 0 || max(recdev[["Value"]]) == 0) {
      if (verbose) message("Skipped SSplotrecdevs - no rec devs estimated")
    } else {
      if (nrow(recdev) > 0) {
        # early
        recdev[["Yr"]] <- as.numeric(substring(recdev[["Label"]], 14))
        if (nrow(recdevEarly) > 0) {
          recdevEarly[["Yr"]] <- as.numeric(substring(
            recdevEarly[["Label"]],
            15
          ))
        } else {
          recdevEarly[["Yr"]] <- integer(0) # empty value to add column to data.frame with 0 rows
        }
        if (nrow(early_initage) > 0) {
          early_initage[["Yr"]] <- startyr -
            as.numeric(substring(early_initage[["Label"]], 15))
          recdevEarly <- rbind(early_initage, recdevEarly)
        }
        # main
        if (nrow(main_initage) > 0) {
          main_initage[["Yr"]] <- startyr -
            as.numeric(substring(main_initage[["Label"]], 14))
          recdev <- rbind(main_initage, recdev)
        }
        # forecast
        if (nrow(recdevFore) > 0) {
          recdevFore[["Yr"]] <- as.numeric(substring(recdevFore[["Label"]], 10))
        } else {
          recdevFore[["Yr"]] <- NULL
        }
        if (nrow(recdevLate) > 0) {
          recdevLate[["Yr"]] <- as.numeric(substring(recdevLate[["Label"]], 14))
          recdevFore <- rbind(recdevLate, recdevFore)
        }

        Yr <- c(recdevEarly[["Yr"]], recdev[["Yr"]], recdevFore[["Yr"]])
        if (forecastplot) {
          goodyrs <- ifelse(Yr >= minyr & Yr <= maxyr, TRUE, FALSE)
        } else {
          # TRUE/FALSE of in range or not
          goodyrs <- Yr <= endyr + 1 & Yr >= minyr & Yr <= maxyr
        }
        xlim <- range(Yr[goodyrs], na.rm = TRUE)
        ylim <- range(
          c(recdevEarly[["Value"]], recdev[["Value"]], recdevFore[["Value"]])[
            goodyrs
          ],
          na.rm = TRUE
        )

        recdevfunc <- function(uncertainty) {
          # recdevs with uncertainty intervals
          alldevs <- rbind(recdevEarly, recdev, recdevFore)[goodyrs, ]

          colvec <- c(
            rep(col2, nrow(recdevEarly)),
            rep(col1, nrow(recdev)),
            rep(col2, nrow(recdevFore))
          )[goodyrs]
          ## alldevs[["Parm_StDev"]][is.na(alldevs[["Parm_StDev"]])] <- 0
          val <- alldevs[["Value"]]
          Yr <- alldevs[["Yr"]]
          if (uncertainty) {
            std <- alldevs[["Parm_StDev"]]
            recdev_hi <- val + 1.96 * std
            recdev_lo <- val - 1.96 * std
            ylim <- range(recdev_hi, recdev_lo, na.rm = TRUE)
          } else {
            ylim <- range(val, na.rm = TRUE)
          }
          plot(
            Yr,
            Yr,
            type = "n",
            xlab = labels[1],
            ylab = labels[3],
            ylim = ylim
          )
          abline(h = 0, col = "grey")
          if (uncertainty) {
            arrows(
              Yr,
              recdev_lo,
              Yr,
              recdev_hi,
              length = 0.03,
              code = 3,
              angle = 90,
              lwd = 1.2,
              col = colvec
            )
          }
          lines(Yr, val, lty = 3)
          points(Yr, val, pch = 16, col = colvec)
        }

        # the following code only applies when uncertainty was computed
        if (uncertainty) {
          recdevfunc3 <- function() {
            # std. dev. of recdevs
            par(mar = par("mar")[c(1:3, 2)])
            ymax <- 1.1 *
              max(
                recdev[["Parm_StDev"]],
                recdevEarly[["Parm_StDev"]],
                recdevFore[["Parm_StDev"]],
                sigma_R_in,
                na.rm = TRUE
              )
            plot(
              recdev[["Yr"]],
              recdev[["Parm_StDev"]],
              xlab = labels[1],
              main = "Recruitment deviation variance",
              cex.main = cex.main,
              ylab = labels[2],
              xlim = xlim,
              ylim = c(0, ymax),
              type = "b"
            )
            if (nrow(recdevEarly) > 0) {
              lines(
                recdevEarly[["Yr"]],
                recdevEarly[["Parm_StDev"]],
                type = "b",
                col = col2
              )
            }
            if (forecastplot & nrow(recdevFore) > 0) {
              lines(
                recdevFore[["Yr"]],
                recdevFore[["Parm_StDev"]],
                type = "b",
                col = col2
              )
            }
            abline(h = 0, col = "grey")
            abline(h = sigma_R_in, col = col4)

            ## # bias correction (2nd axis, scaled by ymax)
            ## lines(recruit[["year"]],ymax*recruit[["biasadj"]],col=col3)
            ## abline(h=ymax*1,col=col3,lty=3)
            ## ypts <- pretty(0:1)
            ## axis(side=4,at=ymax*ypts,label=ypts)
            ## mtext("Bias adjustment fraction",side=4,line=3,cex=par()[["cex"]])
          }
        } # end if uncertainty==TRUE

        if (plot) {
          # if plotting to screen or PDF
          if (1 %in% subplots) {
            recdevfunc(uncertainty = FALSE)
          }
          if (uncertainty) {
            if (2 %in% subplots) {
              recdevfunc(uncertainty = TRUE)
            }
            if (3 %in% subplots) recdevfunc3()
          }
        }
        if (print) {
          # if printing to PNG files
          if (1 %in% subplots) {
            file <- "recdevs1_points.png"
            caption <- "Recruitment deviations"
            plotinfo <- save_png(
              plotinfo = plotinfo,
              file = file,
              plotdir = plotdir,
              pwidth = pwidth,
              pheight = pheight,
              punits = punits,
              res = res,
              ptsize = ptsize,
              caption = caption
            )
            recdevfunc(uncertainty = FALSE)
            dev.off()
          }
          if (uncertainty) {
            if (2 %in% subplots) {
              file <- "recdevs2_withbars.png"
              caption <- "Recruitment deviations with 95% intervals"
              plotinfo <- save_png(
                plotinfo = plotinfo,
                file = file,
                plotdir = plotdir,
                pwidth = pwidth,
                pheight = pheight,
                punits = punits,
                res = res,
                ptsize = ptsize,
                caption = caption
              )
              recdevfunc(uncertainty = TRUE)
              dev.off()
            }
            if (3 %in% subplots) {
              file <- "recdevs3_varcheck.png"
              caption <-
                paste(
                  "Recruitment deviations variance check.<br>",
                  "See later figure of transformed variance values for comparison",
                  "with bias adjustment settings in the model."
                )
              plotinfo <- save_png(
                plotinfo = plotinfo,
                file = file,
                plotdir = plotdir,
                pwidth = pwidth,
                pheight = pheight,
                punits = punits,
                res = res,
                ptsize = ptsize,
                caption = caption
              )
              recdevfunc3()
              dev.off()
            }
          } # end if uncertinaty
        } # end if print
      } # end if nrow(recdevs)>0
    } # end if max(recdev)>0
    if (!is.null(plotinfo)) {
      plotinfo[["category"]] <- "RecDev"
    }
    return(invisible(plotinfo))
  } # end function