R/plot_density.R

In DUToolkit: Visualizing and Quantifying Decision Uncertainty

#' Generate density plots
#'
#' @description For each policy alternative, this function generates
#' probability density plots of the highest (or lowest if the threshold is a
#' minimum) projected outcome across simulation runs. The decision threshold
#' is shown directly on the plot as a vertical line. The area under the
#' probability density curve where the threshold value is exceeded is shaded
#' to visually display the downside risk of the policy alternative.
#'
#' @param max_min_values_list A list generated by [get_max_min_values()]
#' @param D A single threshold value
#' @param Dt_max A logical value indicating whether the decision threshold
#' is a maximum (`TRUE`) or a minimum (`FALSE`). The default is `TRUE`.
#' @param risk_measures A list of risk scores generated by [calculate_risk()].
#' The policy alternatives in the `risk_measures` list must be in the same order
#' as in the `max_min_values_list`.
#'
#' @return A list of ggplots, one for each policy alternative.
#' @export
#'
#' @examples
#' tmin <- "2021-01-01"
#' tmax <- "2021-04-10"
#' Dt <- rep(750, nrow(psa_data$Baseline))
#' D <- 750
#'
#' risk_measures <- calculate_risk(
#'   psa_data,
#'   tmin = tmin,
#'   tmax = tmax,
#'   Dt = Dt,
#'   Dt_max = TRUE
#' )
#'
#' peak_values_list <- get_max_min_values(
#'   psa_data,
#'   tmin = tmin,
#'   tmax = tmax,
#'   Dt_max = TRUE
#' )
#'
#' density_plots <- plot_density(
#'   peak_values_list,
#'   D = D,
#'   Dt_max = TRUE,
#'   risk_measures
#' )
plot_density <- function(
    max_min_values_list,
    D,
    Dt_max = TRUE,
    risk_measures) {
  if (inherits(max_min_values_list, "list")) {
    # create plots
    plot <- lapply(max_min_values_list, plot_density_1, D, Dt_max)
    # create captions
    captions <- caption_density_plot(max_min_values_list, D, Dt_max, risk_measures)
    if (!is.null(names(plot)) && all(names(plot) != "")) {
      # add captions and titles to plots
      plot <- stats::setNames(
        lapply(names(plot), function(name) {
          plot[[name]] + ggplot2::ggtitle(name) +
            ggplot2::labs(subtitle = captions[[name]]) +
            ggplot2::theme(plot.subtitle = ggtext::element_textbox_simple(size = 9))
        }),
        names(plot)
      )
    } else {
      # Assign sequential names if unnamed
      plot_names <- paste0("Plot_", seq_along(plot))
      plot <- stats::setNames(
        lapply(seq_along(plot), function(i) {
          plot[[i]] + ggplot2::ggtitle(plot_names[i]) +
            ggplot2::labs(subtitle = captions[[i]]) +
            ggplot2::theme(plot.subtitle = ggtext::element_textbox_simple(size = 9))
        }),
        plot_names
      )
    }
  } else if (inherits(max_min_values_list, "data.frame")) {
    plot <- plot_density_1(max_min_values_list, D, Dt_max)
  } else {
    rlang::abort("The first argument is not a data.frame or list of data.frames",
      class = "data_type"
    )
  }
  return(plot)
}


#' Plot density plot for a single set of simulations
#'
#' @param max_min_values_list A list (or data.frame) generated by get_max_min_values()
#' @param D A single threshold value
#' @param Dt_max Logical indicating whether decision threshold is a
#' maximum (`TRUE`) or minimum (`FALSE`)
#' @noRd
#' @return ggplot density plot
plot_density_1 <- function(
    max_min_values_list,
    D,
    Dt_max = TRUE) {
  shading_data <- data.frame(
    x = stats::density(max_min_values_list$outcome)$x,
    y = stats::density(max_min_values_list$outcome)$y
  )

  density_plot <- ggplot2::ggplot(
    data = max_min_values_list,
    ggplot2::aes(x = !!rlang::sym("outcome"))
  )

  if (Dt_max) { # if threshold is a maximum
    density_plot <- density_plot +
      ggplot2::geom_ribbon(
        data = shading_data[shading_data$x >= D, ],
        # data = subset(shading_data, !!rlang::sym("x") >= D),
        ggplot2::aes(
          x = !!rlang::sym("x"), ymin = 0,
          ymax = !!rlang::sym("y"),
          fill = "Downside risk"
        ),
        alpha = 0.5
      )
  } else {
    density_plot <- density_plot +
      ggplot2::geom_ribbon(
        data = shading_data[shading_data$x <= D, ],
        # data = subset(shading_data, !!rlang::sym("x") <= D),
        ggplot2::aes(
          x = !!rlang::sym("x"),
          ymin = 0,
          ymax = !!rlang::sym("y"),
          fill = "Downside risk"
        ),
        alpha = 0.5
      )
  }

  density_plot <- density_plot +
    ggplot2::scale_color_manual(values = c("Policy target" = "red")) +
    ggplot2::scale_fill_manual(values = c("Downside risk" = "#9386A6FF"))

  density_plot <- density_plot + ggplot2::geom_density() +
    ggplot2::theme_classic() +
    ggplot2::labs(x = "Outcome at peak", y = "Probability")

  density_plot <- density_plot +
    ggplot2::geom_vline(ggplot2::aes(xintercept = D, color = "Policy target"),
      linetype = "dashed", linewidth = 0.75
    )

  density_plot <- density_plot +
    ggplot2::guides(
      fill = ggplot2::guide_legend(override.aes = list(
        linetype = 0,
        alpha = 0.5, order = 2
      )),
      color = ggplot2::guide_legend(order = 1)
    )


  density_plot <- density_plot + ggplot2::theme(
    legend.position = "bottom",
    legend.title = ggplot2::element_blank(),
    legend.text = ggplot2::element_text(size = 8)
  )

  density_plot <- gen_stand_descr(density_plot,
                                  link = "Probability density plots with risk shading")

  return(density_plot)
}


#' Create caption for density plots
#'
#' @inheritParams plot_density
#' @noRd
#' @return list of text captions (or single caption)
caption_density_plot <- function(max_min_values_list, D, Dt_max = TRUE, risk_measures) {
  if (inherits(max_min_values_list, "list")) {
    caption_all <- vector(mode = "list", length = length(max_min_values_list))
    names(caption_all) <- names(max_min_values_list)

    # calculate RR for risk measure
    rr_all <- calculate_rr(risk_measures)
    # calculate probability of exceeding threshold
    prob <- calculate_threshold_probs(max_min_values_list, Dp = D, Dt_max)
    # calculate risk values at the peak
    risk_peak <- calculate_max_min_risk(max_min_values_list, D, Dt_max)
    # calculate RR for risk measures at peak
    rr_peak <- calculate_rr(risk_peak)

    for (i in 1:length(max_min_values_list)) {
      if (i == 1) {
        caption2 <- NULL
        caption3 <- NULL
      } else {
        caption2 <- paste(
          scales::label_percent(accuracy = 1)(abs(rr_peak[i - 1])),
          ifelse(rr_peak[i - 1] <= 0, "reduction", "increase"),
          "in the risk of",
          ifelse(Dt_max == TRUE, "exceeding", "falling below"),
          "the threshold relative to the baseline scenario at the forecasted",
          ifelse(Dt_max == TRUE, "peak", "low point")
        )
        caption3 <- paste(
          scales::label_percent(accuracy = 1)(abs(rr_all[i - 1])),
          ifelse(rr_all[i - 1] <= 0, "reduction", "increase"),
          "in the risk of",
          ifelse(Dt_max == TRUE, "exceeding", "falling below"),
          "the threshold relative to the baseline scenario for the
                        full time span<br>"
        )
      }

      caption1 <- paste(
        scales::label_percent(accuracy = 1)(prob[[i]]), "probability of",
        ifelse(Dt_max == TRUE, "exceeding", "falling below"),
        "threshold at the forecasted",
        ifelse(Dt_max == TRUE, "peak", "low point")
      )

      caption_all[[i]] <- paste(caption1, caption2, caption3, sep = "<br>")
    }
  } else if (inherits(max_min_values_list, "data.frame")) {
    prob <- calculate_threshold_probs(max_min_values_list, Dp = D, Dt_max)

    caption_all <- paste(
      scales::label_percent(accuracy = 1)(prob), "probability of",
      ifelse(Dt_max == TRUE, "exceeding", "falling below"),
      "threshold at the forecasted",
      ifelse(Dt_max == TRUE, "peak", "low point")
    )
  }

  return(caption_all)
}