R/plotDblAxis.R

In CamFreshwater/samSim: Closed-loop salmon dynamics simulator

#' Double y-axis Walters plot
#'
#' This function generates a double y-axis plot of the variety popularized by C.
#' Walters to visualize tradeoffs between multiple performance metrics
#' simultaneously. It may be untoward to try place so many different variables
#' with divergents units in the same figure, but hey, people seem to want to see
#' them...
#'
#' @param agDat Dataframe generated by \code{buildDataAgg}.
#' @param keyVar A character value corresponding to the harvest control
#' rule-esqure variable that is on the x-axis (defaults to "expRate") that is
#' used to label x-axis and subset data input data Should correspond to the
#' \code{keyVarName} in \code{buildDataAgg}.
#' @param double A logical that determines whether a double or single y-axis
#' plot is generated.
#' @return Returns a base plot object.
#'
#' @examples
#'
#'
#' @export
plotContTradeOffs <- function(agDat, keyVar = "expRate", double = TRUE) {
  # plot format not compatible with multiple types of HCRs; focus on fixed ER
  # for now
  dum <- agDat %>%
    dplyr::filter(hcr == "fixedER",
                  var %in% c("medSpawners", "medCatch", "ppnCUUpper",
                             "ppnCUExtinct"))

  # What is variable along the x-axis (should be some index of harvest effort)
  dum$keyVar <- as.numeric(as.character(dum[ , keyVar]))
  keyVarRange <- dum %>%
    dplyr::select(keyVar) %>%
    unique()
  xLab <- if (keyVar == "expRate") {
    "Exploitation Rate"
  }

  #subset based on plotting variables
  spwn <- dum %>%
    dplyr::filter(var == "medSpawners")
  catch <- dum %>%
    dplyr::filter(var == "medCatch")
  upp <- dum %>%
    dplyr::filter(var == "ppnCUUpper")
  ext <- dum %>%
    dplyr::filter(var == "ppnCUExtinct")

  #colors for plot
  spwnCol <- "#984ea3"
  uppCol <- "#4daf4a"
  catchCol <- "#377eb8"
  extCol <- "#e41a1c"

  par(mfrow = c(1, 1), mar = c(3.75, 4, 1.25, 4.5), oma = c(0, 0, 0, 0),
      cex.lab = 1)

  plot(1, type="n", axes = FALSE, xlim = c(min(keyVarRange), max(keyVarRange)),
       ylim = c(0, max(dum$highQ)), xlab = "", ylab = "")
  lines(avg ~ keyVar, xaxt = "n", yaxt = "n", type = "l", lwd = 1.25,
        col = spwnCol, data = spwn)
  lines(avg ~ keyVar, xaxt = "n", yaxt = "n", type = "l", lwd = 1.25,
        col = catchCol, data = catch)
  rethinking::shade(rbind(spwn$lowQ, spwn$highQ), spwn$keyVar,
                    col = alpha(spwnCol, 0.3))
  rethinking::shade(rbind(catch$lowQ, catch$highQ), catch$keyVar,
                    col = alpha(catchCol, 0.3))
  axis(1, tick = T, at = c(seq(from = 0,
                               to = max(1, 1.1 * max(keyVarRange)),
                               by = 0.2)),
       mgp = c(3, 0.5, 0))
  mtext(side = 1, line = 2.25, xLab, cex = 1.2)
  if(mean(spwn$avg) > 1000) { #adjust plotting based on scale of spawners
    axis(2, tick = T, at = format(c(-100,
                                    pretty(spwn$highQ, n = 3),
                                    2 * max(spwn$highQ)),
                                  digits = 2, nsmall = 0),
         mgp = c(3, 0.6, 0), las = 0)
    mtext(side = 2, line = 2.5, 'Escapement or Catch', cex = 1.2)
  } else if(mean(spwn$avg) < 1000) {
    axis(2, tick = T, at = round(c(-100,
                                   pretty(spwn$highQ, n = 3),
                                   2 * max(spwn$highQ)),
                                 digits = 1),
         mgp = c(3, 0.6, 0), las = 0)
    mtext(side = 2, line = 2.5, 'Escapement or Catch (millions)', cex = 1.2)
  }
  axis(3, tick = T, at = c(0, max(1, 1.1 * max(keyVarRange))))

  # second set of graphics
  if (double == TRUE) {
    par(new = T)
    plot(avg ~ keyVar, type = "l", axes = F, col = extCol, lty = 2, lwd = 1.25,
         xlab = "", ylab="", ylim = c(0, 1), data = ext)
    rethinking::shade(rbind(ext$lowQ, ext$highQ), ext$keyVar,
                      col = alpha(extCol, 0.3))
    lines(avg ~ keyVar, xaxt="n", yaxt="n", type="l", lty = 2, lwd = 1.25,
          col = uppCol, data = upp)
    rethinking::shade(rbind(upp$lowQ, upp$highQ), upp$keyVar,
                      col = alpha(uppCol, 0.3))
    axis(4, tick = T, at = c(-1, 0.0, 0.5, 1, 2), las = 0, mgp = c(3, 0.6, 0))
    mtext(side = 4, line = 3, 'Proportion of CUs Extinct \n or Above Upper BM',
          cex = 1.2)
  } else {
    axis(4, tick = T, at = round(c(-100 * min(spwn$avg), 2 * max(spwn$highQ)),
                                 digits = 1))
  }
}