diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
index 7dc0f9d8..c0683f98 100644
--- a/.github/workflows/R-CMD-check.yaml
+++ b/.github/workflows/R-CMD-check.yaml
@@ -17,7 +17,7 @@ jobs:
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
       R_KEEP_PKG_SOURCE: yes
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
 
       - uses: r-lib/actions/setup-r@v2
         with:
diff --git a/.github/workflows/doc-preview.yaml b/.github/workflows/doc-preview.yaml
index 06818422..22d487f0 100644
--- a/.github/workflows/doc-preview.yaml
+++ b/.github/workflows/doc-preview.yaml
@@ -17,7 +17,7 @@ jobs:
       # Needed to read the PR branch
       contents: read
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
         with:
           # Checkout the PR branch
           ref: refs/pull/${{ github.event.issue.number }}/head
diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml
index 2aae89fa..c8462eea 100644
--- a/.github/workflows/pkgdown.yaml
+++ b/.github/workflows/pkgdown.yaml
@@ -26,7 +26,7 @@ jobs:
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
       DELPHI_EPIDATA_KEY: ${{ secrets.SECRET_EPIPREDICT_GHACTIONS_DELPHI_EPIDATA_KEY }}
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
 
       - uses: r-lib/actions/setup-pandoc@v2
 
diff --git a/.github/workflows/pr-commands.yaml b/.github/workflows/pr-commands.yaml
index 8bb12b9e..5c4ea908 100644
--- a/.github/workflows/pr-commands.yaml
+++ b/.github/workflows/pr-commands.yaml
@@ -57,7 +57,7 @@ jobs:
     permissions:
       contents: write
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
 
       - uses: r-lib/actions/pr-fetch@v2
         with:
diff --git a/DESCRIPTION b/DESCRIPTION
index 2166122a..8aa08d33 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: epipredict
 Title: Basic epidemiology forecasting methods
-Version: 0.2.6
+Version: 0.2.7
 Authors@R: c(
     person("Daniel J.", "McDonald", , "daniel@stat.ubc.ca", role = c("aut", "cre")),
     person("Ryan", "Tibshirani", , "ryantibs@cmu.edu", role = "aut"),
@@ -21,7 +21,7 @@ Description: A forecasting "framework" for creating epidemiological
     are included.
 License: MIT + file LICENSE
 URL: https://github.com/cmu-delphi/epipredict/,
-    https://cmu-delphi.github.io/epipredict
+    https://cmu-delphi.github.io/epipredict/
 BugReports: https://github.com/cmu-delphi/epipredict/issues/
 Depends:
     epidatasets,
@@ -31,7 +31,7 @@ Imports:
     checkmate,
     cli,
     dplyr,
-    epiprocess (>= 0.11.2),
+    epiprocess (>= 0.13.0),
     generics,
     ggplot2,
     glue,
diff --git a/NEWS.md b/NEWS.md
index 4b33df25..77b60459 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,11 +2,21 @@
 
 Pre-1.0.0 numbering scheme: 0.x will indicate releases, while 0.0.x will indicate PR's.
 
+# epipredict 0.2.7
+
+- `autoplot.canned_epipred()` had internal code using
+  `epiprocess:::autoplot.epi_df(.max_facets = Inf)`. This argument no longer
+  exists there, so it caused an error (noted by rebuilding the `README.Rmd`)
+
+
 # epipredict 0.2.6
 
 - `arx_forecaster()` and `flatline_forecaster()` now error early when `quantile_by_key` contains columns that are not keys of the input `epi_df`, rather than silently dropping the invalid keys (#229).
 - `arx_forecaster()` now warns when `quantile_by_key` is supplied with a quantile-output trainer (`quantile_reg()`, `rand_forest()` with engine `"grf_quantiles"`), where the argument would otherwise be silently ignored (#229).
-- Regenerate roxygen-derived `man/step_adjust_latency.Rd` so its recorded example output matches the current `epi_df` print phrasing (`lag` → `latency`) from upstream `epiprocess`. No user-visible behavior change.
+- Regenerate roxygen-derived `man/step_adjust_latency.Rd` so its recorded
+  example output matches the current `epi_df` print phrasing (`lag` → `latency`)
+  from upstream `epiprocess`. No user-visible behavior change.
+- Bump GitHub action checkout version.
 
 # epipredict 0.2.5
 
diff --git a/R/arx_classifier.R b/R/arx_classifier.R
index 57812233..44acca28 100644
--- a/R/arx_classifier.R
+++ b/R/arx_classifier.R
@@ -141,21 +141,33 @@
 #'   )
 #' )
 arx_classifier <- function(
-    epi_data,
-    outcome,
-    predictors,
-    trainer = logistic_reg(),
-    args_list = arx_class_args_list()) {
+  epi_data,
+  outcome,
+  predictors,
+  trainer = logistic_reg(),
+  args_list = arx_class_args_list()
+) {
   if (!is_classification(trainer)) {
-    cli_abort("`trainer` must be a {.pkg parsnip} model of mode 'classification'.")
+    cli_abort(
+      "`trainer` must be a {.pkg parsnip} model of mode 'classification'."
+    )
   }
 
-  wf <- arx_class_epi_workflow(epi_data, outcome, predictors, trainer, args_list)
+  wf <- arx_class_epi_workflow(
+    epi_data,
+    outcome,
+    predictors,
+    trainer,
+    args_list
+  )
   wf <- fit(wf, epi_data)
 
   if (args_list$adjust_latency == "none") {
     forecast_date_default <- max(epi_data$time_value)
-    if (!is.null(args_list$forecast_date) && args_list$forecast_date != forecast_date_default) {
+    if (
+      !is.null(args_list$forecast_date) &&
+        args_list$forecast_date != forecast_date_default
+    ) {
       cli_warn(
         "The specified forecast date {args_list$forecast_date} doesn't match the
         date from which the forecast is occurring {forecast_date}."
@@ -218,23 +230,31 @@ arx_classifier <- function(
 #'   )
 #' )
 arx_class_epi_workflow <- function(
-    epi_data,
-    outcome,
-    predictors,
-    trainer = parsnip::logistic_reg(),
-    args_list = arx_class_args_list()) {
+  epi_data,
+  outcome,
+  predictors,
+  trainer = parsnip::logistic_reg(),
+  args_list = arx_class_args_list()
+) {
   validate_forecaster_inputs(epi_data, outcome, predictors)
   if (!inherits(args_list, c("arx_class", "alist"))) {
     cli_abort("`args_list` was not created using `arx_class_args_list()`.")
   }
   if (!(is.null(trainer) || is_classification(trainer))) {
-    cli_abort("`trainer` must be a {.pkg parsnip} model of mode 'classification'.")
+    cli_abort(
+      "`trainer` must be a {.pkg parsnip} model of mode 'classification'."
+    )
   }
 
   if (args_list$adjust_latency == "none") {
     forecast_date_default <- max(epi_data$time_value)
-    if (!is.null(args_list$forecast_date) && args_list$forecast_date != forecast_date_default) {
-      cli_warn("The specified forecast date {args_list$forecast_date} doesn't match the date from which the forecast is occurring {forecast_date}.")
+    if (
+      !is.null(args_list$forecast_date) &&
+        args_list$forecast_date != forecast_date_default
+    ) {
+      cli_warn(
+        "The specified forecast date {args_list$forecast_date} doesn't match the date from which the forecast is occurring {forecast_date}."
+      )
     }
   } else {
     forecast_date_default <- attributes(epi_data)$metadata$as_of
@@ -257,7 +277,8 @@ arx_class_epi_workflow <- function(
   for (l in seq_along(lags)) {
     pred_names <- predictors[l]
     pred_names <- as.character(glue::glue_data(
-      args_list, "gr_{horizon}_{method}_{pred_names}"
+      args_list,
+      "gr_{horizon}_{method}_{pred_names}"
     ))
     r <- step_epi_lag(r, !!pred_names, lag = lags[[l]])
   }
@@ -293,17 +314,24 @@ arx_class_epi_workflow <- function(
   method_adjust_latency <- args_list$adjust_latency
   if (method_adjust_latency != "none") {
     if (method_adjust_latency != "extend_ahead") {
-      cli_abort("only extend_ahead is currently supported",
+      cli_abort(
+        "only extend_ahead is currently supported",
         class = "epipredict__arx_classifier__adjust_latency_unsupported_method"
       )
     }
-    r <- r %>% step_adjust_latency(!!pre_out_name,
-      fixed_forecast_date = forecast_date,
-      method = method_adjust_latency
-    )
+    r <- r %>%
+      step_adjust_latency(
+        !!pre_out_name,
+        fixed_forecast_date = forecast_date,
+        method = method_adjust_latency
+      )
   }
   r <- r %>%
-    step_epi_ahead(!!pre_out_name, ahead = args_list$ahead, role = "pre-outcome")
+    step_epi_ahead(
+      !!pre_out_name,
+      ahead = args_list$ahead,
+      role = "pre-outcome"
+    )
   r <- r %>%
     step_mutate(
       across(
@@ -345,8 +373,8 @@ arx_class_epi_workflow <- function(
 #'   be created using growth rates (as the predictors are) or lagged
 #'   differences. The second case is closer to the requirements for the
 #'   [2022-23 CDC Flusight Hospitalization Experimental Target](https://github.com/cdcepi/Flusight-forecast-data/blob/745511c436923e1dc201dea0f4181f21a8217b52/data-experimental/README.md).
-#'   See the [Classification chapter from the forecasting book](https://cmu-delphi.github.io/delphi-tooling-book/arx-classifier.html) Vignette for details of how to create a reasonable
-#'   baseline for this case. Selecting `"growth_rate"` (the default) uses
+#'   See the `vignette("epipredict")` for more details.
+#'   Selecting `"growth_rate"` (the default) uses
 #'   [epiprocess::growth_rate()] to create the outcome using some of the
 #'   additional arguments below. Choosing `"lag_difference"` instead simply
 #'   uses the change from the value at the selected `horizon`.
@@ -384,36 +412,48 @@ arx_class_epi_workflow <- function(
 #' # also needs arx_classifier(trainer = parsnip::multinom_reg())
 #' arx_class_args_list(breaks = c(-.2, .25))
 arx_class_args_list <- function(
-    lags = c(0L, 7L, 14L),
-    ahead = 7L,
-    n_training = Inf,
-    forecast_date = NULL,
-    target_date = NULL,
-    adjust_latency = c("none", "extend_ahead", "extend_lags", "locf"),
-    warn_latency = TRUE,
-    outcome_transform = c("growth_rate", "lag_difference"),
-    breaks = 0.25,
-    horizon = 7L,
-    method = c("rel_change", "linear_reg"),
-    log_scale = FALSE,
-    check_enough_data_n = NULL,
-    check_enough_data_epi_keys = NULL,
-    ...) {
+  lags = c(0L, 7L, 14L),
+  ahead = 7L,
+  n_training = Inf,
+  forecast_date = NULL,
+  target_date = NULL,
+  adjust_latency = c("none", "extend_ahead", "extend_lags", "locf"),
+  warn_latency = TRUE,
+  outcome_transform = c("growth_rate", "lag_difference"),
+  breaks = 0.25,
+  horizon = 7L,
+  method = c("rel_change", "linear_reg"),
+  log_scale = FALSE,
+  check_enough_data_n = NULL,
+  check_enough_data_epi_keys = NULL,
+  ...
+) {
   rlang::check_dots_empty()
   .lags <- lags
-  if (is.list(lags)) lags <- unlist(lags)
+  if (is.list(lags)) {
+    lags <- unlist(lags)
+  }
   method <- rlang::arg_match(method)
   outcome_transform <- rlang::arg_match(outcome_transform)
 
   adjust_latency <- rlang::arg_match(adjust_latency)
-  arg_is_scalar(ahead, n_training, horizon, log_scale, adjust_latency, warn_latency)
+  arg_is_scalar(
+    ahead,
+    n_training,
+    horizon,
+    log_scale,
+    adjust_latency,
+    warn_latency
+  )
   arg_is_scalar(forecast_date, target_date, allow_null = TRUE)
   arg_is_date(forecast_date, target_date, allow_null = TRUE)
   arg_is_nonneg_int(ahead, lags, horizon)
   arg_is_numeric(breaks)
   arg_is_lgl(log_scale)
   arg_is_pos(n_training)
-  if (is.finite(n_training)) arg_is_pos_int(n_training)
+  if (is.finite(n_training)) {
+    arg_is_pos_int(n_training)
+  }
   arg_is_pos(check_enough_data_n, allow_null = TRUE)
   arg_is_chr(check_enough_data_epi_keys, allow_null = TRUE)
 
@@ -428,9 +468,12 @@ arx_class_args_list <- function(
   }
 
   breaks <- sort(breaks)
-  if (min(breaks) > -Inf) breaks <- c(-Inf, breaks)
-  if (max(breaks) < Inf) breaks <- c(breaks, Inf)
-
+  if (min(breaks) > -Inf) {
+    breaks <- c(-Inf, breaks)
+  }
+  if (max(breaks) < Inf) {
+    breaks <- c(breaks, Inf)
+  }
 
   max_lags <- max(lags)
   structure(
diff --git a/R/autoplot.R b/R/autoplot.R
index 65db2d5c..08184e28 100644
--- a/R/autoplot.R
+++ b/R/autoplot.R
@@ -86,9 +86,24 @@ autoplot.epi_workflow <- function(
   object,
   predictions = NULL,
   observed_response = NULL,
-  .levels = c(.5, .8, .9), ...,
-  .color_by = c("all_keys", "geo_value", "other_keys", ".response", "all", "none"),
-  .facet_by = c(".response", "other_keys", "all_keys", "geo_value", "all", "none"),
+  .levels = c(.5, .8, .9),
+  ...,
+  .color_by = c(
+    "all_keys",
+    "geo_value",
+    "other_keys",
+    ".response",
+    "all",
+    "none"
+  ),
+  .facet_by = c(
+    ".response",
+    "other_keys",
+    "all_keys",
+    "geo_value",
+    "all",
+    "none"
+  ),
   .base_color = "dodgerblue4",
   .point_pred_color = "orange",
   .facet_filter = NULL
@@ -114,7 +129,9 @@ autoplot.epi_workflow <- function(
   keys <- c("geo_value", "time_value", "key")
   mold_roles <- names(mold$extras$roles)
   # extract the relevant column names for plotting
-  if (starts_with_impl("ahead_", names(y)) || starts_with_impl("lag_", names(y))) {
+  if (
+    starts_with_impl("ahead_", names(y)) || starts_with_impl("lag_", names(y))
+  ) {
     old_name_y <- unlist(strsplit(names(y), "_"))
     new_name_y <- paste(old_name_y[-c(1:2)], collapse = "_")
   } else {
@@ -135,18 +152,25 @@ autoplot.epi_workflow <- function(
     }
     observed_response <- rename(observed_response, !!new_name_y := !!names(y))
     if (!is.null(shift)) {
-      observed_response <- mutate(observed_response, time_value = time_value + shift)
+      observed_response <- mutate(
+        observed_response,
+        time_value = time_value + shift
+      )
     }
     other_keys <- setdiff(key_colnames(object), c("geo_value", "time_value"))
-    observed_response <- as_epi_df(observed_response,
+    observed_response <- as_epi_df(
+      observed_response,
       as_of = object$fit$meta$as_of,
       other_keys = other_keys
     )
   }
   if (is.null(predictions)) {
     return(autoplot(
-      observed_response, new_name_y,
-      .color_by = .color_by, .facet_by = .facet_by, .base_color = .base_color,
+      observed_response,
+      new_name_y,
+      .color_by = .color_by,
+      .facet_by = .facet_by,
+      .base_color = .base_color,
       .facet_filter = {{ .facet_filter }},
       # Avoid subsampling while it is not implemented for this method
       .max_keys = Inf
@@ -159,15 +183,21 @@ autoplot.epi_workflow <- function(
     }
     predictions <- rename(predictions, time_value = target_date)
   }
-  pred_cols_ok <- hardhat::check_column_names(predictions, key_colnames(observed_response))
+  pred_cols_ok <- hardhat::check_column_names(
+    predictions,
+    key_colnames(observed_response)
+  )
   if (!pred_cols_ok$ok) {
     cli_warn(c(
       "`predictions` is missing required variables: {.var {pred_cols_ok$missing_names}}.",
       i = "Plotting the original data."
     ))
     return(autoplot(
-      observed_response, !!new_name_y,
-      .color_by = .color_by, .facet_by = .facet_by, .base_color = .base_color,
+      observed_response,
+      !!new_name_y,
+      .color_by = .color_by,
+      .facet_by = .facet_by,
+      .base_color = .base_color,
       .facet_filter = {{ .facet_filter }},
       # Avoid subsampling while it is not implemented for this method
       .max_keys = Inf
@@ -175,9 +205,13 @@ autoplot.epi_workflow <- function(
   }
 
   # First we plot the history, always faceted by everything
-  bp <- autoplot(observed_response, !!new_name_y,
-    .color_by = "none", .facet_by = "all_keys",
-    .base_color = "black", .facet_filter = {{ .facet_filter }},
+  bp <- autoplot(
+    observed_response,
+    !!new_name_y,
+    .color_by = "none",
+    .facet_by = "all_keys",
+    .base_color = "black",
+    .facet_filter = {{ .facet_filter }},
     # Avoid subsampling while it is not implemented for this method
     .max_keys = Inf
   )
@@ -194,7 +228,6 @@ autoplot.epi_workflow <- function(
       mutate(.facets = droplevels(.facets))
   }
 
-
   if (".pred_distn" %in% names(predictions)) {
     bp <- plot_bands(bp, predictions, .levels, .base_color)
   }
@@ -225,9 +258,25 @@ autoplot.epi_workflow <- function(
 #' @export
 #' @rdname autoplot-epipred
 autoplot.canned_epipred <- function(
-  object, observed_response = NULL, ...,
-  .color_by = c("all_keys", "geo_value", "other_keys", ".response", "all", "none"),
-  .facet_by = c(".response", "other_keys", "all_keys", "geo_value", "all", "none"),
+  object,
+  observed_response = NULL,
+  ...,
+  .color_by = c(
+    "all_keys",
+    "geo_value",
+    "other_keys",
+    ".response",
+    "all",
+    "none"
+  ),
+  .facet_by = c(
+    ".response",
+    "other_keys",
+    "all_keys",
+    "geo_value",
+    "all",
+    "none"
+  ),
   .base_color = "dodgerblue4",
   .point_pred_color = "orange",
   .facet_filter = NULL
@@ -240,9 +289,15 @@ autoplot.canned_epipred <- function(
   predictions <- object$predictions %>%
     rename(time_value = target_date)
 
-  autoplot(ewf, predictions, observed_response, ...,
-    .color_by = .color_by, .facet_by = .facet_by,
-    .base_color = .base_color, .facet_filter = {{ .facet_filter }}
+  autoplot(
+    ewf,
+    predictions,
+    observed_response,
+    ...,
+    .color_by = .color_by,
+    .facet_by = .facet_by,
+    .base_color = .base_color,
+    .facet_filter = {{ .facet_filter }}
   )
 }
 
@@ -265,7 +320,8 @@ starts_with_impl <- function(x, vars) {
 }
 
 plot_bands <- function(
-  base_plot, predictions,
+  base_plot,
+  predictions,
   levels = c(.5, .8, .9),
   fill = "blue4",
   alpha = 0.6,
@@ -281,7 +337,9 @@ plot_bands <- function(
   ntarget_dates <- dplyr::n_distinct(predictions$time_value)
 
   predictions <- predictions %>%
-    mutate(.pred_distn = quantile_pred(quantile(.pred_distn, levels), levels)) %>%
+    mutate(
+      .pred_distn = quantile_pred(quantile(.pred_distn, levels), levels)
+    ) %>%
     pivot_quantiles_wider(.pred_distn)
   qnames <- setdiff(names(predictions), innames)
 
@@ -293,17 +351,29 @@ plot_bands <- function(
         base_plot <- base_plot +
           geom_ribbon(
             data = predictions,
-            aes(x = .data$time_value, ymin = .data[[bottom]], ymax = .data[[top]]),
+            aes(
+              x = .data$time_value,
+              ymin = .data[[bottom]],
+              ymax = .data[[top]]
+            ),
             inherit.aes = FALSE,
-            alpha = 0.2, linewidth = linewidth, fill = fill
+            alpha = 0.2,
+            linewidth = linewidth,
+            fill = fill
           )
       } else {
         base_plot <- base_plot +
           geom_linerange(
             data = predictions,
-            aes(x = .data$time_value, ymin = .data[[bottom]], ymax = .data[[top]]),
+            aes(
+              x = .data$time_value,
+              ymin = .data[[bottom]],
+              ymax = .data[[top]]
+            ),
             inherit.aes = FALSE,
-            alpha = 0.2, linewidth = 2, color = fill
+            alpha = 0.2,
+            linewidth = 2,
+            color = fill
           )
       }
     } else {
@@ -311,17 +381,28 @@ plot_bands <- function(
         base_plot <- base_plot +
           geom_ribbon(
             data = predictions,
-            aes(x = .data$time_value, ymin = .data[[bottom]], ymax = .data[[top]]),
+            aes(
+              x = .data$time_value,
+              ymin = .data[[bottom]],
+              ymax = .data[[top]]
+            ),
             inherit.aes = FALSE,
-            fill = fill, alpha = alpha
+            fill = fill,
+            alpha = alpha
           )
       } else {
         base_plot <- base_plot +
           geom_linerange(
             data = predictions,
-            aes(x = .data$time_value, ymin = .data[[bottom]], ymax = .data[[top]]),
+            aes(
+              x = .data$time_value,
+              ymin = .data[[bottom]],
+              ymax = .data[[top]]
+            ),
             inherit.aes = FALSE,
-            color = fill, alpha = alpha, linewidth = 2
+            color = fill,
+            alpha = alpha,
+            linewidth = 2
           )
       }
     }
diff --git a/R/cdc_baseline_forecaster.R b/R/cdc_baseline_forecaster.R
index a97eece8..34f017d6 100644
--- a/R/cdc_baseline_forecaster.R
+++ b/R/cdc_baseline_forecaster.R
@@ -9,8 +9,8 @@
 #' By default, the predictive intervals are computed separately for each
 #' combination of `geo_value` in the `epi_data` argument.
 #'
-#' This forecaster is meant to produce exactly the CDC Baseline used for
-#' [COVID19ForecastHub](https://covid19forecasthub.org)
+#' This forecaster is meant to produce exactly the CDC Baseline used for, e.g.,
+#' <https://github.com/cdcepi/Flusight-baseline/>
 #'
 #' @param epi_data An [`epiprocess::epi_df`][epiprocess::as_epi_df]
 #' @param outcome A scalar character for the column name we wish to predict.
@@ -56,9 +56,10 @@
 #'   theme_bw() +
 #'   geom_vline(xintercept = forecast_date)
 cdc_baseline_forecaster <- function(
-    epi_data,
-    outcome,
-    args_list = cdc_baseline_args_list()) {
+  epi_data,
+  outcome,
+  args_list = cdc_baseline_args_list()
+) {
   validate_forecaster_inputs(epi_data, outcome, "time_value")
   if (!inherits(args_list, c("cdc_flat_fcast", "alist"))) {
     cli_abort("`args_list` was not created using `cdc_baseline_args_list().")
@@ -67,7 +68,6 @@ cdc_baseline_forecaster <- function(
   ek <- kill_time_value(keys)
   outcome <- rlang::sym(outcome)
 
-
   r <- epi_recipe(epi_data) %>%
     step_epi_ahead(!!outcome, ahead = args_list$data_frequency, skip = TRUE) %>%
     recipes::update_role(!!outcome, new_role = "predictor") %>%
@@ -77,7 +77,6 @@ cdc_baseline_forecaster <- function(
   forecast_date <- args_list$forecast_date %||% max(epi_data$time_value)
   # target_date <- args_list$target_date %||% (forecast_date + args_list$ahead)
 
-
   latest <- get_test_data(epi_recipe(epi_data), epi_data)
 
   f <- frosting() %>%
@@ -93,7 +92,9 @@ cdc_baseline_forecaster <- function(
     layer_add_forecast_date(forecast_date = forecast_date) %>%
     layer_unnest(.pred_distn_all)
   # layer_add_target_date(target_date = target_date)
-  if (args_list$nonneg) f <- layer_threshold(f, ".pred")
+  if (args_list$nonneg) {
+    f <- layer_threshold(f, ".pred")
+  }
 
   eng <- linear_reg(engine = "flatline")
 
@@ -118,7 +119,6 @@ cdc_baseline_forecaster <- function(
 }
 
 
-
 #' CDC baseline forecaster argument constructor
 #'
 #' Constructs a list of arguments for [cdc_baseline_forecaster()].
@@ -157,16 +157,17 @@ cdc_baseline_forecaster <- function(
 #' cdc_baseline_args_list(symmetrize = FALSE)
 #' cdc_baseline_args_list(quantile_levels = c(.1, .3, .7, .9), n_training = 120)
 cdc_baseline_args_list <- function(
-    data_frequency = "1 week",
-    aheads = 1:5,
-    n_training = Inf,
-    forecast_date = NULL,
-    quantile_levels = c(.01, .025, 1:19 / 20, .975, .99),
-    nsims = 1e5L,
-    symmetrize = TRUE,
-    nonneg = TRUE,
-    quantile_by_key = "geo_value",
-    ...) {
+  data_frequency = "1 week",
+  aheads = 1:5,
+  n_training = Inf,
+  forecast_date = NULL,
+  quantile_levels = c(.01, .025, 1:19 / 20, .975, .99),
+  nsims = 1e5L,
+  symmetrize = TRUE,
+  nonneg = TRUE,
+  quantile_by_key = "geo_value",
+  ...
+) {
   rlang::check_dots_empty()
   arg_is_scalar(n_training, nsims, data_frequency)
   data_frequency <- parse_period(data_frequency)
@@ -178,7 +179,9 @@ cdc_baseline_args_list <- function(
   arg_is_lgl(symmetrize, nonneg)
   arg_is_probabilities(quantile_levels, allow_null = TRUE)
   arg_is_pos(n_training)
-  if (is.finite(n_training)) arg_is_pos_int(n_training)
+  if (is.finite(n_training)) {
+    arg_is_pos_int(n_training)
+  }
 
   structure(
     enlist(
@@ -206,10 +209,13 @@ parse_period <- function(x) {
   arg_is_scalar(x)
   if (is.character(x)) {
     x <- unlist(strsplit(x, " "))
-    if (length(x) == 1L) x <- as.numeric(x)
+    if (length(x) == 1L) {
+      x <- as.numeric(x)
+    }
     if (length(x) == 2L) {
       mult <- substr(x[2], 1, 3)
-      mult <- switch(mult,
+      mult <- switch(
+        mult,
         day = 1L,
         wee = 7L,
         cli_abort("incompatible timespan in `aheads`.")
diff --git a/R/flatline_forecaster.R b/R/flatline_forecaster.R
index 2f6d4d2e..8cf9bc84 100644
--- a/R/flatline_forecaster.R
+++ b/R/flatline_forecaster.R
@@ -10,8 +10,8 @@
 #'   combination of key values (`geo_value` + any additional keys) in the
 #'   `epi_data` argument.
 #'
-#' This forecaster is very similar to that used by the
-#'   [COVID19ForecastHub](https://covid19forecasthub.org)
+#' This forecaster is very similar to that used as the
+#'   as [FluSight Baseline](https://github.com/cdcepi/Flusight-baseline/).
 #'
 #' @details
 #'  Here is (roughly) the code for the `flatline_forecaster()` applied to the
@@ -52,9 +52,10 @@
 #'
 #' out <- flatline_forecaster(jhu, "death_rate")
 flatline_forecaster <- function(
-    epi_data,
-    outcome,
-    args_list = flatline_args_list()) {
+  epi_data,
+  outcome,
+  args_list = flatline_args_list()
+) {
   validate_forecaster_inputs(epi_data, outcome, "time_value")
   if (!inherits(args_list, c("flat_fcast", "alist"))) {
     cli_abort("`args_list` was not created using `flatline_args_list()`.")
@@ -77,7 +78,6 @@ flatline_forecaster <- function(
   ek <- kill_time_value(keys)
   outcome <- rlang::sym(outcome)
 
-
   r <- epi_recipe(epi_data) %>%
     step_epi_ahead(!!outcome, ahead = args_list$ahead, skip = TRUE) %>%
     recipes::update_role(!!outcome, new_role = "predictor") %>%
@@ -96,7 +96,9 @@ flatline_forecaster <- function(
     ) %>%
     layer_add_forecast_date(forecast_date = forecast_date) %>%
     layer_add_target_date(target_date = target_date)
-  if (args_list$nonneg) f <- layer_threshold(f, dplyr::starts_with(".pred"))
+  if (args_list$nonneg) {
+    f <- layer_threshold(f, dplyr::starts_with(".pred"))
+  }
 
   eng <- linear_reg(engine = "flatline")
 
@@ -120,7 +122,6 @@ flatline_forecaster <- function(
 }
 
 
-
 #' Flatline forecaster argument constructor
 #'
 #' Constructs a list of arguments for [flatline_forecaster()].
@@ -141,15 +142,16 @@ flatline_forecaster <- function(
 #' flatline_args_list(symmetrize = FALSE)
 #' flatline_args_list(quantile_levels = c(.1, .3, .7, .9), n_training = 120)
 flatline_args_list <- function(
-    ahead = 7L,
-    n_training = Inf,
-    forecast_date = NULL,
-    target_date = NULL,
-    quantile_levels = c(0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95),
-    symmetrize = TRUE,
-    nonneg = TRUE,
-    quantile_by_key = character(0L),
-    ...) {
+  ahead = 7L,
+  n_training = Inf,
+  forecast_date = NULL,
+  target_date = NULL,
+  quantile_levels = c(0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95),
+  symmetrize = TRUE,
+  nonneg = TRUE,
+  quantile_by_key = character(0L),
+  ...
+) {
   rlang::check_dots_empty()
   arg_is_scalar(ahead, n_training)
   arg_is_chr(quantile_by_key, allow_empty = TRUE)
@@ -159,7 +161,9 @@ flatline_args_list <- function(
   arg_is_lgl(symmetrize, nonneg)
   arg_is_probabilities(quantile_levels, allow_null = TRUE)
   arg_is_pos(n_training)
-  if (is.finite(n_training)) arg_is_pos_int(n_training)
+  if (is.finite(n_training)) {
+    arg_is_pos_int(n_training)
+  }
 
   if (!is.null(forecast_date) && !is.null(target_date)) {
     if (forecast_date + ahead != target_date) {
diff --git a/R/pivot_quantiles.R b/R/pivot_quantiles.R
index 7ae0b461..fa7cb604 100644
--- a/R/pivot_quantiles.R
+++ b/R/pivot_quantiles.R
@@ -82,7 +82,12 @@ pivot_quantiles_longer <- function(.data, ...) {
   names(long_tib)[1:2] <- c(glue::glue("{col}_value"), glue::glue("{col}_quantile_level"))
   out <- left_join(.data, long_tib, by = ".row") %>% select(!.row)
   if (inherits(.data, "epi_df")) {
-    attr(out, "metadata")$other_keys <- c(attr(.data, "metadata")$other_keys, glue::glue("{col}_quantile_level"))
+    # register the new quantile_level column as a key
+    out <- as_epi_df(
+      out,
+      as_of = attr(.data, "metadata")$as_of,
+      other_keys = c(attr(.data, "metadata")$other_keys, glue::glue("{col}_quantile_level"))
+    )
   }
   out
 }
diff --git a/README.Rmd b/README.Rmd
index 28ae1409..31956753 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -11,7 +11,7 @@ knitr::opts_chunk$set(
   digits = 3,
   comment = "#>",
   collapse = TRUE,
-  cache = TRUE,
+  cache = FALSE,
   dev.args = list(bg = "transparent"),
   dpi = 300,
   cache.lazy = FALSE,
@@ -51,15 +51,20 @@ tertiary <- "#177245"
 fourth_colour <- "#A393BF"
 fifth_colour <- "#2e8edd"
 colvec <- c(
-  base = base, primary = primary, secondary = secondary,
-  tertiary = tertiary, fourth_colour = fourth_colour,
+  base = base,
+  primary = primary,
+  secondary = secondary,
+  tertiary = tertiary,
+  fourth_colour = fourth_colour,
   fifth_colour = fifth_colour
 )
 library(epiprocess)
 suppressMessages(library(tidyverse))
 theme_update(legend.position = "bottom", legend.title = element_blank())
 delphi_pal <- function(n) {
-  if (n > 6L) warning("Not enough colors in this palette!")
+  if (n > 6L) {
+    warning("Not enough colors in this palette!")
+  }
   unname(colvec)[1:n]
 }
 scale_fill_delphi <- function(..., aesthetics = "fill") {
@@ -87,16 +92,17 @@ Both of the packages are meant to work well with the panel data provided by
 Pre-compiled example datasets are also available in
 [`{epidatasets}`](https://cmu-delphi.github.io/epidatasets/).
 
+<!--
 If you are looking for detail beyond the package documentation, see our
 [forecasting book](https://cmu-delphi.github.io/delphi-tooling-book/).
-
+-->
 
 ## Installation
 
 Unless you’re planning on contributing to package development, we suggest using the stable version.
 To install, run:
 
-```r
+```{r echo=TRUE, eval=FALSE}
 # Stable version
 pak::pkg_install("cmu-delphi/epipredict@main")
 
@@ -105,8 +111,8 @@ pak::pkg_install("cmu-delphi/epipredict@dev")
 ```
 
 The documentation for the stable version is at
-<https://cmu-delphi.github.io/epipredict>, while the development version is at
-<https://cmu-delphi.github.io/epipredict/dev>.
+<https://cmu-delphi.github.io/epipredict/>, while the development version is at
+<https://cmu-delphi.github.io/epipredict/dev/>.
 
 
 ## Motivating example
@@ -151,7 +157,7 @@ First we pull both `jhu-csse` cases and deaths data from the
 [Delphi API](https://cmu-delphi.github.io/delphi-epidata/api/covidcast.html) using the
 [`{epidatr}`](https://cmu-delphi.github.io/epidatr/) package:
 
-```{r case_death, warning = FALSE, eval = TRUE}
+```{r case_death, warning = FALSE, eval = FALSE}
 cases <- pub_covidcast(
   source = "jhu-csse",
   signals = "confirmed_7dav_incidence_prop",
@@ -177,6 +183,15 @@ cases_deaths <-
   as_epi_df(as_of = as.Date("2022-01-01"))
 ```
 
+```{r case_death_epidatasets, echo = FALSE}
+cases_deaths <- epidatasets::covid_case_death_rates_extended |>
+  filter(
+    geo_value %in% used_locations,
+    between(time_value, ymd("2020-06-01"), ymd("2021-12-31"))
+  ) |>
+  as_epi_df(as_of = as.Date("2022-01-01"))
+```
+
 Since visualizing the results on every geography is somewhat overwhelming,
 we’ll only train on a subset of locations.
 
@@ -208,11 +223,13 @@ cases_deaths <-
   group_by(geo_value) |>
   mutate(
     outlr_death_rate = detect_outlr_rm(
-      time_value, death_rate,
+      time_value,
+      death_rate,
       detect_negatives = TRUE
     ),
     outlr_case_rate = detect_outlr_rm(
-      time_value, case_rate,
+      time_value,
+      case_rate,
       detect_negatives = TRUE
     )
   ) |>
@@ -257,7 +274,8 @@ processed_data_plot <-
   geom_text(
     data = forecast_date_label,
     aes(x = dates, label = "forecast\ndate", y = heights),
-    size = 3, hjust = "right"
+    size = 3,
+    hjust = "right"
   ) +
   scale_x_date(date_breaks = "3 months", date_labels = "%Y %b") +
   theme(axis.text.x = element_text(angle = 90, hjust = 1))
@@ -315,7 +333,8 @@ forecast_plot <-
   geom_text(
     data = forecast_date_label %>% filter(.response_name == "death_rate"),
     aes(x = dates, label = "forecast\ndate", y = heights),
-    size = 3, hjust = "right"
+    size = 3,
+    hjust = "right"
   ) +
   scale_x_date(date_breaks = "3 months", date_labels = "%y %b") +
   theme(axis.text.x = element_text(angle = 90, hjust = 1))
@@ -331,7 +350,13 @@ And as a tibble of quantile level-value pairs:
 four_week_ahead$predictions |>
   select(-.pred) |>
   pivot_quantiles_longer(.pred_distn) |>
-  select(geo_value, forecast_date, target_date, quantile = .pred_distn_quantile_level, value = .pred_distn_value)
+  select(
+    geo_value,
+    forecast_date,
+    target_date,
+    quantile = .pred_distn_quantile_level,
+    value = .pred_distn_value
+  )
 ```
 
 The orange dot gives the point prediction, while the blue intervals give the
diff --git a/README.md b/README.md
index 97ee8121..0ebe963d 100644
--- a/README.md
+++ b/README.md
@@ -23,8 +23,10 @@ provided by [`{epidatr}`](https://cmu-delphi.github.io/epidatr/).
 Pre-compiled example datasets are also available in
 [`{epidatasets}`](https://cmu-delphi.github.io/epidatasets/).
 
+<!--
 If you are looking for detail beyond the package documentation, see our
 [forecasting book](https://cmu-delphi.github.io/delphi-tooling-book/).
+-->
 
 ## Installation
 
@@ -40,13 +42,15 @@ pak::pkg_install("cmu-delphi/epipredict@dev")
 ```
 
 The documentation for the stable version is at
-<https://cmu-delphi.github.io/epipredict>, while the development version
-is at <https://cmu-delphi.github.io/epipredict/dev>.
+<https://cmu-delphi.github.io/epipredict/>, while the development
+version is at <https://cmu-delphi.github.io/epipredict/dev/>.
 
 ## Motivating example
 
 <details>
+
 <summary>
+
 Required packages
 </summary>
 
@@ -81,7 +85,9 @@ data](https://cmu-delphi.github.io/delphi-epidata/api/covidcast-signals/jhu-csse
 Below the fold, we pull the dataset from the epidata API and clean it.
 
 <details>
+
 <summary>
+
 Creating the dataset using `{epidatr}` and `{epiprocess}`
 </summary>
 
@@ -141,7 +147,7 @@ cases_deaths |>
   theme(axis.text.x = element_text(angle = 90, hjust = 1))
 ```
 
-<img src="man/figures/README-date-1.png" width="90%" style="display: block; margin: auto;" />
+<img src="man/figures/README-date-1.png" alt="" width="90%" style="display: block; margin: auto;" />
 
 As with the typical dataset, we will need to do some cleaning to make it
 actually usable; we’ll use some utilities from
@@ -154,11 +160,13 @@ cases_deaths <-
   group_by(geo_value) |>
   mutate(
     outlr_death_rate = detect_outlr_rm(
-      time_value, death_rate,
+      time_value,
+      death_rate,
       detect_negatives = TRUE
     ),
     outlr_case_rate = detect_outlr_rm(
-      time_value, case_rate,
+      time_value,
+      case_rate,
       detect_negatives = TRUE
     )
   ) |>
@@ -178,7 +186,9 @@ After downloading and cleaning deaths per capita, as well as cases per
 forecast date with a vertical line:
 
 <details>
+
 <summary>
+
 Plot
 </summary>
 
@@ -207,7 +217,8 @@ processed_data_plot <-
   geom_text(
     data = forecast_date_label,
     aes(x = dates, label = "forecast\ndate", y = heights),
-    size = 3, hjust = "right"
+    size = 3,
+    hjust = "right"
   ) +
   scale_x_date(date_breaks = "3 months", date_labels = "%Y %b") +
   theme(axis.text.x = element_text(angle = 90, hjust = 1))
@@ -215,7 +226,7 @@ processed_data_plot <-
 
 </details>
 
-<img src="man/figures/README-show-processed-data-1.png" width="90%" style="display: block; margin: auto;" />
+<img src="man/figures/README-show-processed-data-1.png" alt="" width="90%" style="display: block; margin: auto;" />
 
 To make a forecast, we will use a simple “canned” auto-regressive
 forecaster to predict the death rate four weeks into the future using
@@ -235,12 +246,12 @@ four_week_ahead <- arx_forecaster(
 four_week_ahead
 #> ══ A basic forecaster of type ARX Forecaster ════════════════════════════════
 #> 
-#> This forecaster was fit on 2025-05-22 11:56:44.
+#> This forecaster was fit on 2026-06-30 18:35:46.
 #> 
 #> Training data was an <epi_df> with:
 #> • Geography: state,
 #> • Time type: day,
-#> • Using data up-to-date as of: 2022-01-01.
+#> • Using data up-to-date as of: 2023-03-10.
 #> • With the last data available on 2021-08-01
 #> 
 #> ── Predictions ──────────────────────────────────────────────────────────────
@@ -263,7 +274,9 @@ Plotting the prediction intervals on the true values for our location
 subset[^2]:
 
 <details>
+
 <summary>
+
 Plot
 </summary>
 
@@ -280,7 +293,8 @@ forecast_plot <-
   geom_text(
     data = forecast_date_label %>% filter(.response_name == "death_rate"),
     aes(x = dates, label = "forecast\ndate", y = heights),
-    size = 3, hjust = "right"
+    size = 3,
+    hjust = "right"
   ) +
   scale_x_date(date_breaks = "3 months", date_labels = "%y %b") +
   theme(axis.text.x = element_text(angle = 90, hjust = 1))
@@ -288,7 +302,7 @@ forecast_plot <-
 
 </details>
 
-<img src="man/figures/README-show-single-forecast-1.png" width="90%" style="display: block; margin: auto;" />
+<img src="man/figures/README-show-single-forecast-1.png" alt="" width="90%" style="display: block; margin: auto;" />
 
 And as a tibble of quantile level-value pairs:
 
@@ -296,7 +310,13 @@ And as a tibble of quantile level-value pairs:
 four_week_ahead$predictions |>
   select(-.pred) |>
   pivot_quantiles_longer(.pred_distn) |>
-  select(geo_value, forecast_date, target_date, quantile = .pred_distn_quantile_level, value = .pred_distn_value)
+  select(
+    geo_value,
+    forecast_date,
+    target_date,
+    quantile = .pred_distn_quantile_level,
+    value = .pred_distn_value
+  )
 #> # A tibble: 20 × 5
 #>   geo_value forecast_date target_date quantile  value
 #>   <chr>     <date>        <date>         <dbl>  <dbl>
diff --git a/man/arx_class_args_list.Rd b/man/arx_class_args_list.Rd
index 7359c876..919894b3 100644
--- a/man/arx_class_args_list.Rd
+++ b/man/arx_class_args_list.Rd
@@ -67,8 +67,8 @@ latency is large. If this is \code{FALSE}, that warning is turned off.}
 be created using growth rates (as the predictors are) or lagged
 differences. The second case is closer to the requirements for the
 \href{https://github.com/cdcepi/Flusight-forecast-data/blob/745511c436923e1dc201dea0f4181f21a8217b52/data-experimental/README.md}{2022-23 CDC Flusight Hospitalization Experimental Target}.
-See the \href{https://cmu-delphi.github.io/delphi-tooling-book/arx-classifier.html}{Classification chapter from the forecasting book} Vignette for details of how to create a reasonable
-baseline for this case. Selecting \code{"growth_rate"} (the default) uses
+See the \code{vignette("epipredict")} for more details.
+Selecting \code{"growth_rate"} (the default) uses
 \code{\link[epiprocess:growth_rate]{epiprocess::growth_rate()}} to create the outcome using some of the
 additional arguments below. Choosing \code{"lag_difference"} instead simply
 uses the change from the value at the selected \code{horizon}.}
diff --git a/man/cdc_baseline_forecaster.Rd b/man/cdc_baseline_forecaster.Rd
index c201bd99..08e741e2 100644
--- a/man/cdc_baseline_forecaster.Rd
+++ b/man/cdc_baseline_forecaster.Rd
@@ -33,8 +33,8 @@ forward over all available training data.
 By default, the predictive intervals are computed separately for each
 combination of \code{geo_value} in the \code{epi_data} argument.
 
-This forecaster is meant to produce exactly the CDC Baseline used for
-\href{https://covid19forecasthub.org}{COVID19ForecastHub}
+This forecaster is meant to produce exactly the CDC Baseline used for, e.g.,
+\url{https://github.com/cdcepi/Flusight-baseline/}
 }
 \examples{
 library(dplyr)
diff --git a/man/figures/README-date-1.png b/man/figures/README-date-1.png
index b66ec04c..15eb4f10 100644
Binary files a/man/figures/README-date-1.png and b/man/figures/README-date-1.png differ
diff --git a/man/figures/README-show-processed-data-1.png b/man/figures/README-show-processed-data-1.png
index e312c9fd..0f7973d7 100644
Binary files a/man/figures/README-show-processed-data-1.png and b/man/figures/README-show-processed-data-1.png differ
diff --git a/man/figures/README-show-single-forecast-1.png b/man/figures/README-show-single-forecast-1.png
index 8ec4f90a..3d880600 100644
Binary files a/man/figures/README-show-single-forecast-1.png and b/man/figures/README-show-single-forecast-1.png differ
diff --git a/man/flatline_forecaster.Rd b/man/flatline_forecaster.Rd
index cc789bac..bbe0576b 100644
--- a/man/flatline_forecaster.Rd
+++ b/man/flatline_forecaster.Rd
@@ -29,8 +29,8 @@ By default, the predictive intervals are computed separately for each
 combination of key values (\code{geo_value} + any additional keys) in the
 \code{epi_data} argument.
 
-This forecaster is very similar to that used by the
-\href{https://covid19forecasthub.org}{COVID19ForecastHub}
+This forecaster is very similar to that used as the
+as \href{https://github.com/cdcepi/Flusight-baseline/}{FluSight Baseline}.
 }
 \details{
 Here is (roughly) the code for the \code{flatline_forecaster()} applied to the
diff --git a/man/layer_predict.Rd b/man/layer_predict.Rd
index efaff62e..d678c9fa 100644
--- a/man/layer_predict.Rd
+++ b/man/layer_predict.Rd
@@ -15,18 +15,10 @@ layer_predict(
 \arguments{
 \item{frosting}{a frosting object}
 
-\item{type}{A single character value or \code{NULL}. Possible values are:
-\itemize{
-\item regression: "\code{numeric}"
-\item classification: "\code{class}", "\code{prob}"
-\item censored regression: "\code{survival}", "\code{time}", "\code{hazard}", "\code{linear_pred}"
-\item quantile regression: "\code{quantile}"
-\item interval estimates: "\code{conf_int}", "\code{pred_int}"
-\item other:  "\code{raw}"
-}
-
-When \code{NULL}, \code{predict()} will choose an appropriate value based on the
-model's mode.}
+\item{type}{A single character value or \code{NULL}. Possible values
+are \code{"numeric"}, \code{"class"}, \code{"prob"}, \code{"conf_int"}, \code{"pred_int"},
+\code{"quantile"}, \code{"time"}, \code{"hazard"}, \code{"survival"}, or \code{"raw"}. When \code{NULL},
+\code{predict()} will choose an appropriate value based on the model's mode.}
 
 \item{opts}{A list of optional arguments to the underlying
 predict function that will be used when \code{type = "raw"}. The
@@ -38,7 +30,7 @@ value of \code{type}. Arguments to the underlying model's prediction
 function cannot be passed here (use the \code{opts} argument instead).
 Possible arguments are:
 \itemize{
-\item \code{interval}: for \code{type} equal to \code{"survival"}, should
+\item \code{interval}: for \code{type} equal to \code{"survival"} or \code{"quantile"}, should
 interval estimates be added, if available? Options are \code{"none"}
 and \code{"confidence"}.
 \item \code{level}: for \code{type} equal to \code{"conf_int"}, \code{"pred_int"}, or \code{"survival"},
@@ -48,6 +40,8 @@ Default value is \code{0.95}.
 \item \code{std_error}: for \code{type} equal to \code{"conf_int"} or \code{"pred_int"}, add
 the standard error of fit or prediction (on the scale of the
 linear predictors). Default value is \code{FALSE}.
+\item \code{quantile}: for \code{type} equal to \code{quantile}, the quantiles of the
+distribution. Default is \code{(1:9)/10}.
 \item \code{eval_time}: for \code{type} equal to \code{"survival"} or \code{"hazard"}, the
 time points at which the survival probability or hazard is estimated.
 }}
diff --git a/man/predict-epi_workflow.Rd b/man/predict-epi_workflow.Rd
index f9a37784..577b20ab 100644
--- a/man/predict-epi_workflow.Rd
+++ b/man/predict-epi_workflow.Rd
@@ -14,18 +14,10 @@
 \item{new_data}{A data frame containing the new predictors to preprocess
 and predict on}
 
-\item{type}{A single character value or \code{NULL}. Possible values are:
-\itemize{
-\item regression: "\code{numeric}"
-\item classification: "\code{class}", "\code{prob}"
-\item censored regression: "\code{survival}", "\code{time}", "\code{hazard}", "\code{linear_pred}"
-\item quantile regression: "\code{quantile}"
-\item interval estimates: "\code{conf_int}", "\code{pred_int}"
-\item other:  "\code{raw}"
-}
-
-When \code{NULL}, \code{predict()} will choose an appropriate value based on the
-model's mode.}
+\item{type}{A single character value or \code{NULL}. Possible values
+are \code{"numeric"}, \code{"class"}, \code{"prob"}, \code{"conf_int"}, \code{"pred_int"},
+\code{"quantile"}, \code{"time"}, \code{"hazard"}, \code{"survival"}, or \code{"raw"}. When \code{NULL},
+\code{predict()} will choose an appropriate value based on the model's mode.}
 
 \item{opts}{A list of optional arguments to the underlying
 predict function that will be used when \code{type = "raw"}. The
@@ -37,7 +29,7 @@ value of \code{type}. Arguments to the underlying model's prediction
 function cannot be passed here (use the \code{opts} argument instead).
 Possible arguments are:
 \itemize{
-\item \code{interval}: for \code{type} equal to \code{"survival"}, should
+\item \code{interval}: for \code{type} equal to \code{"survival"} or \code{"quantile"}, should
 interval estimates be added, if available? Options are \code{"none"}
 and \code{"confidence"}.
 \item \code{level}: for \code{type} equal to \code{"conf_int"}, \code{"pred_int"}, or \code{"survival"},
@@ -47,6 +39,8 @@ Default value is \code{0.95}.
 \item \code{std_error}: for \code{type} equal to \code{"conf_int"} or \code{"pred_int"}, add
 the standard error of fit or prediction (on the scale of the
 linear predictors). Default value is \code{FALSE}.
+\item \code{quantile}: for \code{type} equal to \code{quantile}, the quantiles of the
+distribution. Default is \code{(1:9)/10}.
 \item \code{eval_time}: for \code{type} equal to \code{"survival"} or \code{"hazard"}, the
 time points at which the survival probability or hazard is estimated.
 }}
diff --git a/man/step_adjust_latency.Rd b/man/step_adjust_latency.Rd
index a77a86d2..d0d6c3ef 100644
--- a/man/step_adjust_latency.Rd
+++ b/man/step_adjust_latency.Rd
@@ -141,8 +141,6 @@ toy_recipe \%>\%
 #> * geo_type  = state
 #> * time_type = day
 #> * as_of     = 2015-01-14
-#> Latency (time between last available observation and epi_df's as_of, by time series):
-#> * latency across all time series = 0 days
 #> 
 #> # A tibble: 8 x 4
 #>   geo_value time_value     a     b
@@ -179,8 +177,6 @@ toy_recipe \%>\%
 #> * geo_type  = state
 #> * time_type = day
 #> * as_of     = 2015-01-14
-#> Latency (time between last available observation and epi_df's as_of, by time series):
-#> * latency  = -2–3 days
 #> 
 #> # A tibble: 21 x 7
 #>    geo_value time_value     a     b lag_3_a lag_4_b ahead_1_a
@@ -229,8 +225,6 @@ toy_recipe \%>\%
 #> * geo_type  = state
 #> * time_type = day
 #> * as_of     = 2015-01-14
-#> Latency (time between last available observation and epi_df's as_of, by time series):
-#> * latency  = 1–5 days
 #> 
 #> # A tibble: 10 x 6
 #>    geo_value time_value     a     b lag_0_a ahead_3_a
diff --git a/tests/testthat/test-autoplot.R b/tests/testthat/test-autoplot.R
new file mode 100644
index 00000000..3de90989
--- /dev/null
+++ b/tests/testthat/test-autoplot.R
@@ -0,0 +1,51 @@
+
+train_data <- epidatasets::cases_deaths_subset
+
+test_that("autoplot renders a single-target-date forecast (geom_linerange band)", {
+  fc <- arx_forecaster(
+    train_data,
+    "death_rate_7d_av",
+    c("death_rate_7d_av", "case_rate_7d_av"),
+    args_list = arx_args_list(ahead = 7L, quantile_levels = c(.1, .5, .9))
+  )
+  p <- autoplot(fc, observed_response = train_data)
+  expect_s3_class(p, "ggplot")
+  expect_silent(ggplot2::ggplot_build(p))
+})
+
+test_that("autoplot renders a multi-target-date forecast (geom_ribbon band)", {
+
+  forecast_date <- as.Date("2021-08-01")
+
+  # arx_forecaster only takes a scalar `ahead`, so a real multi-horizon
+  # forecast means forecasting each ahead separately and combining them.
+  all_canned_results <- lapply(
+    seq(0, 28),
+    \(days_ahead) {
+      arx_forecaster(
+        train_data |>
+          filter(time_value <= forecast_date),
+        outcome = "death_rate_7d_av",
+        predictors = c("case_rate_7d_av", "death_rate_7d_av"),
+        trainer = quantile_reg(),
+        args_list = arx_args_list(
+          lags = list(c(0, 1, 2, 3, 7, 14), c(0, 7, 14)),
+          ahead = days_ahead
+        )
+      )
+    }
+  )
+  # pull out the workflow and the predictions to be able to use autoplot
+  workflow <- all_canned_results[[1]]$epi_workflow
+  results <- all_canned_results |>
+    purrr::map(~ `$`(., "predictions")) |>
+    purrr::list_rbind()
+  p <- autoplot(
+    object = workflow,
+    predictions = results,
+    observed_response = train_data |>
+      filter(time_value > "2021-07-01")
+  )
+  expect_s3_class(p, "ggplot")
+  expect_silent(ggplot2::ggplot_build(p))
+})
diff --git a/vignettes/backtesting.Rmd b/vignettes/backtesting.Rmd
index f355f551..4039f80a 100644
--- a/vignettes/backtesting.Rmd
+++ b/vignettes/backtesting.Rmd
@@ -24,7 +24,7 @@ new data is constantly being collected and added to the dataset, and old data po
 Training and making
 predictions only on finalized data can lead to overly optimistic estimates of accuracy
 (see, for example, [McDonald et al.
-2021](https://www.pnas.org/content/118/51/e2111453118/) and the references
+2021](https://doi.org/10.1073/pnas.2111453118) and the references
 therein).
 
 In the `{epiprocess}` package, we provide the function `epix_slide()` to help conviently perform version-faithful forecasting by only using the data as
@@ -95,7 +95,8 @@ geo_choose <- "ca"
 forecast_dates <- seq(
   from = as.Date("2020-08-01"),
   to = as.Date("2021-11-01"),
-  by = "1 month")
+  by = "1 month"
+)
 percent_cli_data <- bind_rows(
   # Snapshotted data for the version-faithful forecasts
   map(
@@ -110,18 +111,22 @@ percent_cli_data <- bind_rows(
   doctor_visits |>
     epix_as_of(doctor_visits$versions_end) |>
     mutate(version_faithful = "Version un-faithful")
-) |> as_tibble()
+) |>
+  as_tibble()
 p0 <- autoplot(
-  archive_cases_dv_subset, percent_cli, 
-  .versions = forecast_dates, 
+  archive_cases_dv_subset,
+  percent_cli,
+  .versions = forecast_dates,
   .mark_versions = TRUE,
   .facet_filter = (geo_value == "ca")
 ) +
   scale_x_date(minor_breaks = "month", date_labels = "%b %Y") +
-  labs(x = "", y = "% of doctor's visits with\n Covid-like illness") + 
+  labs(x = "", y = "% of doctor's visits with\n Covid-like illness") +
   scale_color_viridis_c(
     option = "viridis",
-    guide = guide_legend(reverse=TRUE), direction = -1) +
+    guide = guide_legend(reverse = TRUE),
+    direction = -1
+  ) +
   scale_y_continuous(limits = c(0, NA), expand = expansion(c(0, 0.05))) +
   theme(legend.position = "none")
 ```
@@ -246,21 +251,29 @@ We will also create the helper function `forecast_wrapper()` to let us easily ma
 
 ```{r arx-kweek-preliminaries, warning = FALSE}
 forecast_wrapper <- function(
-    epi_data, aheads, outcome, predictors, process_data = identity
-    ) {
+  epi_data,
+  aheads,
+  outcome,
+  predictors,
+  process_data = identity
+) {
   map(
     aheads,
     \(ahead) {
       arx_forecaster(
-        process_data(epi_data), outcome, predictors,
+        process_data(epi_data),
+        outcome,
+        predictors,
         args_list = arx_args_list(
           ahead = ahead,
           lags = c(0:7, 14, 21),
           adjust_latency = "extend_ahead"
         )
-      )$predictions |> pivot_quantiles_wider(.pred_distn)
+      )$predictions |>
+        pivot_quantiles_wider(.pred_distn)
     }
-  ) |> bind_rows()
+  ) |>
+    bind_rows()
 }
 ```
 
@@ -332,11 +345,14 @@ plotting_data <- bind_rows(
     mutate(version_faithful = "Version faithful")
 )
 
-p1 <- ggplot(data = forecasts_filtered,
-             aes(x = target_date, group = time_value)) +
+p1 <- ggplot(
+  data = forecasts_filtered,
+  aes(x = target_date, group = time_value)
+) +
   geom_ribbon(
     aes(ymin = `0.05`, ymax = `0.95`, fill = (time_value)),
-    alpha = 0.4) +
+    alpha = 0.4
+  ) +
   geom_line(aes(y = .pred, color = (time_value)), linetype = 2L) +
   geom_point(aes(y = .pred, color = (time_value)), size = 0.75) +
   # the forecast date
@@ -351,13 +367,16 @@ p1 <- ggplot(data = forecasts_filtered,
   geom_line(
     data = plotting_data |> filter(geo_value == geo_choose),
     aes(x = time_value, y = percent_cli, color = (version), group = version),
-    inherit.aes = FALSE, na.rm = TRUE
+    inherit.aes = FALSE,
+    na.rm = TRUE
   ) +
   facet_grid(version_faithful ~ geo_value, scales = "free") +
   scale_x_date(breaks = "2 months", date_labels = "%b %Y") +
   scale_y_continuous(expand = expansion(c(0, 0.05))) +
-  labs(x = "Date",
-       y = "smoothed, day of week adjusted covid-like doctors visits") +
+  labs(
+    x = "Date",
+    y = "smoothed, day of week adjusted covid-like doctors visits"
+  ) +
   scale_color_viridis_c(option = "viridis", direction = -1) +
   scale_fill_viridis_c(option = "viridis", direction = -1) +
   theme(legend.position = "none")
@@ -372,12 +391,17 @@ forecasts_filtered <- forecasts |>
 forecasts_filtered %>% names
 p2 <-
   ggplot(data = forecasts_filtered, aes(x = target_date, group = time_value)) +
-  geom_ribbon(aes(ymin = `0.05`, ymax = `0.95`, fill = (time_value)), alpha = 0.4) +
+  geom_ribbon(
+    aes(ymin = `0.05`, ymax = `0.95`, fill = (time_value)),
+    alpha = 0.4
+  ) +
   geom_line(aes(y = .pred, color = (time_value)), linetype = 2L) +
   geom_point(aes(y = .pred, color = (time_value)), size = 0.75) +
   # the forecast date
   geom_vline(
-    data = percent_cli_data |> filter(geo_value == geo_choose) |> select(-version_faithful),
+    data = percent_cli_data |>
+      filter(geo_value == geo_choose) |>
+      select(-version_faithful),
     aes(color = version, xintercept = version, group = version),
     lty = 2
   ) +
@@ -385,12 +409,16 @@ p2 <-
   geom_line(
     data = plotting_data |> filter(geo_value == geo_choose),
     aes(x = time_value, y = percent_cli, color = (version), group = version),
-    inherit.aes = FALSE, na.rm = TRUE
+    inherit.aes = FALSE,
+    na.rm = TRUE
   ) +
   facet_grid(version_faithful ~ geo_value, scales = "free") +
   scale_x_date(breaks = "2 months", date_labels = "%b %Y") +
   scale_y_continuous(expand = expansion(c(0, 0.05))) +
-  labs(x = "Date", y = "smoothed, day of week adjusted covid-like doctors visits") +
+  labs(
+    x = "Date",
+    y = "smoothed, day of week adjusted covid-like doctors visits"
+  ) +
   scale_color_viridis_c(option = "viridis", direction = -1) +
   scale_fill_viridis_c(option = "viridis", direction = -1) +
   theme(legend.position = "none")
diff --git a/vignettes/epipredict.Rmd b/vignettes/epipredict.Rmd
index 3786bc35..ef9a36c6 100644
--- a/vignettes/epipredict.Rmd
+++ b/vignettes/epipredict.Rmd
@@ -272,14 +272,15 @@ results <- all_flatlines |>
 autoplot(
   object = workflow,
   predictions = results,
-  observed_response = covid_case_death_rates |> filter(geo_value %in% used_locations, time_value > "2021-07-01")
+  observed_response = covid_case_death_rates |>
+    filter(geo_value %in% used_locations, time_value > "2021-07-01")
 )
 ```
 
 ### `cdc_baseline_forecaster()`
 
 This is a different method of generating a flatline forecast, used as a baseline
-for [the CDC COVID-19 Forecasting Hub](https://covid19forecasthub.org).
+for various CDC forecasting challenges, e.g., <https://github.com/cdcepi/Flusight-baseline/>.
 
 ```{r make-cdc-forecast, warning=FALSE}
 all_cdc_flatline <-
@@ -298,7 +299,8 @@ results <- all_cdc_flatline$predictions
 autoplot(
   object = workflow,
   predictions = results,
-  observed_response = covid_case_death_rates |> filter(geo_value %in% used_locations, time_value > "2021-07-01")
+  observed_response = covid_case_death_rates |>
+    filter(geo_value %in% used_locations, time_value > "2021-07-01")
 )
 ```
 
@@ -327,14 +329,15 @@ We'll predict the 2023/24 season using all previous data, including 2020-2022, t
 ```{r make-climatological-forecast, warning=FALSE}
 fluview_hhs <- pub_fluview(
   regions = paste0("hhs", 1:10),
-  epiweeks = epirange(100001,222201)
+  epiweeks = epirange(100001, 222201)
 )
 fluview <- fluview_hhs %>%
   select(
     geo_value = region,
     time_value = epiweek,
     issue,
-    ili) %>%
+    ili
+  ) %>%
   as_epi_archive() %>%
   epix_as_of_current()
 
@@ -451,7 +454,8 @@ autoplot(
   grad_forecast$epi_workflow,
   grad_forecast$predictions,
   observed_response = grad_employ,
-) + geom_vline(aes(xintercept = 2016))
+) +
+  geom_vline(aes(xintercept = 2016))
 ```
 
 The 8 graphs represent all combinations of the `geo_values` (`"Quebec"` and `"British Columbia"`), `edu_quals` (`"Undergraduate degree"` and `"Professional degree"`), and age brackets (`"15 to 34 years"` and `"35 to 64 years"`).