R/SSplotIndices.R

In r4ss: R Code for Stock Synthesis

#' 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)
#' }
#'
#' @param plot plot to active plot device?
#' @param print print to PNG files?
#' @param fleets optional vector to subset fleets for which plots will be made
#' @param fleetnames optional replacement for fleenames used in data file
#' @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?
#' @param labels vector of labels for plots (titles and axis 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)
#' @param legendloc add a legend to seasonal colors (default is "topright")
#' @param seasnames optional vector of names for each season to replace
#' defaults if a legend is used
#' @param pwidth width of plot
#' @param pheight height of plot
#' @param punits units for PNG file
#' @template res
#' @param ptsize point size for PNG file
#' @param cex.main character expansion for plot titles
#' @param mainTitle switch which allows the plot title to be left off
#' @param plotdir directory where PNG files will be written. by default it will
#' be the directory where the model was run.
#' @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).
#' @param verbose report progress to R GUI?
#' @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)
            ))
          } 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)
              )
            )
          }
        }
        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"),
          ...
        )
      }

      # 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_", 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_", 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_", 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_", 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_", 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_", 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_", 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_", 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_", 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_", 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_", 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("top",
          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