R/plotTradeoffs.R

In CamFreshwater/samSim: Closed-loop salmon dynamics simulator

#' CU-specific trade off plot
#'
#' This function generates Kobe style plots for a series of conservation- and
#' catch-based performance metrics.
#'
#' @importFrom dplyr everything filter mutate select
#' @importFrom ggplot2 aes facet_wrap geom_errorbar geom_errorbarh geom_point
#' ggplot guides labs scale_alpha_discrete scale_shape_manual scale_fill_manual
#' theme
#' @importFrom tidyr gather spread unite
#'
#' @param cuDat Dataframe generated by \code{buildCUDat}.
#' @param consVar A character value corresponding to a conservation based PM
#' in cuDat$vars.
#' @param catchVar A character value corresponding to a catch based PM in
#' cuDat$vars.
#' @param facet A character value that can take the values:
#' \code{"cu", "mp", "om"} and specifies along which categorical variable
#' the plot should be faceted.
#' @param panel  A character value that can take the values:
#' \code{"mp", "om"} and specifies along which categorical variable new
#' pages in the output PDF will be generated.
#' @param showUncertainty A logical specifying whether whiskers for each
#' variable's credible interval should be plotted.
#' @param hotColors A logical (default \code{TRUE}) that specifies whether
#' symbols should be filled with \code{viridis} palette or grey scale.
#' @param legendLab A character representing the legend title.
#' @param xLab A character representing the x axis label.
#' @param yLab A character representing the y axis label.
#' @param main A logical specifying whether a plot title should be added.
#' @param scaleAxis A character vector that can take values `c("fixed", "free",
#' "free_x", "free_y")` that determines which axes, if any, have variable axes
#' dimensions across facets.
#' @return Returns a ggplot object.
#'
#' @examples
#' plotCUTradeoff(cuPlottingDF, consVar = "medSpawners", catchVar = "medCatch",
#' facet = "cu", panel = "om", showUncertainty = FALSE,
#' legendLab = "Prop. TAC in mixed stock fishery", xLab = "Median Catch",
#' yLab = "Median Spawners", main = FALSE)
#'
#' @export
plotCUTradeoff <- function(cuDat, consVar = "medSpawners", catchVar = "medCatch",
                           facet = "cu", panel = NULL, showUncertainty = FALSE,
                           hotColors = TRUE, legendLab = NULL, xLab = NULL,
                           yLab = NULL, main = TRUE,
                           axisSize = 14, dotSize = 4, lineSize = 1.25,
                           legendSize = 14, freeY = TRUE, scaleAxis = "free") {
  xLab <- ifelse(is.null(xLab), catchVar, xLab)
  yLab <- ifelse(is.null(yLab), consVar, yLab)
  #save index variables
  nCU <- length(unique(cuDat$cuName))

  #identify whether second dimension of plots should be by om or MP
  #(first dimension is by keyvariable, faceting is by CU/MU)
  panels <- if (is.null(panel)) {
    NA
  } else if (panel == "om") {
    unique(cuDat$om)
  } else if (panel == "mp") {
    unique(cuDat$mp)
  }

  # Plot
  plotList <- lapply(seq_along(panels), function(h) { #iterate across catch variables
    dum <- if(is.null(panel)) {
      cuDat
    } else if (panel == "om") {
      cuDat %>%
        filter(om == panels[h])
    } else if (panel == "mp") {
      cuDat %>%
        filter(mp == panels[h])
    }

    plotTitle <- ifelse(main == TRUE, paste(panels[h], "Plot", sep = ""), "")

    dum <- dum %>% filter(var == catchVar | var == consVar)
    dum$var <- plyr::mapvalues(dum$var, from = c(consVar, catchVar), #change factor names to make plotting universal
                               to = c("consVar", "catchVar"))
    #necessary to spread for tradeoff plots; NOTE: if errors, check indexing correct)
    wideDum <- dum %>%
      gather(temp, value, avg, lowQ, highQ) %>%
      unite(temp1, var, temp, sep = "_") %>%
      spread(temp1, value) %>%
      dplyr::select(keyVar = 1, everything()) %>%
      mutate(keyVar = as.factor(keyVar))

    #identify faceting
    if (facet == "cu") {
      wideDum <- wideDum %>%
        mutate(facetVar = as.factor(cuName))
    }
    if (facet == "mp") {
      wideDum <- wideDum %>%
        mutate(facetVar = as.factor(mp))
    }
    if (facet == "om") {
      wideDum <- wideDum %>%
        mutate(facetVar = as.factor(om))
    }

    #groupDat for axis break limits
    axBreaks <-
      if (hotColors == TRUE) {
        colPal <- viridis::viridis(length(levels(wideDum$keyVar)),
                                         begin = 0, end = 1, option = "plasma")
        names(colPal) <- levels(wideDum$keyVar)
      } else if (hotColors == FALSE) {
        colPal <- grDevices::gray.colors(n = length(levels(wideDum$keyVar)),
                                   start = 0.9, end = 0.05)
        names(colPal) <- levels(wideDum$keyVar)
      }
    p <- ggplot(wideDum, aes(x = catchVar_avg, y = consVar_avg, shape = hcr,
                             fill = keyVar)) +
      geom_point(size = dotSize) +
      scale_fill_manual(values = colPal, name = legendLab) +
      guides(fill = guide_legend(override.aes = list(shape = 21))) +
      theme_sleekX() +
      theme(strip.text = element_text(size = axisSize),
            axis.text = element_text(size = 0.9 * axisSize),
            axis.title = element_text(size = axisSize),
            legend.text = element_text(size = 0.9 * legendSize),
            legend.title = element_text(size = legendSize)) +
      labs(x = xLab, y = yLab, title = plotTitle) +
      scale_shape_manual(values = c(21, 25), name = "Control Rule") +
      facet_wrap(~ facetVar, scales = scaleAxis) +
      scale_x_continuous(breaks = scales::pretty_breaks(n = 3)) +
      scale_y_continuous(breaks = scales::pretty_breaks(n = 3))

    if (length(unique(wideDum$hcr)) < 2) {
      p <- p +
        guides(shape = "none")
    } else {
      p <- p +
        guides(shape = guide_legend(override.aes = list(fill = "black")))
    }
    if (showUncertainty == FALSE) {
      return(p)
    }
    if (showUncertainty==TRUE) {
      q <- p +
        geom_errorbar(aes(ymin = consVar_lowQ, ymax = consVar_highQ),
                      alpha = 0.3, width = 0, size = lineSize) +
        geom_errorbarh(aes(xmin = catchVar_lowQ, xmax = catchVar_highQ),
                       alpha = 0.3, height = 0, size = lineSize)
      return(q)
    }
  }) #end panel lapply subset
  names(plotList) <- sapply(panels, function(x) paste(x, "Plot", sep = ""))
  return(plotList)
}

#______________________________________________________________________________

#' Aggregate trade off plot
#'
#' This function generates Kobe style plots for a series of conservation- and
#' catch-based performance metrics.
#'
#' @importFrom dplyr everything filter mutate select
#' @importFrom ggplot2 aes facet_wrap geom_errorbar geom_errorbarh geom_point
#' ggplot guides labs scale_alpha_discrete scale_shape_manual scale_fill_manual
#' theme
#' @importFrom tidyr gather spread unite
#'
#' @param agDat Dataframe generated by \code{buildAgDat}.
#' @param consVar A character value corresponding to a conservation based PM
#' in cuDat$vars.
#' @param catchVar A character value corresponding to a catch based PM in
#' cuDat$vars.
#' @param facet A character value that can take the values: \code{"mp", "om"}
#' and specifies along which categorical variable the plot should be faceted.
#' @param shape A character value that defaults to \code{NULL}, but can take
#' values \code{"mp"} or \code{"om"}, and specifies along which categorical
#' variable shapes should be plotted. Note maximum number of levels is 5.
#' @param hotColors A logical (default \code{TRUE}) that specifies whether
#' symbols should be filled with \code{viridis} palette or solid black.
#' @param showUncertainty A logical specifying whether whiskers for each
#' variables credible interval should be plotted.
#' @param legendLab A character representing the legend title.
#' @param xLab A character representing the x axis label.
#' @param yLab A character representing the y axis label.
#' @param mainLab A character specifying a plot title (defaults to NULL).
#' @param scaleAxis A character vector that can take values `c("fixed", "free",
#' "free_x", "free_y")` that determines which axes, if any, have variable axes dimensions
#' across facets.
#' @param facetLetter A logical that determines facets are labelled with letters
#' for referencing in text.
#' @return Returns a ggplot object.
#'
#' @examples
#' plotAgTradeoff(agPlottingDF, consVar = "medSpawners", catchVar = "medCatch",
#' facet = "om", showUncertainty = TRUE,
#' legendLab = "Prop. TAC in mixed stock fishery", xLab = "Median Catch",
#' yLab = "Median Spawners")
#'
#' @export
plotAgTradeoff <- function(agDat, consVar = "medSpawners",
                           catchVar = "medCatch", facet = "om", shape = NULL,
                           hotColors = TRUE, showUncertainty = FALSE,
                           legendLab = NULL, xLab = NULL, yLab = NULL,
                           mainLab = NULL, axisSize = 14, dotSize = 4,
                           lineSize = 1.25, legendSize = 14, scaleAxis = "free",
                           facetLetter = FALSE) {
  xLab <- ifelse(is.null(xLab), catchVar, xLab)
  yLab <- ifelse(is.null(yLab), consVar, yLab)

  dum <- agDat %>%
    dplyr::filter(var == catchVar | var == consVar)
  #change factor names to make plotting universal
  dum$var <- plyr::mapvalues(dum$var, from = c(consVar, catchVar),
                             to = c("consVar", "catchVar"))

  #necessary to spread for tradeoff plots; NOTE: if errors, check indexing correct)
  wideDum <- dum %>%
    gather(temp, value, avg, lowQ, highQ) %>%
    unite(temp1, var, temp, sep = "_") %>%
    spread(temp1, value) %>%
    dplyr::select(keyVar = 1, everything()) %>%
    mutate(keyVar = as.factor(keyVar))

  #identify faceting and shape variables
  if (length(unique(wideDum$mp)) == 1 & is.null(wideDum$hcr)) {
    wideDum <- wideDum %>%
      mutate(hcr = mp)
  }

  if (!is.null(facet)) {
    if (facet == "mp") {
      wideDum <- wideDum %>%
        mutate(facetVar = as.factor(mp))
    } else if (facet == "om") {
      wideDum <- wideDum %>%
        mutate(facetVar = as.factor(om))
    }
  }

  if (is.null(shape)) {
    wideDum <- wideDum %>%
      mutate(shapeVar = as.factor(hcr))
    secLegendLab = "Harvest\nControl Rule"
  } else if (shape == "hcr") {
    wideDum <- wideDum %>%
      mutate(shapeVar = as.factor(hcr))
    secLegendLab = "Harvest\nControl Rule"
  } else if (shape == "om") {
    wideDum <- wideDum %>%
      mutate(shapeVar = as.factor(om))
    secLegendLab = "Operating Model"
  } else if (shape == "mp") {
    wideDum <- wideDum %>%
      mutate(shapeVar = as.factor(mp))
    secLegendLab = "Fixed\nExploitation Rate"
  }

  if (length(levels(wideDum$shapeVar)) > 5) {
    warning("Too many factor levels to plot as shapes, switch to facet")
  } else {
    shapePalette <- c(21,25,23,22,24)
    names(shapePalette) <- levels(wideDum$shapeVar)
  }

  if (hotColors == TRUE) {
    colPal <- viridis::viridis(length(unique(wideDum$keyVar)), begin = 0,
                                  end = 1, option = "plasma")
    names(colPal) <- unique(wideDum$keyVar)
  } else {
    colPal <- grDevices::gray.colors(n = length(levels(wideDum$keyVar)),
                                     start = 0.9, end = 0.05)
    names(colPal) <- levels(wideDum$keyVar)
  }

  p <- ggplot(wideDum, aes(x = catchVar_avg, y = consVar_avg,
                                    shape = shapeVar, fill = keyVar)) +
    geom_point(size = dotSize) +
    scale_shape_manual(values = shapePalette, name = secLegendLab) +
    scale_fill_manual(values = colPal, name = legendLab) +
    guides(fill = guide_legend(override.aes = list(shape = 21)),
           shape = guide_legend(override.aes = list(fill = "black"))) +
    theme_sleekX() +
    theme(strip.text = element_text(size = axisSize),
          axis.text = element_text(size = 0.9 * axisSize),
          axis.title = element_text(size = axisSize),
          legend.text = element_text(size = 0.9 * legendSize),
          legend.title = element_text(size = legendSize)) +
    labs(x = xLab, y = yLab, title = mainLab)

  if (!is.null(facet)) {
    p <- p +
      facet_wrap(~ facetVar, scales = scaleAxis)
    if (facetLetter == TRUE) {
      nLetters <- length(unique(wideDum$facetVar))
      if (scaleAxis %in% c("fixed", "free_y")) {
        labDat <- data.frame(facetVar = unique(wideDum$facetVar),
                             lab = paste(letters[1:nLetters], ")", sep = ""),
                             maxX = ifelse(showUncertainty == TRUE,
                                           max(wideDum$catchVar_highQ),
                                           max(wideDum$catchVar_avg)))
      } else {
        maxX <- wideDum %>%
          group_by(facetVar) %>%
          summarize(maxX = ifelse(showUncertainty == TRUE,
                                  max(catchVar_highQ),
                                  max(catchVar_avg)))
        labDat <- data.frame(facetVar = unique(wideDum$facetVar),
                             lab = paste(letters[1:nLetters], ")", sep = "")) %>%
          left_join(., maxX)
      }
      p <- p +
        geom_text(data = labDat,
                  mapping = aes(x = 0.95 * maxX, y = Inf, label = lab,
                                vjust = 1.75), size = dotSize,
                  show.legend = FALSE, inherit.aes = FALSE)
    } #end if(facetLett == TRUE)
  }
  if (length(unique(wideDum$shapeVar)) < 2) {
    p <- p +
      guides(shape = "none")
  } else {
    p <- p +
      guides(shape = guide_legend(override.aes = list(fill = "black")))
  }
  if (showUncertainty == FALSE) {
    return(p)
  }
  if (showUncertainty==TRUE) {
    q <- p +
      geom_errorbar(aes(ymin = consVar_lowQ, ymax = consVar_highQ),
                    alpha = 0.3, width = 0, size = lineSize) +
      geom_errorbarh(aes(xmin = catchVar_lowQ, xmax = catchVar_highQ),
                     alpha = 0.3, height = 0, size = lineSize)
    return(q)
  }
}