Add basic functionality for pointinterval plots from hubverse quantile tables

DESCRIPTION

@@ -25,6 +25,7 @@ Imports:
     dplyr,
     fs,
     gert,
+    ggdist,
     ggplot2,
     glue,
     gridExtra,

NAMESPACE

@@ -16,7 +16,9 @@ export(get_hubverse_table)
 export(inferencedata_to_tidy_draws)
 export(location_lookup)
 export(nhsn_soda_query)
+export(pivot_hubverse_quantiles_wider)
 export(plot_forecast_quantiles)
+export(plot_hubverse_pointintervals)
 export(plot_hubverse_quantiles)
 export(plot_hubverse_quantiles_loc)
 export(plots_to_pdf)
@@ -25,6 +27,7 @@ export(rank_sampled_trajectories)
 export(sample_aggregated_trajectories)
 export(soql_is_in)
 export(target_end_dates_from_horizons)
+export(theme_forecasttools)
 export(to_location_table_column)
 export(trajectories_to_quantiles)
 export(update_hub)

R/pivot_hubverse_quantiles_wider.R

@@ -0,0 +1,70 @@
+#' Pivot a hubverse quantile table wider with columns representing
+#' individual quantile levels.
+#'
+#' @param hubverse_table hubverse-format forecast table to pivot,
+#' as a [`tibble`][tibble::tibble()]
+#' @param pivot_quantiles quantiles to pivot to columns, as a vector or
+#' named vector. Default `c("point" = 0.5, "lower" = 0.025, "upper" = 0.975)`,
+#' i.e. get the median and the central 95% interquantile interval,
+#' with names "point" for the median, "lower" for the 0.025th quantile,
+#' and "upper" for the 0.975th quantile.
+#' @return A pivoted version of the hubverse table in which each
+#' forecast for a given target, horizon, reference_date, and
+#' location corresponds to a single row with multiple
+#' value columns, one for each of the quantiles in `pivot_quantiles`,
+#' and named according to the corresponding names given in that vector,
+#' or generically as `q_<quantile_level>` if an unnamed numeric
+#' vector is provided.
+#' So with the default `pivot_quantiles`, the output will have three
+#' value columns named `lower"`, `"point"`, and `"upper"`
+#' @export
+pivot_hubverse_quantiles_wider <- function(hubverse_table,
+                                           pivot_quantiles = c(
+                                             "point" = 0.5,
+                                             "lower" = 0.025,
+                                             "upper" = 0.975
+                                           )) {
+  if (!("quantile" %in% hubverse_table$output_type)) {
+    cli::cli_abort(message = paste0(
+      "Hubverse table must contain at least ",
+      "one quantile forecast."
+    ))
+  }
+
+  dat <- hubverse_table |>
+    dplyr::filter(.data$output_type == "quantile") |>
+    dplyr::mutate("output_type_id" = as.numeric(.data$output_type_id))
+
+  pivot_quantiles_present <- pivot_quantiles %in% hubverse_table$output_type_id
+
+  if (!all(pivot_quantiles_present)) {
+    missing_pivot_quantiles <- pivot_quantiles[!pivot_quantiles_present]
+    cli::cli_abort(message = paste0(
+      "Hubverse table is missing one or more of ",
+      "the requested pivot quantiles for all forecasts. ",
+      "The following requested pivot quantiles ",
+      "could not be found: {missing_pivot_quantiles}."
+    ))
+  }
+
+  if (is.null(names(pivot_quantiles))) {
+    names(pivot_quantiles) <- paste("q", pivot_quantiles, sep = "")
+  }
+
+  pivot_quant_map <- setNames(names(pivot_quantiles), pivot_quantiles)
+
+  dat <- dat |>
+    dplyr::filter(.data$output_type_id %in% !!pivot_quantiles) |>
+    dplyr::mutate(
+      "which_quantile" = dplyr::recode(
+        .data$output_type_id,
+        !!!pivot_quant_map
+      )
+    ) |>
+    dplyr::select(-"output_type", -"output_type_id") |>
+    tidyr::pivot_wider(
+      names_from = "which_quantile",
+      values_from = "value"
+    )
+  return(dat)
+}

R/plot_hubverse_pointinterval.R

@@ -0,0 +1,84 @@
+#' Plot a table of hubverse-formatted forecasts
+#' as pointintervals
+#'
+#' @param hubverse_table Hubverse table, as a [`tibble`][tibble::tibble()]
+#' @param horizons horizons to plot, as a vector of integers. If `NULL`,
+#' plot all available, each in its own facet. Default `NULL`.
+#' @param point_estimate_quantile Quantile to plot as the point estimate.
+#' Default `0.5`, the median.
+#' @param lower_limit_quantile Quantile to plot as the lower bound of
+#' the interval. Default `0.025`.
+#' @param upper_limit_quantile Quantile to plot as the upper bound of
+#' the interval. Default `0.975`.
+#' @param location_input_format Format of the hubverse table
+#' `location` column.
+#' Permitted formats are `"abbr"` (state/territory
+#' or nation two letter USPS abbreviation), `"hub"` (
+#' legacy 2-digit FIPS code for states and territories, `US`
+#' for the USA as a whole), and `"long_name"` (full English
+#' jurisdiction names; not recommended). Default `"hub"`.
+#' @param location_output_format How to code locations in the output plot.
+#' Permitted formats are `"abbr"` (state/territory
+#' or nation two letter USPS abbreviation), `"hub"` (
+#' legacy 2-digit FIPS code for states and territories, `US`
+#' #' for the USA as a whole), and `"long_name"` (full English
+#' jurisdiction names). Default `"abbr"`.
+#' @return A ggplot2 plot of the forecasts as pointintervals
+#' @export
+plot_hubverse_pointintervals <- function(hubverse_table,
+                                         horizons = NULL,
+                                         point_estimate_quantile = 0.5,
+                                         lower_limit_quantile = 0.025,
+                                         upper_limit_quantile = 0.975,
+                                         location_input_format = "hub",
+                                         location_output_format = "abbr") {
+  if (is.null(horizons)) {
+    horizons <- unique(hubverse_table$horizon)
+  }
+
+  pivoted <- hubverse_table |>
+    pivot_hubverse_quantiles_wider(
+      pivot_quantiles = c(
+        "point" = point_estimate_quantile,
+        "lower" = lower_limit_quantile,
+        "upper" = upper_limit_quantile
+      )
+    ) |>
+    dplyr::mutate(
+      "location" = location_lookup(
+        .data$location,
+        location_input_format,
+        location_output_format
+      )
+    ) |>
+    dplyr::filter(.data$horizon %in% horizons)
+
+  ## order by the median at the first horizon plotted
+  loc_levels <- pivoted |>
+    dplyr::filter(
+      .data$horizon == min(.data$horizon)
+    ) |>
+    dplyr::arrange(.data$point) |>
+    dplyr::pull("location")
+
+  pivoted <- pivoted |>
+    dplyr::mutate("location" = factor(
+      location,
+      levels = loc_levels,
+      ordered = TRUE
+    ))
+
+  plot <- pivoted |>
+    ggplot2::ggplot(ggplot2::aes_string(
+      y = "location",
+      x = "point",
+      xmin = "lower",
+      xmax = "upper"
+    )) +
+    ggdist::geom_pointinterval() +
+    ggplot2::facet_wrap(~horizon ~ target) +
+    ggplot2::labs(x = "target value") +
+    theme_forecasttools()
+
+  return(plot)
+}

R/plot_quantiles.R

@@ -169,12 +169,12 @@ plot_hubverse_quantiles_loc <- function(location,
       ggplot2::scale_y_continuous(
         trans = ytrans
       ) +
-      ggplot2::theme_minimal(base_size = 15) +
       ggplot2::labs(
         y = target_name,
         x = "Date"
       ) +
-      ggplot2::ggtitle(plot_title)
+      ggplot2::ggtitle(plot_title) +
+      theme_forecasttools()
 
   )
 

R/recode_locations.R

@@ -83,29 +83,44 @@ to_location_table_column <- function(location_format) {
 #' corresponding to a given location vector
 #' and format, with repeats possible
 #' @param location_vector vector of location values
-#' @param location_format how the vector is coded.
-#' @param location_format format in which the location
+#' @param location_input_format format in which the location
 #' vector is coded.
 #' Permitted formats are `"abbr"` (state/territory
 #' or nation two letter USPS abbreviation), `"hub"`
 #' (legacy 2-digit FIPS code for states and territories, `US`
 #' for the USA as a whole), and `"long_name"` (full English
 #' name for jurisdiction).
-#' @return the rows of the [us_location_table]
+#' @param location_output_format Return only this column of the
+#' output table, if it is provided. Otherwise return the whole
+#' table. Default `NULL` (return all columns).
+#  Permitted formats are `"abbr"` (state/territory
+#' or nation two letter USPS abbreviation), `"hub"`
+#' (legacy 2-digit FIPS code for states and territories, `US`
+#' for the USA as a whole), and `"long_name"` (full English
+#' name for jurisdiction).
+#' @return the corresponding rows of the [us_location_table]
 #' matching the location vector, with repeats possible.
 #' @export
 location_lookup <- function(location_vector,
-                            location_format) {
+                            location_input_format,
+                            location_output_format = NULL) {
   ## coerce location vector to character
   ## (handles case of just states with no US,
   ## by FIPS, which R will default to treating
   ## as numeric)
   location_vector <- as.character(location_vector)
-  join_key <- to_location_table_column(location_format)
+  join_key <- to_location_table_column(location_input_format)
   locs <- tibble::tibble({{ join_key }} := !!location_vector) |>
     dplyr::inner_join(forecasttools::us_location_table,
       by = join_key
     )
 
-  return(locs)
+  if (is.null(location_output_format)) {
+    result <- locs
+  } else {
+    output_column <- to_location_table_column(location_output_format)
+    result <- locs[[output_column]]
+  }
+
+  return(result)
 }

R/theme_forecasttools.R

@@ -0,0 +1,13 @@
+#' forecasttools ggplot2 theme.
+#'
+#' Built on top of [ggplot2::theme_minimal()]
+#' @param ... arguments passed to [ggplot2::theme()]
+#' @return The theme, which can be added to a ggplot
+#' object.
+#' @export
+theme_forecasttools <- function(...) {
+  return(
+    ggplot2::theme_minimal(base_size = 15) +
+      ggplot2::theme(...)
+  )
+}