r4ss: R Code for Stock Synthesis

Documented in SSplotIndices

#' Plot indices of abundance and associated quantities.
#'
#' Plot indices of abundance with or without model fit as well as other diagnostic
#' plots such as observed vs. expected index and plots related to time-varying
#' catchability (if present).
#'
#'
#' @template replist
#' @param subplots vector controlling which subplots to create
#' Numbering of subplots is as follows, where subplot 9 (comparison of all indices) is
#' provided first:
#' \itemize{
#'   \item 1  index data by fleet
#'   \item 2  index data with fit by fleet
#'   \item 3  observed vs expected index values with smoother
#'   \item 4  index data by fleet on a log scale (lognormal error only)
#'   \item 5  index data with fit by fleet on a log scale (lognormal error only)
#'   \item 6  log(observed) vs log(expected) with smoother (lognormal error only)
#'   \item 7  time series of time-varying catchability (only if actually time-varying)
#'   \item 8  catchability vs. vulnerable biomass (if catchability is not constant)
#'   \item 9  comparison of all indices
#'   \item 10  index residuals based on total uncertainty
#'   \item 11  index residuals based on input uncertainty (not currently provided)
#'   \item 12  index deviations (independent of index uncertainty)
#' }
#'
#' @template plot
#' @template print
#' @template fleets
#' @template fleetnames
#' @param smooth add smoothed line to plots of observed vs. expected sample
#' sizes
#' @param add add to existing plot (not yet implemented)
#' @param datplot make plot of data only?
#' @template labels
#' @param fleetcols vector of colors for all fleets (including those
#' with no index data)
#' @param col1 vector of colors for points in each season for time series plot.
#' Default is red for single season models and a rainbow using the
#' rich.colors.short function for multiple seasons.
#' @param col2 vector of colors for points in each season for obs. vs. exp.
#' plot.  Default is blue for single season models and a rainbow using the
#' rich.colors.short function for multiple seasons.
#' @param col3 color of line showing expected index in time series plot.
#' Default is blue.
#' @param col4 color of smoother shown in obs. vs. exp. plots. Default is red.
#' @param pch1 single value or vector of plotting characters (pch parameter)
#' for time-series plots of index fit. Default=21.
#' @param pch2 single value or vector of plotting characters (pch parameter)
#' for sample size plots of index fit. Default=16.
#' @param cex character expansion factor for points showing observed values.
#' Default=1.
#' @param bg Background color for points with pch=21.
#' @param legend add a legend to seasonal colors (only for seasonal models)
#' @template legendloc
#' @param seasnames optional vector of names for each season to replace
#' defaults if a legend is used
#' @template pwidth
#' @template pheight
#' @template punits
#' @template res
#' @template ptsize
#' @template cex.main
#' @template mainTitle
#' @template plotdir
#' @param minyr First year to show in plot (for zooming in on a subset of
#' values)
#' @param maxyr Last year to show in plot (for zooming in on a subset of
#' values)
#' @param maximum_ymax_ratio Maximum allowed value for ymax (specified
#' as ratio of y), which overrides any
#' value of ymax that is greater (default = Inf)
#' @param show_input_uncertainty Switch controlling whether to add thicker
#' uncertainty interval lines indicating the input uncertainty relative to
#' the total uncertainty which may result from estimating a parameter for
#' extra standard deviations. This is only added for the plots with index
#' fit included (the data-only plots only show the input uncertainty).
#' @template verbose
#' @param \dots Extra arguments to pass to calls to `plot`
#' @author Ian Stewart, Ian Taylor, James Thorson
#' @export
#' @seealso [SS_plots()], [SS_output()]
SSplotIndices <-
  function(
    replist,
    subplots = c(1:10, 12), # IGT 2021/4/15: not sure why 11 is skipped
    plot = TRUE,
    print = FALSE,
    fleets = "all",
    fleetnames = "default",
    smooth = TRUE,
    add = FALSE,
    datplot = TRUE,
    labels = c(
      "Year", # 1
      "Index", # 2
      "Observed index", # 3
      "Expected index", # 4
      "Log index", # 5
      "Log observed index", # 6
      "Log expected index", # 7
      "Standardized index", # 8
      "Catchability (Q)", # 9
      "Time-varying catchability", # 10
      "Vulnerable biomass", # 11
      "Catchability vs. vulnerable biomass", # 12
      "Residual", # 13
      "Deviation"
    ), # 14
    fleetcols = NULL,
    col1 = "default",
    col2 = "default",
    col3 = "blue",
    col4 = "red",
    pch1 = 21,
    pch2 = 16,
    cex = 1,
    bg = "white",
    legend = TRUE,
    legendloc = "topright",
    seasnames = NULL,
    pwidth = 6.5,
    pheight = 5.0,
    punits = "in",
    res = 300,
    ptsize = 10,
    cex.main = 1,
    mainTitle = FALSE,
    plotdir = "default",
    minyr = NULL,
    maxyr = NULL,
    maximum_ymax_ratio = Inf,
    show_input_uncertainty = TRUE,
    verbose = TRUE,
    ...
  ) {
    # get some quantities from replist
    cpue <- replist[["cpue"]]
    SS_versionNumeric <- replist[["SS_versionNumeric"]]

    # confirm that some CPUE values are present
    if (is.null(dim(cpue))) {
      message("skipping index plots: no index data in this model")
      return()
    }

    # table to store information on each plot
    plotinfo <- NULL

    # define a bunch of internal functions

    index.fn <- function(addexpected = TRUE, log = FALSE, ...) {
      # plot of time series of observed values with fit (if requested)

      # don't do anything for log-scale plot if normal error structure is used
      if (error == -1 & log == TRUE) {
        return()
      }

      # function to get uncertainty intervals around points
      # (with or without extra standard error included)
      get_intervals <- function(total = TRUE) {
        if (total) {
          colname <- "SE"
        } else {
          colname <- "SE_input"
        }
        if (error == 0) {
          if (!log) {
            lower <- qlnorm(
              .025,
              meanlog = log(y[include]),
              sdlog = cpueuse[[colname]][include]
            )
            upper <- qlnorm(
              .975,
              meanlog = log(y[include]),
              sdlog = cpueuse[[colname]][include]
            )
          } else {
            lower <- qnorm(
              .025,
              mean = log(y[include]),
              sd = cpueuse[[colname]][include]
            )
            upper <- qnorm(
              .975,
              mean = log(y[include]),
              sd = cpueuse[[colname]][include]
            )
          }
        }
        # normal error interval
        if (error == -1) {
          lower <- qnorm(
            .025,
            mean = y[include],
            sd = cpueuse[[colname]][include]
          )
          upper <- qnorm(
            .975,
            mean = y[include],
            sd = cpueuse[[colname]][include]
          )
        }

        # T-distribution interval
        if (error > 0) {
          lower <- log(y[include]) +
            qt(.025, df = error) * cpueuse[[colname]][include]
          upper <- log(y[include]) +
            qt(.975, df = error) * cpueuse[[colname]][include]
          if (!log) {
            lower <- exp(lower)
            upper <- exp(upper)
          }
        }
        return(data.frame(lower = lower, upper = upper))
      }

      # get total uncertainty
      total_intervals <- get_intervals(total = TRUE)
      lower_total <- total_intervals[["lower"]]
      upper_total <- total_intervals[["upper"]]

      # get input uncertainty
      # SE_input column was first available in SS version 3.30.15 but is calculated
      # elsewhere in this function for models that didn't report it
      input_intervals <- get_intervals(total = FALSE)
      lower_input <- input_intervals[["lower"]]
      upper_input <- input_intervals[["upper"]]

      if (max(upper_total) == Inf) {
        warning(
          "Removing upper interval on indices with infinite upper quantile values.\n",
          "Check the uncertainty inputs for the indices."
        )
        upper_total[upper_total == Inf] <- 100 *
          max(cpueuse[["Obs"]][upper_total == Inf])
      }

      # plot title
      main <- paste0(labels[2], Fleet)
      if (log) {
        main <- paste0(labels[5], Fleet)
      }

      # no title
      if (!mainTitle) {
        main <- ""
      }

      xlim <- c(max(minyr, min(x)), min(maxyr, max(x)))
      if (legend & length(colvec1) > 1) {
        xlim[2] <- xlim[2] + 0.25 * diff(xlim)
      }
      if (!add) {
        # get range for expected values
        zrange <- NULL
        if (addexpected) {
          zrange <- range(z, na.rm = TRUE)
        }
        logzrange <- range(log(z))

        # y-limits for non-log plot
        if (!log) {
          # ylim for standard scale (if lognormal)
          if (error != -1) {
            ylim <- c(
              0,
              1.05 *
                min(
                  max(upper_total, zrange, na.rm = TRUE),
                  max(maximum_ymax_ratio * y, na.rm = TRUE)
                )
            )
          } else {
            ylim <- 1.05 *
              c(
                min(lower_total, zrange, na.rm = TRUE),
                min(
                  max(upper_total, zrange, na.rm = TRUE),
                  max(maximum_ymax_ratio * y, na.rm = TRUE)
                )
              )
          }
        }
        if (log) {
          # ylim for log scale plot
          ylim <- range(c(lower_total, upper_total), na.rm = TRUE)
        }

        plot(
          x = x[include],
          y = y[include],
          type = "n",
          xlab = labels[1],
          ylab = ifelse(!log, labels[2], labels[5]),
          main = main,
          cex.main = cex.main,
          xlim = xlim,
          ylim = ylim,
          yaxs = ifelse(log, "r", "i"),
          ...
        )

        # add line at 0 if it's not a log-scale plot
        # and the axes include zero
        if (!log & min(ylim) < 0) {
          abline(h = 0, lty = 3)
        }
      }

      # set bounds for arrows at total uncertainty
      lower <- lower_total
      upper <- upper_total

      if (addexpected) {
        # show thicker lines behind final lines for input uncertainty
        # only in plot with expected value as well
        if (
          show_input_uncertainty &
            !all(lower_input == lower_total, na.rm = TRUE)
        ) {
          segments(
            x[include],
            lower_input,
            x[include],
            upper_input,
            col = colvec1[s],
            lwd = 3,
            lend = 1
          )
        }
      } else {
        # change bounds for arrows at input uncertainty for plots without fit
        # if values are available
        lower <- lower_input
        upper <- upper_input
      }

      # add intervals
      arrows(
        x0 = x[include],
        y0 = lower,
        x1 = x[include],
        y1 = upper,
        length = 0.03,
        angle = 90,
        code = 3,
        col = colvec1[s]
      )

      # add points and expected values on standard scale
      if (!log) {
        points(
          x = x[include],
          y = y[include],
          pch = pch1,
          cex = cex,
          bg = bg,
          col = colvec1[s]
        )
        if (addexpected) {
          lines(x, z, lwd = 2, col = col3)
          if (length(x) == 1) {
            points(x, z, pch = 23, col = col3)
          }
        }
      } else {
        # add points and expected values on log scale
        points(
          x = x[include],
          y = log(y[include]),
          pch = pch1,
          cex = cex,
          bg = bg,
          col = colvec1[s]
        )
        if (addexpected) {
          lines(x, log(z), lwd = 2, col = col3)
          if (length(x) == 1) {
            points(x, log(z), pch = 23, col = col3)
          }
        }
      }
      if (legend & length(colvec1) > 1) {
        legend(
          x = legendloc,
          legend = seasnames,
          pch = pch1,
          col = colvec1,
          cex = cex
        )
      }
    }

    index_resids.fn <- function(option = 1, ...) {
      # plot of time series of residuals
      # options
      # 1: residuals based on total SE
      # 2: residuals based on input SE
      # 3: deviations (independent of index variability)

      # choose y value and y-axis label

      if (option == 1) {
        # residuals based on total SE
        ylab <- labels[13]
        y <- (log(cpueuse[["Obs"]]) - log(cpueuse[["Exp"]])) / cpueuse[["SE"]]
      }
      if (error == 0 & option == 2) {
        # residuals based on input SE
        ylab <- labels[13]
        # manually calculating residual based on SE_input
        y <- (log(cpueuse[["Obs"]]) - log(cpueuse[["Exp"]])) /
          cpueuse[["SE_input"]]
      }
      if (option == 3) {
        # deviations
        ylab <- labels[14]
        # Dev should be equal to log(Obs/Exp)
        y <- cpueuse[["Dev"]]
      }

      # plot title
      main <- paste(ylab, Fleet)

      # no plot title
      if (!mainTitle) {
        main <- ""
      }
      # xlim (maybe reduced by inputs minyr and maxyr)
      xlim <- c(max(minyr, min(x)), min(maxyr, max(x)))
      if (legend & length(colvec1) > 1) {
        xlim[2] <- xlim[2] + 0.25 * diff(xlim)
      }

      # ylim is symetrical around 0
      ylim <- c(-1.05, 1.05) * max(abs(y[include]))
      if (!add) {
        plot(
          x = x[include],
          y = y[include],
          type = "n",
          xlab = labels[1],
          xlim = xlim,
          ylab = ylab,
          ylim = ylim,
          yaxs = "i",
          main = main,
          cex.main = cex.main,
          ...
        )
      }
      # add points
      points(
        x = x[include],
        y = y[include],
        pch = pch1,
        cex = cex,
        bg = adjustcolor(colvec1[s], alpha.f = 0.7),
        col = adjustcolor(colvec1[s], alpha.f = 0.7)
      )

      # add line at 0
      abline(h = 0, lty = 3)

      # add legend if more than one color (indicating season) was used
      if (legend & length(colvec1) > 1) {
        legend(
          x = legendloc,
          legend = seasnames,
          pch = pch1,
          pt.bg = colvec1,
          col = colvec1,
          cex = cex
        )
      }
    }

    obs_vs_exp.fn <- function(log = FALSE, ...) {
      # plot of observed vs. expected with smoother

      # plot title
      main <- paste(labels[2], Fleet, sep = " ")
      # no title
      if (!mainTitle) {
        main <- ""
      }

      if (!add) {
        if (!log) {
          # standard plot
          plot(
            y[include],
            z[include],
            type = "n",
            xlab = labels[3],
            ylab = labels[4],
            main = main,
            cex.main = cex.main,
            ylim = c(0, 1.05 * max(z)),
            xlim = c(0, 1.05 * max(y)),
            xaxs = "i",
            yaxs = "i",
            ...
          )
        } else {
          # log-scale plot doesn't specificy y limits
          plot(
            log(y[include]),
            log(z[include]),
            type = "n",
            xlab = labels[6],
            ylab = labels[7],
            main = main,
            cex.main = cex.main
          )
        }
      }
      if (!log) {
        points(y[include], z[include], col = colvec2[s], pch = pch2, cex = cex)
      } else {
        points(
          log(y[include]),
          log(z[include]),
          col = colvec2[s],
          pch = pch2,
          cex = cex
        )
      }
      abline(a = 0, b = 1, lty = 3)
      if (smooth && npoints > 6 && diff(range(y)) > 0) {
        if (!log) {
          psmooth <- loess(z[include] ~ y[include], degree = 1)
          lines(
            psmooth[["x"]][order(psmooth[["x"]])],
            psmooth[["fitted"]][order(psmooth[["x"]])],
            lwd = 1.2,
            col = col4,
            lty = "dashed"
          )
        } else {
          psmooth <- loess(log(z[include]) ~ log(y[include]), degree = 1)
          lines(
            psmooth[["x"]][order(psmooth[["x"]])],
            psmooth[["fitted"]][order(psmooth[["x"]])],
            lwd = 1.2,
            col = col4,
            lty = "dashed"
          )
        }
      }
      if (legend & length(colvec2) > 1) {
        legend(
          x = legendloc,
          legend = seasnames,
          pch = pch2,
          col = colvec2,
          cex = cex
        )
      }
    }

    timevarying_q.fn <- function() {
      # plot of time-varying catchability (if present)
      main <- paste(labels[10], Fleet, sep = " ")
      if (!mainTitle) {
        main <- ""
      }
      q <- cpueuse[["Calc_Q"]]
      if (!add) {
        plot(
          x,
          q,
          type = "o",
          xlab = labels[1],
          main = main,
          cex.main = cex.main,
          ylab = labels[9],
          col = colvec2[1],
          pch = pch2
        )
      }
    }

    q_vs_vuln_bio.fn <- function() {
      # plot of time-varying catchability (if present)
      main <- paste(labels[12], Fleet, sep = " ")
      if (!mainTitle) {
        main <- ""
      }
      v <- cpueuse[["Vuln_bio"]]
      q1 <- cpueuse[["Calc_Q"]]
      q2 <- cpueuse[["Eff_Q"]]
      if (all(q1 == q2)) {
        ylab <- labels[9]
      } else {
        ylab <- "Effective catchability"
      }
      if (!add) {
        plot(
          v,
          q2,
          type = "o",
          xlab = labels[11],
          main = main,
          cex.main = cex.main,
          ylab = ylab,
          col = colvec2[1],
          pch = pch2
        )
      }
    }

    # check for super periods
    if (length(grep("supr_per", cpue[["Supr_Per"]]))) {
      warning(
        "Some indices have superperiods. Values will be plotted\n",
        "in year/season associated with data in report file."
      )
      cpue <- cpue[!is.na(cpue[["Dev"]]), ]
    }

    FleetNames <- replist[["FleetNames"]]
    nfleets <- replist[["nfleets"]]
    nseasons <- replist[["nseasons"]]

    # find any extra SD parameters
    parameters <- replist[["parameters"]]
    Q_extraSD_info <- parameters[grep("Q_extraSD", parameters[["Label"]]), ]
    # calculate how many of these parameters there are
    nSDpars <- nrow(Q_extraSD_info)
    if (nSDpars > 0) {
      # parse the parameter label to get the fleet number
      Q_extraSD_info[["Fleet"]] <- NA
      for (ipar in 1:nSDpars) {
        if (SS_versionNumeric >= 3.3) {
          # parsing label with ending like "(2)" assuming only one set of parentheses
          num <- strsplit(
            Q_extraSD_info[["Label"]][ipar],
            split = "[()]",
            fixed = FALSE
          )[[1]][2]
        } else {
          num <- strsplit(
            substring(Q_extraSD_info[["Label"]][ipar], nchar("Q_extraSD_") + 1),
            split = "_",
            fixed = TRUE
          )[[1]][1]
        }
        Q_extraSD_info[["Fleet"]][ipar] <- as.numeric(num)
      }
      # NOTE: important columns in Q_extraSD_info to use below are $Value and $Fleet
    }
    if (nseasons > 1) {
      # if seasons, put CPUE at season midpoint
      cpue[["YrSeas"]] <- cpue[["Yr"]] + (cpue[["Seas"]] - 0.5) / nseasons
    } else {
      # if no seasons, put at integer year value
      cpue[["YrSeas"]] <- cpue[["Yr"]]
    }
    if (plotdir == "default") {
      plotdir <- replist[["inputs"]][["dir"]]
    }

    if (fleetnames[1] == "default") {
      fleetnames <- FleetNames
    }
    if (fleets[1] == "all") {
      fleets <- 1:nfleets
    } else {
      if (length(intersect(fleets, 1:nfleets)) != length(fleets)) {
        return(
          "Input 'fleets' should be 'all' or a vector of values between 1 and nfleets."
        )
      }
    }

    # subset fleets as requested
    fleetvec <- intersect(fleets, unique(as.numeric(cpue[["Fleet"]])))

    # empty data.frame to store data for comparison among indices
    allcpue <- data.frame()
    # keep track of whether any indices with negative observations is excluded
    any_negative <- FALSE

    # loop over fleets
    for (ifleet in fleetvec) {
      # use fancy colors only if the individual index spans more than one season
      usecol <- FALSE
      if (length(unique(cpue[["Seas"]][cpue[["Fleet"]] == ifleet])) > 1) {
        usecol <- TRUE
      }

      # turn off use of legend if there's never more than 1 season per index
      if (!usecol) {
        legend <- FALSE
      }

      if (col1[1] == "default") {
        colvec1 <- "black"
        if (usecol & nseasons == 4) {
          colvec1 <- c("blue4", "green3", "orange2", "red3")
        }
        if (usecol & !nseasons %in% c(1, 4)) {
          colvec1 <- rich.colors.short(nseasons)
        }
      } else {
        colvec1 <- col1
        # if user provides single value (or vector of length less than nseasons)
        # make sure it's adequate to cover all seasons
        if (length(colvec1) < nseasons) {
          colvec1 <- rep(col1, nseasons)
        }
      }
      if (col2[1] == "default") {
        colvec2 <- "blue"
        if (usecol & nseasons == 4) {
          colvec2 <- c("blue4", "green3", "orange2", "red3")
        }
        if (usecol & !nseasons %in% c(1, 4)) {
          colvec2 <- rich.colors.short(nseasons)
        }
      } else {
        colvec2 <- col2
        # if user provides single value (or vector of length less than nseasons)
        # make sure it's adequate to cover all seasons
        if (length(colvec1) < nseasons) {
          colvec1 <- rep(col1, nseasons)
        }
      }
      if (is.null(seasnames)) {
        seasnames <- paste("Season", 1:nseasons, sep = "")
      }

      Fleet <- fleetnames[ifleet]
      error <- replist[["survey_error"]][ifleet]
      if (error == 0) {
        error_caption <- "lognormal error"
      }
      if (error == -1) {
        error_caption <- "normal error"
      }
      if (error == 1) {
        error_caption <- paste0(
          "T-distributed error with ",
          error,
          " degree of freedom"
        )
      }
      if (error > 1) {
        error_caption <- paste0(
          "T-distributed error with ",
          error,
          " degrees of freedom"
        )
      }

      cpueuse <- cpue[cpue[["Fleet"]] == ifleet, ]
      cpueuse <- cpueuse[order(cpueuse[["YrSeas"]]), ]

      # look for time-vary
      time <- diff(range(cpueuse[["Calc_Q"]])) > 0
      # look for time-varying effective Q
      time2 <- diff(range(cpueuse[["Eff_Q"]])) > 0
      # Teresa's model had NA values in Eff_Q for unknown reasons
      # line below will allow model to play on
      if (is.na(time2)) {
        time2 <- FALSE
      }
      # if "SE_input" column not available, look for extra SD and
      # calculate input SD (if different from final value)
      if (!"SE_input" %in% names(cpue)) {
        if (exists("Q_extraSD_info") && ifleet %in% Q_extraSD_info[["Fleet"]]) {
          # input uncertainty is final value minus extra SD parameter (if present)
          cpueuse[["SE_input"]] <- cpueuse[["SE"]] -
            Q_extraSD_info[["Value"]][Q_extraSD_info[["Fleet"]] == ifleet]
        } else {
          cpueuse[["SE_input"]] <- cpueuse[["SE"]]
        }
      }
      # use short variable names for often-used quantities
      x <- cpueuse[["YrSeas"]]
      y <- cpueuse[["Obs"]]
      z <- cpueuse[["Exp"]]
      npoints <- length(z)
      include <- !is.na(cpueuse[["Like"]])
      if (any(include)) {
        if (usecol) {
          s <- cpueuse[["Seas"]][which(include)]
        } else {
          s <- 1 # only use colorvector if more than 1 season
        }
        if (datplot) {
          # add index data to data frame which is used to compare all indices
          if (min(cpueuse[["Obs"]] >= 0)) {
            cpueuse[["Index"]] <- rep(ifleet, length(cpueuse[["YrSeas"]]))
            cpueuse[["stdvalue"]] <- cpueuse[["Obs"]] / mean(cpueuse[["Obs"]])
            tempcpue <- cbind(
              cpueuse[["Index"]],
              cpueuse[["YrSeas"]],
              cpueuse[["Obs"]],
              cpueuse[["stdvalue"]]
            )
            colnames(tempcpue) <- c("Index", "year", "value", "stdvalue")
            allcpue <- rbind(allcpue, tempcpue)
          } else {
            if (verbose & 9 %in% subplots & datplot) {
              message(
                "Excluding fleet ",
                ifleet,
                " from index comparison figure because it has negative values"
              )
            }
            any_negative <- TRUE
          }
        }

        addlegend <- function(pch, colvec) {
          names <- paste(seasnames, "observations")
        }

        if (plot) {
          if (1 %in% subplots & datplot) {
            index.fn(addexpected = FALSE)
          }
          if (2 %in% subplots) {
            index.fn()
          }
          if (3 %in% subplots) {
            obs_vs_exp.fn()
          }
        }
        if (print) {
          if (1 %in% subplots & datplot) {
            file <- paste0("index1_cpuedata_", gsub(" ", "", Fleet), ".png")
            caption <- paste0(
              "Index data for ",
              Fleet,
              ". ",
              "Lines indicate 95% uncertainty interval around index values ",
              "based on the model assumption of ",
              error_caption,
              ". ",
              "Thicker lines (if present) indicate input uncertainty before addition of ",
              "estimated additional uncertainty parameter."
            )
            plotinfo <- save_png(
              plotinfo = plotinfo,
              file = file,
              plotdir = plotdir,
              pwidth = pwidth,
              pheight = pheight,
              punits = punits,
              res = res,
              ptsize = ptsize,
              caption = caption
            )
            index.fn(addexpected = FALSE)
            dev.off()
          }
          if (2 %in% subplots) {
            file <- paste0("index2_cpuefit_", gsub(" ", "", Fleet), ".png")
            caption <- paste0(
              "Fit to index data for ",
              Fleet,
              ". ",
              "Lines indicate 95% uncertainty interval around index values ",
              "based on the model assumption of ",
              error_caption,
              ". ",
              "Thicker lines (if present) indicate input uncertainty before addition of ",
              "estimated additional uncertainty parameter."
            )
            plotinfo <- save_png(
              plotinfo = plotinfo,
              file = file,
              plotdir = plotdir,
              pwidth = pwidth,
              pheight = pheight,
              punits = punits,
              res = res,
              ptsize = ptsize,
              caption = caption
            )
            index.fn()
            dev.off()
          }
          if (3 %in% subplots) {
            file <- paste0("index3_obs_vs_exp_", gsub(" ", "", Fleet), ".png")
            caption <- paste(
              "Observed vs. expected index values with smoother for",
              Fleet
            )
            plotinfo <- save_png(
              plotinfo = plotinfo,
              file = file,
              plotdir = plotdir,
              pwidth = pwidth,
              pheight = pheight,
              punits = punits,
              res = res,
              ptsize = ptsize,
              caption = caption
            )
            obs_vs_exp.fn()
            dev.off()
          }
        }

        # same plots again in log space
        # check for lognormal error
        if (error != -1) {
          # plot subplots 4-6 to current device
          if (plot) {
            if (4 %in% subplots & datplot) {
              index.fn(log = TRUE, addexpected = FALSE)
            }
            if (5 %in% subplots) {
              index.fn(log = TRUE)
            }
            if (6 %in% subplots) {
              obs_vs_exp.fn(log = TRUE)
            }
          }

          # print subplots 4-6 to PNG files
          if (print) {
            if (4 %in% subplots & datplot) {
              file <- paste0(
                "index4_logcpuedata_",
                gsub(" ", "", Fleet),
                ".png"
              )
              caption <- paste0(
                "Log index data for ",
                Fleet,
                ". ",
                "Lines indicate 95% uncertainty interval around index values ",
                "based on the model assumption of ",
                error_caption,
                ". ",
                "Thicker lines (if present) indicate input uncertainty before addition of ",
                "estimated additional uncertainty parameter."
              )
              plotinfo <- save_png(
                plotinfo = plotinfo,
                file = file,
                plotdir = plotdir,
                pwidth = pwidth,
                pheight = pheight,
                punits = punits,
                res = res,
                ptsize = ptsize,
                caption = caption
              )
              index.fn(log = TRUE, addexpected = FALSE)
              dev.off()
            }
            if (5 %in% subplots) {
              file <- paste0("index5_logcpuefit_", gsub(" ", "", Fleet), ".png")
              caption <- paste0(
                "Fit to log index data on log scale for ",
                Fleet,
                ". ",
                "Lines indicate 95% uncertainty interval around index values ",
                "based on the model assumption of ",
                error_caption,
                ". ",
                "Thicker lines (if present) indicate input uncertainty before addition of ",
                "estimated additional uncertainty parameter."
              )
              plotinfo <- save_png(
                plotinfo = plotinfo,
                file = file,
                plotdir = plotdir,
                pwidth = pwidth,
                pheight = pheight,
                punits = punits,
                res = res,
                ptsize = ptsize,
                caption = caption
              )
              index.fn(log = TRUE)
              dev.off()
            }
            if (6 %in% subplots) {
              file <- paste0(
                "index6_log_obs_vs_exp_",
                gsub(" ", "", Fleet),
                ".png"
              )
              caption <- paste(
                "log(observed) vs. log(expected) index values with smoother for",
                Fleet
              )
              plotinfo <- save_png(
                plotinfo = plotinfo,
                file = file,
                plotdir = plotdir,
                pwidth = pwidth,
                pheight = pheight,
                punits = punits,
                res = res,
                ptsize = ptsize,
                caption = caption
              )
              obs_vs_exp.fn(log = TRUE)
              dev.off()
            }
          }
        } # end plots that require lognormal error

        # plots 7 and 8 related to time-varying catchability
        if (plot) {
          if (7 %in% subplots & time) {
            timevarying_q.fn()
          }
          if (8 %in% subplots & time2) {
            q_vs_vuln_bio.fn()
          }
        } # end plot to graphics device

        if (print) {
          if (7 %in% subplots & time) {
            file <- paste0(
              "index7_timevarying_q_",
              gsub(" ", "", Fleet),
              ".png"
            )
            caption <- paste("Timeseries of catchability for", Fleet)
            plotinfo <- save_png(
              plotinfo = plotinfo,
              file = file,
              plotdir = plotdir,
              pwidth = pwidth,
              pheight = pheight,
              punits = punits,
              res = res,
              ptsize = ptsize,
              caption = caption
            )
            timevarying_q.fn()
            dev.off()
          }
          if (8 %in% subplots & time2) {
            file <- paste0(
              "index8_q_vs_vuln_bio_",
              gsub(" ", "", Fleet),
              ".png"
            )
            caption <-
              paste0(
                "Catchability vs. vulnerable biomass for fleet ",
                Fleet,
                "<br> \n",
                "This plot should illustrate curvature of nonlinear catchability relationship<br> \n",
                "or reveal patterns associated with random-walk catchability."
              )
            plotinfo <- save_png(
              plotinfo = plotinfo,
              file = file,
              plotdir = plotdir,
              pwidth = pwidth,
              pheight = pheight,
              punits = punits,
              res = res,
              ptsize = ptsize,
              caption = caption
            )
            q_vs_vuln_bio.fn()
            dev.off()
          }
        } # end print to PNG

        # residual/deviation plots
        if (plot) {
          if (10 %in% subplots & all(cpueuse[["Obs"]] >= 0)) {
            index_resids.fn(option = 1)
          }
          if (11 %in% subplots) {
            index_resids.fn(option = 2)
          }
          if (12 %in% subplots) {
            index_resids.fn(option = 3)
          }
        }
        if (print) {
          #### residuals based on total uncertainty
          if (10 %in% subplots & all(cpueuse[["Obs"]] >= 0)) {
            file <- paste0(
              "index10_resids_SE_total_",
              gsub(" ", "", Fleet),
              ".png"
            )
            caption <- paste0("Residuals of fit to index for ", Fleet, ".")
            if (error == 0) {
              caption <- paste0(
                caption,
                "<br>Values are (log(Obs) - log(Exp))/SE ",
                "where SE is the total standard error including any ",
                "estimated additional uncertainty."
              )
            } else {
              caption <- paste0(
                caption,
                "<br>Values are based on the total standard error ",
                "including any estimated additional uncertainty."
              )
            }
            plotinfo <- save_png(
              plotinfo = plotinfo,
              file = file,
              plotdir = plotdir,
              pwidth = pwidth,
              pheight = pheight,
              punits = punits,
              res = res,
              ptsize = ptsize,
              caption = caption
            )
            index_resids.fn(option = 1)
            dev.off()
          }
          #### residuals based on input uncertainty
          if (
            11 %in%
              subplots &
              show_input_uncertainty &&
              any(!is.null(cpueuse[["SE_input"]][include])) &&
              any(cpueuse[["SE_input"]] > cpueuse[["SE"]])
          ) {
            file <- paste0(
              "index11_resids_SE_input_",
              gsub(" ", "", Fleet),
              ".png"
            )
            caption <- paste0("Residuals for fit to index for ", Fleet, ".")
            if (error == 0) {
              caption <- paste0(
                caption,
                "<br>Values are (log(Obs) - log(Exp))/SE_input ",
                "where SE_input is the input standard error",
                "excluding any estimated additional uncertainty."
              )
            } else {
              caption <- paste0(
                caption,
                "<br>Values are based on the input standard error ",
                "excluding any estimated additional uncertainty."
              )
            }
            plotinfo <- save_png(
              plotinfo = plotinfo,
              file = file,
              plotdir = plotdir,
              pwidth = pwidth,
              pheight = pheight,
              punits = punits,
              res = res,
              ptsize = ptsize,
              caption = caption
            )
            index_resids.fn(option = 2)
            dev.off()
          }
          #### simple deviation plot
          if (12 %in% subplots) {
            file <- paste0(
              "index12_resids_SE_total_",
              gsub(" ", "", Fleet),
              ".png"
            )
            caption <- paste0("Deviations for fit to index for ", Fleet, ".")
            if (error != -1) {
              # lognormal or T-distributed error
              caption <- paste0(
                caption,
                "<br>Values are log(Obs) - log(Exp) ",
                "and thus independent of index uncertainty."
              )
            }
            if (error == -1) {
              # normal error
              caption <- paste0(
                caption,
                "<br>Values are Obs - Exp ",
                "and thus independent of index uncertainty."
              )
            }
            plotinfo <- save_png(
              plotinfo = plotinfo,
              file = file,
              plotdir = plotdir,
              pwidth = pwidth,
              pheight = pheight,
              punits = punits,
              res = res,
              ptsize = ptsize,
              caption = caption
            )
            index_resids.fn(option = 3)
            dev.off()
          }
        } # end if(print)
      } # end check for any values to include
    } # end loop over fleets

    ### standardized plot of all CPUE indices
    if (datplot == TRUE & nrow(allcpue) > 0) {
      all_index.fn <- function() {
        main <- "All index plot"
        if (!mainTitle) {
          main <- ""
        }
        xlim <- c(
          min(allcpue[["year"]], na.rm = TRUE) - 1,
          max(allcpue[["year"]], na.rm = TRUE) + 1
        )

        # change year range if requested
        xlim[1] <- max(xlim[1], minyr)
        xlim[2] <- min(xlim[2], maxyr)

        # set y limits
        ylim <- c(0, 1.05 * max(allcpue[["stdvalue"]], na.rm = TRUE))
        # set colors
        if (!is.null(fleetcols) & length(fleetcols) >= nfleets) {
          usecols <- fleetcols
        } else {
          usecols <- rich.colors.short(
            max(allcpue[["Index"]], na.rm = TRUE),
            alpha = 0.7
          )
          if (max(allcpue[["Index"]], na.rm = TRUE) >= 2) {
            usecols <- rich.colors.short(
              max(allcpue[["Index"]], na.rm = TRUE) + 1,
              alpha = 0.7
            )[-1]
          }
        }
        # make empty plot
        if (!add) {
          plot(
            0,
            type = "n",
            xlab = labels[1],
            main = main,
            cex.main = cex.main,
            col = usecols[1],
            ylab = labels[8],
            xlim = xlim,
            ylim = ylim,
            yaxs = "i"
          )
        }
        # add points and lines for each fleet
        for (ifleet in fleetvec) {
          lines(
            x = allcpue[["year"]][allcpue[["Index"]] == ifleet],
            y = allcpue[["stdvalue"]][allcpue[["Index"]] == ifleet],
            col = adjustcolor(usecols[ifleet], alpha.f = 0.7),
            lwd = 2
          )
          points(
            x = allcpue[["year"]][allcpue[["Index"]] == ifleet],
            y = allcpue[["stdvalue"]][allcpue[["Index"]] == ifleet],
            pch = pch1,
            bg = adjustcolor(usecols[ifleet], alpha.f = 0.7),
            col = gray(0, alpha = 0.7),
            cex = cex
          )
        }
        legend(
          legendloc,
          legend = fleetnames[fleetvec],
          ncol = 2,
          bty = "n",
          pch = pch1,
          col = gray(0, alpha = 0.7),
          pt.bg = usecols[fleetvec]
        )
      } # end all_index.fn
      if (plot & (9 %in% subplots)) {
        all_index.fn()
      }
      if (print & (9 %in% subplots)) {
        file <- paste0("index9_standcpueall", ".png")
        caption <- paste(
          "Standardized indices overlaid.",
          "Each index is rescaled to have mean observation = 1.0."
        )
        if (any_negative) {
          caption <- paste(
            caption,
            "Indices with negative observations have been excluded."
          )
        }
        plotinfo <- save_png(
          plotinfo = plotinfo,
          file = file,
          plotdir = plotdir,
          pwidth = pwidth,
          pheight = pheight,
          punits = punits,
          res = res,
          ptsize = ptsize,
          caption = caption
        )
        all_index.fn()
        dev.off()
      }
    } # end datplot

    if (!is.null(plotinfo)) {
      plotinfo[["category"]] <- "Index"
    }
    return(invisible(plotinfo))
  } # end function