diff --git a/NAMESPACE b/NAMESPACE
index a24e5844b..b3e0bc8ae 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -167,7 +167,7 @@ import(distributional)
 import(epiprocess)
 import(parsnip)
 import(recipes)
-import(vctrs)
+importFrom(cli,cli_abort)
 importFrom(epiprocess,growth_rate)
 importFrom(generics,augment)
 importFrom(generics,fit)
@@ -202,3 +202,12 @@ importFrom(stats,residuals)
 importFrom(tibble,as_tibble)
 importFrom(tibble,is_tibble)
 importFrom(tibble,tibble)
+importFrom(vctrs,as_list_of)
+importFrom(vctrs,field)
+importFrom(vctrs,new_rcrd)
+importFrom(vctrs,new_vctr)
+importFrom(vctrs,vec_cast)
+importFrom(vctrs,vec_data)
+importFrom(vctrs,vec_ptype_abbr)
+importFrom(vctrs,vec_ptype_full)
+importFrom(vctrs,vec_recycle_common)
diff --git a/R/arx_classifier.R b/R/arx_classifier.R
index 9370da423..5b4dc0477 100644
--- a/R/arx_classifier.R
+++ b/R/arx_classifier.R
@@ -250,7 +250,8 @@ arx_class_args_list <- function(
     method = c("rel_change", "linear_reg", "smooth_spline", "trend_filter"),
     log_scale = FALSE,
     additional_gr_args = list(),
-    nafill_buffer = Inf) {
+    nafill_buffer = Inf,
+    ...) {
   .lags <- lags
   if (is.list(lags)) lags <- unlist(lags)
   method <- match.arg(method)
@@ -305,3 +306,5 @@ print.arx_class <- function(x, ...) {
   name <- "ARX Classifier"
   NextMethod(name = name, ...)
 }
+
+# this is a trivial change to induce a check
diff --git a/R/arx_forecaster.R b/R/arx_forecaster.R
index 2e242d770..3dc54286e 100644
--- a/R/arx_forecaster.R
+++ b/R/arx_forecaster.R
@@ -33,7 +33,7 @@
 #' out <- arx_forecaster(jhu, "death_rate",
 #'   c("case_rate", "death_rate"),
 #'   trainer = quantile_reg(),
-#'   args_list = arx_args_list(levels = 1:9 / 10)
+#'   args_list = arx_args_list(quantile_levels = 1:9 / 10)
 #' )
 arx_forecaster <- function(epi_data,
                            outcome,
@@ -99,7 +99,7 @@ arx_forecaster <- function(epi_data,
 #' arx_fcast_epi_workflow(jhu, "death_rate",
 #'   c("case_rate", "death_rate"),
 #'   trainer = quantile_reg(),
-#'   args_list = arx_args_list(levels = 1:9 / 10)
+#'   args_list = arx_args_list(quantile_levels = 1:9 / 10)
 #' )
 arx_fcast_epi_workflow <- function(
     epi_data,
@@ -134,18 +134,20 @@ arx_fcast_epi_workflow <- function(
   # --- postprocessor
   f <- frosting() %>% layer_predict() # %>% layer_naomit()
   if (inherits(trainer, "quantile_reg")) {
-    # add all levels to the forecaster and update postprocessor
-    tau <- sort(compare_quantile_args(
-      args_list$levels,
-      rlang::eval_tidy(trainer$args$tau)
+    # add all quantile_level to the forecaster and update postprocessor
+    quantile_levels <- sort(compare_quantile_args(
+      args_list$quantile_levels,
+      rlang::eval_tidy(trainer$args$quantile_levels)
     ))
-    args_list$levels <- tau
-    trainer$args$tau <- rlang::enquo(tau)
-    f <- layer_quantile_distn(f, levels = tau) %>% layer_point_from_distn()
+    args_list$quantile_levels <- quantile_levels
+    trainer$args$quantile_levels <- rlang::enquo(quantile_levels)
+    f <- layer_quantile_distn(f, quantile_levels = quantile_levels) %>%
+      layer_point_from_distn()
   } else {
     f <- layer_residual_quantiles(
       f,
-      probs = args_list$levels, symmetrize = args_list$symmetrize,
+      quantile_levels = args_list$quantile_levels,
+      symmetrize = args_list$symmetrize,
       by_key = args_list$quantile_by_key
     )
   }
@@ -173,7 +175,7 @@ arx_fcast_epi_workflow <- function(
 #'   The default `NULL` will attempt to determine this automatically.
 #' @param target_date Date. The date for which the forecast is intended.
 #'   The default `NULL` will attempt to determine this automatically.
-#' @param levels Vector or `NULL`. A vector of probabilities to produce
+#' @param quantile_levels Vector or `NULL`. A vector of probabilities to produce
 #'   prediction intervals. These are created by computing the quantiles of
 #'   training residuals. A `NULL` value will result in point forecasts only.
 #' @param symmetrize Logical. The default `TRUE` calculates
@@ -197,6 +199,7 @@ arx_fcast_epi_workflow <- function(
 #'   create a prediction. For this reason, setting `nafill_buffer < min(lags)`
 #'   will be treated as _additional_ allowed recent data rather than the
 #'   total amount of recent data to examine.
+#' @param ... Space to handle future expansions (unused).
 #'
 #'
 #' @return A list containing updated parameter choices with class `arx_flist`.
@@ -205,18 +208,19 @@ arx_fcast_epi_workflow <- function(
 #' @examples
 #' arx_args_list()
 #' arx_args_list(symmetrize = FALSE)
-#' arx_args_list(levels = c(.1, .3, .7, .9), n_training = 120)
+#' arx_args_list(quantile_levels = c(.1, .3, .7, .9), n_training = 120)
 arx_args_list <- function(
     lags = c(0L, 7L, 14L),
     ahead = 7L,
     n_training = Inf,
     forecast_date = NULL,
     target_date = NULL,
-    levels = c(0.05, 0.95),
+    quantile_levels = c(0.05, 0.95),
     symmetrize = TRUE,
     nonneg = TRUE,
     quantile_by_key = character(0L),
-    nafill_buffer = Inf) {
+    nafill_buffer = Inf,
+    ...) {
   # error checking if lags is a list
   .lags <- lags
   if (is.list(lags)) lags <- unlist(lags)
@@ -227,7 +231,7 @@ arx_args_list <- function(
   arg_is_date(forecast_date, target_date, allow_null = TRUE)
   arg_is_nonneg_int(ahead, lags)
   arg_is_lgl(symmetrize, nonneg)
-  arg_is_probabilities(levels, allow_null = TRUE)
+  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(nafill_buffer)) arg_is_pos_int(nafill_buffer, allow_null = TRUE)
@@ -238,7 +242,7 @@ arx_args_list <- function(
       lags = .lags,
       ahead,
       n_training,
-      levels,
+      quantile_levels,
       forecast_date,
       target_date,
       symmetrize,
@@ -259,8 +263,8 @@ print.arx_fcast <- function(x, ...) {
 }
 
 compare_quantile_args <- function(alist, tlist) {
-  default_alist <- eval(formals(arx_args_list)$levels)
-  default_tlist <- eval(formals(quantile_reg)$tau)
+  default_alist <- eval(formals(arx_args_list)$quantile_level)
+  default_tlist <- eval(formals(quantile_reg)$quantile_level)
   if (setequal(alist, default_alist)) {
     if (setequal(tlist, default_tlist)) {
       return(sort(unique(union(alist, tlist))))
diff --git a/R/canned-epipred.R b/R/canned-epipred.R
index bf99d74c7..7458655e8 100644
--- a/R/canned-epipred.R
+++ b/R/canned-epipred.R
@@ -8,13 +8,13 @@ validate_forecaster_inputs <- function(epi_data, outcome, predictors) {
   arg_is_chr(predictors)
   arg_is_chr_scalar(outcome)
   if (!outcome %in% names(epi_data)) {
-    cli::cli_abort("{outcome} was not found in the training data.")
+    cli::cli_abort("{.var {outcome}} was not found in the training data.")
   }
   check <- hardhat::check_column_names(epi_data, predictors)
   if (!check$ok) {
     cli::cli_abort(c(
       "At least one predictor was not found in the training data.",
-      "!" = "The following required columns are missing: {check$missing_names}."
+      "!" = "The following required columns are missing: {.val {check$missing_names}}."
     ))
   }
   invisible(TRUE)
@@ -41,8 +41,8 @@ arx_lags_validator <- function(predictors, lags) {
         predictors_miss <- setdiff(predictors, names(lags))
         cli::cli_abort(c(
           "If lags is a named list, then all predictors must be present.",
-          i = "The predictors are '{predictors}'.",
-          i = "So lags is missing '{predictors_miss}'."
+          i = "The predictors are {.var {predictors}}.",
+          i = "So lags is missing {.var {predictors_miss}}'."
         ))
       }
     }
diff --git a/R/dist_quantiles.R b/R/dist_quantiles.R
index ff14d6733..750e9560d 100644
--- a/R/dist_quantiles.R
+++ b/R/dist_quantiles.R
@@ -1,24 +1,28 @@
-new_quantiles <- function(q = double(), tau = double()) {
-  arg_is_probabilities(tau)
-
-  vec_cast(q, double())
-  vec_cast(tau, double())
-  stopifnot(length(q) == length(tau))
-  stopifnot(!vec_duplicate_any(tau))
-  if (is.unsorted(tau)) {
-    o <- vec_order(tau)
-    q <- q[o]
-    tau <- tau[o]
+#' @importFrom vctrs field vec_cast new_rcrd
+new_quantiles <- function(values = double(), quantile_levels = double()) {
+  arg_is_probabilities(quantile_levels)
+
+  vec_cast(values, double())
+  vec_cast(quantile_levels, double())
+  stopifnot(length(values) == length(quantile_levels))
+  stopifnot(!vctrs::vec_duplicate_any(quantile_levels))
+  if (is.unsorted(quantile_levels)) {
+    o <- vctrs::vec_order(quantile_levels)
+    values <- values[o]
+    quantile_levels <- quantile_levels[o]
   }
-  if (is.unsorted(q, na.rm = TRUE)) {
-    rlang::abort("`q[order(tau)]` produces unsorted quantiles.")
+  if (is.unsorted(values, na.rm = TRUE)) {
+    cli::cli_abort("`values[order(quantile_levels)]` produces unsorted quantiles.")
   }
 
-  new_rcrd(list(q = q, tau = tau),
+  new_rcrd(list(values = values, quantile_levels = quantile_levels),
     class = c("dist_quantiles", "dist_default")
   )
 }
 
+
+
+#' @importFrom vctrs vec_ptype_abbr vec_ptype_full
 #' @export
 vec_ptype_abbr.dist_quantiles <- function(x, ...) "dist_qntls"
 #' @export
@@ -26,46 +30,67 @@ vec_ptype_full.dist_quantiles <- function(x, ...) "dist_quantiles"
 
 #' @export
 format.dist_quantiles <- function(x, digits = 2, ...) {
-  q <- field(x, "q")
-  tau <- field(x, "tau")
-  rng <- range(tau, na.rm = TRUE)
-  paste0("[", round(rng[1], digits), ", ", round(rng[2], digits), "]<q-rng>")
+  m <- suppressWarnings(median(x))
+  paste0("quantiles(", round(m, digits), ")[", vctrs::vec_size(x), "]")
 }
 
 
-
-
 #' A distribution parameterized by a set of quantiles
 #'
-#' @param x A vector of values
-#' @param tau A vector of probabilities corresponding to `x`
+#' @param values A vector of values
+#' @param quantile_levels A vector of probabilities corresponding to `values`
 #'
 #' @export
 #'
-#' @import vctrs
 #' @examples
 #' dstn <- dist_quantiles(list(1:4, 8:11), list(c(.2, .4, .6, .8)))
 #' quantile(dstn, p = c(.1, .25, .5, .9))
 #' median(dstn)
 #'
 #' # it's a bit annoying to inspect the data
-#' vctrs::vec_data(vctrs::vec_data(dstn[1])[[1]])
-dist_quantiles <- function(x, tau) {
-  if (!is.list(x)) x <- list(x)
-  if (!is.list(tau)) tau <- list(tau)
-
-  x <- as_list_of(x, .ptype = double())
-  tau <- as_list_of(tau, .ptype = double())
-  args <- vec_recycle_common(x = x, tau = tau)
-  qntls <- as_list_of(map2(args$x, args$tau, new_quantiles))
+#' distributional::parameters(dstn[1])
+#' nested_quantiles(dstn[1])[[1]]
+#'
+#' dist_quantiles(1:4, 1:4 / 5)
+#' @importFrom vctrs as_list_of vec_recycle_common new_vctr
+dist_quantiles <- function(values, quantile_levels) {
+  if (!is.list(values)) values <- list(values)
+  if (!is.list(quantile_levels)) quantile_levels <- list(quantile_levels)
+
+  values <- as_list_of(values, .ptype = double())
+  quantile_levels <- as_list_of(quantile_levels, .ptype = double())
+  args <- vec_recycle_common(values = values, quantile_levels = quantile_levels)
+  qntls <- as_list_of(map2(args$values, args$quantile_levels, new_quantiles))
   new_vctr(qntls, class = "distribution")
 }
 
+validate_dist_quantiles <- function(values, quantile_levels) {
+  map(quantile_levels, arg_is_probabilities)
+  common_length <- vctrs::vec_size_common( # aborts internally
+    values = values,
+    quantile_levels = quantile_levels
+  )
+  length_diff <- vctrs::list_sizes(values) != vctrs::list_sizes(quantile_levels)
+  if (any(length_diff)) {
+    cli::cli_abort(c(
+      "`values` and `quantile_levels` must have common length.",
+      i = "Mismatches found at position(s): {.val {which(length_diff)}}."
+    ))
+  }
+  level_duplication <- map_lgl(quantile_levels, vctrs::vec_duplicate_any)
+  if (any(level_duplication)) {
+    cli::cli_abort(c(
+      "`quantile_levels` must not be duplicated.",
+      i = "Duplicates found at position(s): {.val {which(level_duplication)}}."
+    ))
+  }
+}
+
 
 #' Summarize a distribution with a set of quantiles
 #'
 #' @param x a `distribution` vector
-#' @param p a vector of probabilities at which to calculate quantiles
+#' @param probs a vector of probabilities at which to calculate quantiles
 #' @param ... additional arguments passed on to the `quantile` method
 #'
 #' @return a `distribution` vector containing `dist_quantiles`
@@ -74,80 +99,95 @@ dist_quantiles <- function(x, tau) {
 #' @examples
 #' library(distributional)
 #' dstn <- dist_normal(c(10, 2), c(5, 10))
-#' extrapolate_quantiles(dstn, p = c(.25, 0.5, .75))
+#' extrapolate_quantiles(dstn, probs = c(.25, 0.5, .75))
 #'
 #' dstn <- dist_quantiles(list(1:4, 8:11), list(c(.2, .4, .6, .8)))
 #' # because this distribution is already quantiles, any extra quantiles are
 #' # appended
-#' extrapolate_quantiles(dstn, p = c(.25, 0.5, .75))
+#' extrapolate_quantiles(dstn, probs = c(.25, 0.5, .75))
 #'
 #' dstn <- c(
 #'   dist_normal(c(10, 2), c(5, 10)),
 #'   dist_quantiles(list(1:4, 8:11), list(c(.2, .4, .6, .8)))
 #' )
-#' extrapolate_quantiles(dstn, p = c(.25, 0.5, .75))
-extrapolate_quantiles <- function(x, p, ...) {
+#' extrapolate_quantiles(dstn, probs = c(.25, 0.5, .75))
+extrapolate_quantiles <- function(x, probs, ...) {
   UseMethod("extrapolate_quantiles")
 }
 
 #' @export
-extrapolate_quantiles.distribution <- function(x, p, ...) {
-  arg_is_probabilities(p)
-  dstn <- lapply(vec_data(x), extrapolate_quantiles, p = p, ...)
-  distributional:::wrap_dist(dstn)
+#' @importFrom vctrs vec_data
+extrapolate_quantiles.distribution <- function(x, probs, ...) {
+  arg_is_probabilities(probs)
+  dstn <- lapply(vec_data(x), extrapolate_quantiles, probs = probs, ...)
+  new_vctr(dstn, vars = NULL, class = "distribution")
 }
 
 #' @export
-extrapolate_quantiles.dist_default <- function(x, p, ...) {
-  q <- quantile(x, p, ...)
-  new_quantiles(q = q, tau = p)
+extrapolate_quantiles.dist_default <- function(x, probs, ...) {
+  values <- quantile(x, probs, ...)
+  new_quantiles(values = values, quantile_levels = probs)
 }
 
 #' @export
-extrapolate_quantiles.dist_quantiles <- function(x, p, ...) {
-  q <- quantile(x, p, ...)
-  tau <- field(x, "tau")
-  qvals <- field(x, "q")
-  new_quantiles(q = c(qvals, q), tau = c(tau, p))
+extrapolate_quantiles.dist_quantiles <- function(x, probs, ...) {
+  new_values <- quantile(x, probs, ...)
+  quantile_levels <- field(x, "quantile_levels")
+  values <- field(x, "values")
+  new_quantiles(
+    values = c(values, new_values),
+    quantile_levels = c(quantile_levels, probs)
+  )
 }
 
 is_dist_quantiles <- function(x) {
-  is_distribution(x) && all(stats::family(x) == "quantiles")
+  is_distribution(x) & all(stats::family(x) == "quantiles")
 }
 
 
 
-
 #' @export
 #' @importFrom stats median qnorm family
-median.dist_quantiles <- function(x, ..., middle = c("cubic", "linear")) {
+median.dist_quantiles <- function(x, na.rm = FALSE, ..., middle = c("cubic", "linear")) {
+  quantile_levels <- field(x, "quantile_levels")
+  values <- field(x, "values")
+  if (0.5 %in% quantile_levels) {
+    return(values[match(0.5, quantile_levels)])
+  }
+  if (length(quantile_levels) < 2 || min(quantile_levels) > 0.5 || max(quantile_levels) < 0.5) {
+    return(NA)
+  }
+  if (length(quantile_levels) < 3 || min(quantile_levels) > .25 || max(quantile_levels) < .75) {
+    return(stats::approx(quantile_levels, values, xout = 0.5)$y)
+  }
   quantile(x, 0.5, ..., middle = middle)
 }
 
 # placeholder to avoid errors, but not ideal
 #' @export
-mean.dist_quantiles <- function(x, ..., middle = c("cubic", "linear")) {
+mean.dist_quantiles <- function(x, na.rm = FALSE, ..., middle = c("cubic", "linear")) {
   median(x, ..., middle = middle)
 }
 
 #' @export
 #' @importFrom stats quantile
 #' @import distributional
-quantile.dist_quantiles <- function(x, probs, ...,
-                                    middle = c("cubic", "linear"),
-                                    left_tail = c("normal", "exponential"),
-                                    right_tail = c("normal", "exponential")) {
-  arg_is_probabilities(probs)
+quantile.dist_quantiles <- function(
+    x, p, ...,
+    middle = c("cubic", "linear"),
+    left_tail = c("normal", "exponential"),
+    right_tail = c("normal", "exponential")) {
+  arg_is_probabilities(p)
   middle <- match.arg(middle)
   left_tail <- match.arg(left_tail)
   right_tail <- match.arg(right_tail)
-  quantile_extrapolate(x, probs, middle, left_tail, right_tail)
+  quantile_extrapolate(x, p, middle, left_tail, right_tail)
 }
 
 
 quantile_extrapolate <- function(x, tau_out, middle, left_tail, right_tail) {
-  tau <- field(x, "tau")
-  qvals <- field(x, "q")
+  tau <- field(x, "quantile_levels")
+  qvals <- field(x, "values")
   r <- range(tau, na.rm = TRUE)
   qvals_out <- rep(NA, length(tau_out))
 
@@ -157,7 +197,7 @@ quantile_extrapolate <- function(x, tau_out, middle, left_tail, right_tail) {
     return(qvals[match(tau_out, tau)])
   }
   if (length(qvals) < 3 || r[1] > .25 || r[2] < .75) {
-    rlang::warn(c(
+    cli::cli_warn(c(
       "Quantile extrapolation is not possible with fewer than",
       "3 quantiles or when the probs don't span [.25, .75]"
     ))
@@ -258,10 +298,10 @@ norm_tail_q <- function(p, q, target) {
 #' @method Math dist_quantiles
 #' @export
 Math.dist_quantiles <- function(x, ...) {
-  tau <- field(x, "tau")
-  q <- field(x, "q")
-  q <- vctrs::vec_math(.Generic, q, ...)
-  new_quantiles(q = q, tau = tau)
+  quantile_levels <- field(x, "quantile_levels")
+  values <- field(x, "values")
+  values <- vctrs::vec_math(.Generic, values, ...)
+  new_quantiles(values = values, quantile_levels = quantile_levels)
 }
 
 #' @method Ops dist_quantiles
@@ -274,16 +314,16 @@ Ops.dist_quantiles <- function(e1, e2) {
   is_dist <- c(inherits(e1, "dist_default"), inherits(e2, "dist_default"))
   tau1 <- tau2 <- NULL
   if (is_quantiles[1]) {
-    q1 <- field(e1, "q")
-    tau1 <- field(e1, "tau")
+    q1 <- field(e1, "values")
+    tau1 <- field(e1, "quantile_levels")
   }
   if (is_quantiles[2]) {
-    q2 <- field(e2, "q")
-    tau2 <- field(e2, "tau")
+    q2 <- field(e2, "values")
+    tau2 <- field(e2, "quantile_levels")
   }
   tau <- union(tau1, tau2)
   if (all(is_dist)) {
-    rlang::abort(
+    cli::cli_abort(
       "You can't perform arithmetic between two distributions like this."
     )
   } else {
@@ -294,18 +334,18 @@ Ops.dist_quantiles <- function(e1, e2) {
     }
   }
   q <- vctrs::vec_arith(.Generic, q1, q2)
-  new_quantiles(q = q, tau = tau)
+  new_quantiles(values = q, quantile_levels = tau)
 }
 
 #' @method is.na distribution
 #' @export
 is.na.distribution <- function(x) {
-  sapply(vctrs::vec_data(x), is.na)
+  sapply(vec_data(x), is.na)
 }
 
 #' @method is.na dist_quantiles
 #' @export
 is.na.dist_quantiles <- function(x) {
-  q <- field(x, "q")
+  q <- field(x, "values")
   all(is.na(q))
 }
diff --git a/R/epipredict-package.R b/R/epipredict-package.R
index da4991feb..11e2ec833 100644
--- a/R/epipredict-package.R
+++ b/R/epipredict-package.R
@@ -2,6 +2,7 @@
 #' @importFrom tibble tibble
 #' @importFrom rlang := !!
 #' @importFrom stats poly predict lm residuals quantile
+#' @importFrom cli cli_abort
 #' @importFrom lifecycle deprecated
 #' @import epiprocess parsnip
 ## usethis namespace: end
diff --git a/R/extract.R b/R/extract.R
index 574cc40cc..e227b59b1 100644
--- a/R/extract.R
+++ b/R/extract.R
@@ -13,7 +13,7 @@
 #' @examples
 #' f <- frosting() %>%
 #'   layer_predict() %>%
-#'   layer_residual_quantiles(probs = c(0.0275, 0.975), symmetrize = FALSE) %>%
+#'   layer_residual_quantiles(quantile_levels = c(0.0275, 0.975), symmetrize = FALSE) %>%
 #'   layer_naomit(.pred)
 #'
 #' extract_argument(f, "layer_residual_quantiles", "symmetrize")
@@ -27,10 +27,10 @@ extract_argument.layer <- function(x, name, arg, ...) {
   arg_is_chr_scalar(name, arg)
   in_layer_name <- class(x)[1]
   if (name != in_layer_name) {
-    cli_stop("Requested {name} not found. This is a(n) {in_layer_name}.")
+    cli_abort("Requested {.val {name}} not found. This is a(n) {.cls {in_layer_name}}.")
   }
   if (!arg %in% names(x)) {
-    cli_stop("Requested argument {arg} not found in {name}.")
+    cli_abort("Requested argument {.val {arg}} not found in {.val {name}}.")
   }
   x[[arg]]
 }
@@ -41,10 +41,10 @@ extract_argument.step <- function(x, name, arg, ...) {
   arg_is_chr_scalar(name, arg)
   in_step_name <- class(x)[1]
   if (name != in_step_name) {
-    cli_stop("Requested {name} not found. This is a {in_step_name}.")
+    cli_abort("Requested {.val {name}} not found. This is a {.cls {in_step_name}}.")
   }
   if (!arg %in% names(x)) {
-    cli_stop("Requested argument {arg} not found in {name}.")
+    cli_abort("Requested argument {.val {arg}} not found in {.val {name}}.")
   }
   x[[arg]]
 }
@@ -55,7 +55,7 @@ extract_argument.recipe <- function(x, name, arg, ...) {
   step_names <- map_chr(x$steps, ~ class(.x)[1])
   has_step <- name %in% step_names
   if (!has_step) {
-    cli_stop("recipe object does not contain a {name}.")
+    cli_abort("recipe object does not contain a {.val {name}}.")
   }
   step_locations <- which(name == step_names)
   out <- map(x$steps[step_locations], extract_argument, name = name, arg = arg)
@@ -69,7 +69,7 @@ extract_argument.frosting <- function(x, name, arg, ...) {
   layer_names <- map_chr(x$layers, ~ class(.x)[1])
   has_layer <- name %in% layer_names
   if (!has_layer) {
-    cli_stop("frosting object does not contain a {name} layer.")
+    cli_abort("frosting object does not contain a {.val {name}}.")
   }
   layer_locations <- which(name == layer_names)
   out <- map(x$layers[layer_locations], extract_argument, name = name, arg = arg)
@@ -83,7 +83,7 @@ extract_argument.epi_workflow <- function(x, name, arg, ...) {
   type <- sub("_.*", "", name)
   if (type %in% c("check", "step")) {
     if (!workflows:::has_preprocessor_recipe(x)) {
-      cli_stop("The workflow must have a recipe preprocessor.")
+      cli_abort("The workflow must have a recipe preprocessor.")
     }
     out <- extract_argument(x$pre$actions$recipe$recipe, name, arg)
   }
@@ -91,7 +91,7 @@ extract_argument.epi_workflow <- function(x, name, arg, ...) {
     out <- extract_argument(extract_frosting(x), name, arg)
   }
   if (!type %in% c("check", "step", "layer")) {
-    cli_stop("{name} must begin with one of step, check, or layer")
+    cli_abort("{.val {name}} must begin with one of step, check, or layer")
   }
   return(out)
 }
diff --git a/R/flatline_forecaster.R b/R/flatline_forecaster.R
index 197c8cca5..9342bd740 100644
--- a/R/flatline_forecaster.R
+++ b/R/flatline_forecaster.R
@@ -58,7 +58,7 @@ flatline_forecaster <- function(
   f <- frosting() %>%
     layer_predict() %>%
     layer_residual_quantiles(
-      probs = args_list$levels,
+      quantile_levels = args_list$quantile_levels,
       symmetrize = args_list$symmetrize,
       by_key = args_list$quantile_by_key
     ) %>%
@@ -107,24 +107,25 @@ flatline_forecaster <- function(
 #' @examples
 #' flatline_args_list()
 #' flatline_args_list(symmetrize = FALSE)
-#' flatline_args_list(levels = c(.1, .3, .7, .9), n_training = 120)
+#' 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,
-    levels = c(0.05, 0.95),
+    quantile_levels = c(0.05, 0.95),
     symmetrize = TRUE,
     nonneg = TRUE,
     quantile_by_key = character(0L),
-    nafill_buffer = Inf) {
+    nafill_buffer = Inf,
+    ...) {
   arg_is_scalar(ahead, n_training)
   arg_is_chr(quantile_by_key, allow_empty = TRUE)
   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)
   arg_is_lgl(symmetrize, nonneg)
-  arg_is_probabilities(levels, allow_null = TRUE)
+  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(nafill_buffer)) arg_is_pos_int(nafill_buffer, allow_null = TRUE)
@@ -135,7 +136,7 @@ flatline_args_list <- function(
       n_training,
       forecast_date,
       target_date,
-      levels,
+      quantile_levels,
       symmetrize,
       nonneg,
       quantile_by_key,
diff --git a/R/layer_point_from_distn.R b/R/layer_point_from_distn.R
index 9c7b0eb3e..8aece79e0 100644
--- a/R/layer_point_from_distn.R
+++ b/R/layer_point_from_distn.R
@@ -24,7 +24,7 @@
 #'   step_epi_ahead(death_rate, ahead = 7) %>%
 #'   step_epi_naomit()
 #'
-#' wf <- epi_workflow(r, quantile_reg(tau = c(.25, .5, .75))) %>% fit(jhu)
+#' wf <- epi_workflow(r, quantile_reg(quantile_levels = c(.25, .5, .75))) %>% fit(jhu)
 #'
 #' latest <- get_test_data(recipe = r, x = jhu)
 #'
diff --git a/R/layer_quantile_distn.R b/R/layer_quantile_distn.R
index 2b63206b2..a99eed326 100644
--- a/R/layer_quantile_distn.R
+++ b/R/layer_quantile_distn.R
@@ -3,11 +3,11 @@
 #' This function calculates quantiles when the prediction was _distributional_.
 #' Currently, the only distributional engine is `quantile_reg()`.
 #' If this engine is used, then this layer will grab out estimated (or extrapolated)
-#' quantiles at the requested levels.
+#' quantiles at the requested quantile values.
 #'
 #' @param frosting a `frosting` postprocessor
 #' @param ... Unused, include for consistency with other layers.
-#' @param levels a vector of probabilities (quantiles) to extract
+#' @param quantile_levels a vector of probabilities to extract
 #' @param truncate Do we truncate the distribution to an interval
 #' @param name character. The name for the output column.
 #' @param id a random id string
@@ -25,7 +25,8 @@
 #'   step_epi_ahead(death_rate, ahead = 7) %>%
 #'   step_epi_naomit()
 #'
-#' wf <- epi_workflow(r, quantile_reg(tau = c(.25, .5, .75))) %>% fit(jhu)
+#' wf <- epi_workflow(r, quantile_reg(quantile_levels = c(.25, .5, .75))) %>%
+#'   fit(jhu)
 #'
 #' latest <- get_test_data(recipe = r, x = jhu)
 #'
@@ -39,13 +40,13 @@
 #' p
 layer_quantile_distn <- function(frosting,
                                  ...,
-                                 levels = c(.25, .75),
+                                 quantile_levels = c(.25, .75),
                                  truncate = c(-Inf, Inf),
                                  name = ".pred_distn",
                                  id = rand_id("quantile_distn")) {
   rlang::check_dots_empty()
   arg_is_chr_scalar(name, id)
-  arg_is_probabilities(levels)
+  arg_is_probabilities(quantile_levels)
   stopifnot(
     length(truncate) == 2L, is.numeric(truncate), truncate[1] < truncate[2]
   )
@@ -53,7 +54,7 @@ layer_quantile_distn <- function(frosting,
   add_layer(
     frosting,
     layer_quantile_distn_new(
-      levels = levels,
+      quantile_levels = quantile_levels,
       truncate = truncate,
       name = name,
       id = id
@@ -61,9 +62,9 @@ layer_quantile_distn <- function(frosting,
   )
 }
 
-layer_quantile_distn_new <- function(levels, truncate, name, id) {
+layer_quantile_distn_new <- function(quantile_levels, truncate, name, id) {
   layer("quantile_distn",
-    levels = levels,
+    quantile_levels = quantile_levels,
     truncate = truncate,
     name = name,
     id = id
@@ -75,14 +76,15 @@ slather.layer_quantile_distn <-
   function(object, components, workflow, new_data, ...) {
     dstn <- components$predictions$.pred
     if (!inherits(dstn, "distribution")) {
-      rlang::abort(
-        c(
-          "`layer_quantile_distn` requires distributional predictions.",
-          "These are of class {class(dstn)}."
-        )
-      )
+      cli_abort(c(
+        "`layer_quantile_distn()` requires distributional predictions.",
+        "These are of class {.cls {class(dstn)}}."
+      ))
     }
-    dstn <- dist_quantiles(quantile(dstn, object$levels), object["levels"])
+    dstn <- dist_quantiles(
+      quantile(dstn, object$quantile_levels),
+      object$quantile_levels
+    )
 
     truncate <- object$truncate
     if (!all(is.infinite(truncate))) {
@@ -100,9 +102,9 @@ print.layer_quantile_distn <- function(
   title <- "Creating predictive quantiles"
   td <- "<calculated>"
   td <- rlang::enquos(td)
-  ext <- x$levels
+  ext <- x$quantile_levels
   print_layer(td,
-    title = title, width = width, conjunction = "levels",
+    title = title, width = width, conjunction = "quantile_levels",
     extra_text = ext
   )
 }
diff --git a/R/layer_residual_quantiles.R b/R/layer_residual_quantiles.R
index b9a71e265..932f73246 100644
--- a/R/layer_residual_quantiles.R
+++ b/R/layer_residual_quantiles.R
@@ -2,7 +2,7 @@
 #'
 #' @param frosting a `frosting` postprocessor
 #' @param ... Unused, include for consistency with other layers.
-#' @param probs numeric vector of probabilities with values in (0,1)
+#' @param quantile_levels numeric vector of probabilities with values in (0,1)
 #'   referring to the desired quantile.
 #' @param symmetrize logical. If `TRUE` then interval will be symmetric.
 #' @param by_key A character vector of keys to group the residuals by before
@@ -28,7 +28,7 @@
 #'
 #' f <- frosting() %>%
 #'   layer_predict() %>%
-#'   layer_residual_quantiles(probs = c(0.0275, 0.975), symmetrize = FALSE) %>%
+#'   layer_residual_quantiles(quantile_levels = c(0.0275, 0.975), symmetrize = FALSE) %>%
 #'   layer_naomit(.pred)
 #' wf1 <- wf %>% add_frosting(f)
 #'
@@ -36,27 +36,28 @@
 #'
 #' f2 <- frosting() %>%
 #'   layer_predict() %>%
-#'   layer_residual_quantiles(probs = c(0.3, 0.7), by_key = "geo_value") %>%
+#'   layer_residual_quantiles(quantile_levels = c(0.3, 0.7), by_key = "geo_value") %>%
 #'   layer_naomit(.pred)
 #' wf2 <- wf %>% add_frosting(f2)
 #'
 #' p2 <- predict(wf2, latest)
-layer_residual_quantiles <- function(frosting, ...,
-                                     probs = c(0.05, 0.95),
-                                     symmetrize = TRUE,
-                                     by_key = character(0L),
-                                     name = ".pred_distn",
-                                     id = rand_id("residual_quantiles")) {
+layer_residual_quantiles <- function(
+    frosting, ...,
+    quantile_levels = c(0.05, 0.95),
+    symmetrize = TRUE,
+    by_key = character(0L),
+    name = ".pred_distn",
+    id = rand_id("residual_quantiles")) {
   rlang::check_dots_empty()
   arg_is_scalar(symmetrize)
   arg_is_chr_scalar(name, id)
   arg_is_chr(by_key, allow_empty = TRUE)
-  arg_is_probabilities(probs)
+  arg_is_probabilities(quantile_levels)
   arg_is_lgl(symmetrize)
   add_layer(
     frosting,
     layer_residual_quantiles_new(
-      probs = probs,
+      quantile_levels = quantile_levels,
       symmetrize = symmetrize,
       by_key = by_key,
       name = name,
@@ -65,9 +66,10 @@ layer_residual_quantiles <- function(frosting, ...,
   )
 }
 
-layer_residual_quantiles_new <- function(probs, symmetrize, by_key, name, id) {
+layer_residual_quantiles_new <- function(
+    quantile_levels, symmetrize, by_key, name, id) {
   layer("residual_quantiles",
-    probs = probs, symmetrize = symmetrize,
+    quantile_levels = quantile_levels, symmetrize = symmetrize,
     by_key = by_key, name = name, id = id
   )
 }
@@ -77,7 +79,7 @@ slather.layer_residual_quantiles <-
   function(object, components, workflow, new_data, ...) {
     the_fit <- workflows::extract_fit_parsnip(workflow)
 
-    if (is.null(object$probs)) {
+    if (is.null(object$quantile_levels)) {
       return(components)
     }
 
@@ -93,19 +95,19 @@ slather.layer_residual_quantiles <-
       common <- intersect(object$by_key, names(key_cols))
       excess <- setdiff(object$by_key, names(key_cols))
       if (length(excess) > 0L) {
-        rlang::warn(
-          "Requested residual grouping key(s) {excess} are unavailable ",
-          "in the original data. Grouping by the remainder: {common}."
-        )
+        cli::cli_warn(c(
+          "Requested residual grouping key(s) {.val {excess}} are unavailable ",
+          "in the original data. Grouping by the remainder: {.val {common}}."
+        ))
       }
       if (length(common) > 0L) {
         r <- r %>% dplyr::select(tidyselect::any_of(c(common, ".resid")))
         common_in_r <- common[common %in% names(r)]
         if (length(common_in_r) != length(common)) {
-          rlang::warn(
+          cli::cli_warn(c(
             "Some grouping keys are not in data.frame returned by the",
             "`residuals()` method. Groupings may not be correct."
-          )
+          ))
         }
         r <- dplyr::bind_cols(key_cols, r) %>%
           dplyr::group_by(!!!rlang::syms(common))
@@ -114,15 +116,15 @@ slather.layer_residual_quantiles <-
 
     r <- r %>%
       dplyr::summarize(
-        q = list(quantile(
+        dstn = list(quantile(
           c(.resid, s * .resid),
-          probs = object$probs, na.rm = TRUE
+          probs = object$quantile_levels, na.rm = TRUE
         ))
       )
 
     estimate <- components$predictions$.pred
     res <- tibble::tibble(
-      .pred_distn = dist_quantiles(map2(estimate, r$q, "+"), object$probs)
+      .pred_distn = dist_quantiles(map2(estimate, r$dstn, "+"), object$quantile_levels)
     )
     res <- check_pname(res, components$predictions, object)
     components$predictions <- dplyr::mutate(components$predictions, !!!res)
@@ -142,7 +144,7 @@ grab_residuals <- function(the_fit, components) {
         return(r)
       } else { # failure
         cli::cli_warn(c(
-          "The `residuals()` method for {.cls cl} objects results in",
+          "The `residuals()` method for objects of class {.cls {cl}} results in",
           "a data frame without a column named `.resid`.",
           i = "Residual quantiles will be calculated directly from the",
           i = "difference between predictions and observations.",
@@ -153,7 +155,7 @@ grab_residuals <- function(the_fit, components) {
       return(tibble(.resid = drop(r)))
     } else { # failure
       cli::cli_warn(c(
-        "The `residuals()` method for {.cls cl} objects results in an",
+        "The `residuals()` method for objects of class {.cls {cl}} results in an",
         "object that is neither a data frame with a column named `.resid`,",
         "nor something coercible to a vector.",
         i = "Residual quantiles will be calculated directly from the",
@@ -176,9 +178,9 @@ print.layer_residual_quantiles <- function(
   title <- "Resampling residuals for predictive quantiles"
   td <- "<calculated>"
   td <- rlang::enquos(td)
-  ext <- x$probs
+  ext <- x$quantile_levels
   print_layer(td,
-    title = title, width = width, conjunction = "levels",
+    title = title, width = width, conjunction = "quantile_levels",
     extra_text = ext
   )
 }
diff --git a/R/layer_threshold_preds.R b/R/layer_threshold_preds.R
index 4107504a9..957ac2419 100644
--- a/R/layer_threshold_preds.R
+++ b/R/layer_threshold_preds.R
@@ -91,10 +91,10 @@ snap.dist_default <- function(x, lower, upper, ...) {
 
 #' @export
 snap.dist_quantiles <- function(x, lower, upper, ...) {
-  q <- field(x, "q")
-  tau <- field(x, "tau")
-  q <- snap(q, lower, upper)
-  new_quantiles(q = q, tau = tau)
+  values <- field(x, "values")
+  quantile_levels <- field(x, "quantile_levels")
+  values <- snap(values, lower, upper)
+  new_quantiles(values = values, quantile_levels = quantile_levels)
 }
 
 #' @export
diff --git a/R/make_quantile_reg.R b/R/make_quantile_reg.R
index eef4d4c97..bf17e7202 100644
--- a/R/make_quantile_reg.R
+++ b/R/make_quantile_reg.R
@@ -9,7 +9,7 @@
 #'   The only possible value for this model is "regression".
 #' @param engine Character string naming the fitting function. Currently, only
 #'   "rq" is supported.
-#' @param tau A scalar or vector of values in (0, 1) to determine which
+#' @param quantile_levels A scalar or vector of values in (0, 1) to determine which
 #'   quantiles to estimate (default is 0.5).
 #'
 #' @export
@@ -19,23 +19,23 @@
 #' @importFrom quantreg rq
 #' @examples
 #' tib <- data.frame(y = rnorm(100), x1 = rnorm(100), x2 = rnorm(100))
-#' rq_spec <- quantile_reg(tau = c(.2, .8)) %>% set_engine("rq")
+#' rq_spec <- quantile_reg(quantile_levels = c(.2, .8)) %>% set_engine("rq")
 #' ff <- rq_spec %>% fit(y ~ ., data = tib)
 #' predict(ff, new_data = tib)
-quantile_reg <- function(mode = "regression", engine = "rq", tau = 0.5) {
+quantile_reg <- function(mode = "regression", engine = "rq", quantile_levels = 0.5) {
   # Check for correct mode
   if (mode != "regression") {
-    rlang::abort("`mode` should be 'regression'")
+    cli_abort("`mode` must be 'regression'")
   }
 
   # Capture the arguments in quosures
-  if (any(tau > 1)) rlang::abort("All `tau` must be less than 1.")
-  if (any(tau < 0)) rlang::abort("All `tau` must be greater than 0.")
-  if (is.unsorted(tau)) {
-    rlang::warn("Sorting tau to increasing order.")
-    tau <- sort(tau)
+  if (any(quantile_levels > 1)) cli_abort("All `quantile_levels` must be less than 1.")
+  if (any(quantile_levels < 0)) cli_abort("All `quantile_levels` must be greater than 0.")
+  if (is.unsorted(quantile_levels)) {
+    cli::cli_warn("Sorting `quantile_levels` to increasing order.")
+    quantile_levels <- sort(quantile_levels)
   }
-  args <- list(tau = rlang::enquo(tau))
+  args <- list(quantile_levels = rlang::enquo(quantile_levels))
 
   # Save some empty slots for future parts of the specification
   parsnip::new_model_spec(
@@ -60,7 +60,7 @@ make_quantile_reg <- function() {
   parsnip::set_model_arg(
     model = "quantile_reg",
     eng = "rq",
-    parsnip = "tau",
+    parsnip = "quantile_levels",
     original = "tau",
     func = list(pkg = "quantreg", fun = "rq"),
     has_submodel = FALSE
@@ -101,13 +101,14 @@ make_quantile_reg <- function() {
 
     # can't make a method because object is second
     out <- switch(type,
-      rq = dist_quantiles(unname(as.list(x)), object$tau), # one quantile
+      rq = dist_quantiles(unname(as.list(x)), object$quantile_levels), # one quantile
       rqs = {
         x <- lapply(unname(split(x, seq(nrow(x)))), function(q) sort(q))
         dist_quantiles(x, list(object$tau))
       },
-      rlang::abort(c("Prediction not implemented for this `rq` type.",
-        i = "See `?quantreg::rq`."
+      cli_abort(c(
+        "Prediction is not implemented for this `rq` type.",
+        i = "See {.fun quantreg::rq}."
       ))
     )
     return(data.frame(.pred = out))
diff --git a/R/make_smooth_quantile_reg.R b/R/make_smooth_quantile_reg.R
index 7d25a8c6b..49b7b4e36 100644
--- a/R/make_smooth_quantile_reg.R
+++ b/R/make_smooth_quantile_reg.R
@@ -9,7 +9,7 @@
 #'   The only possible value for this model is "regression".
 #' @param engine Character string naming the fitting function. Currently, only
 #'   "smooth_qr" is supported.
-#' @param tau A scalar or vector of values in (0, 1) to determine which
+#' @param quantile_levels A scalar or vector of values in (0, 1) to determine which
 #'   quantiles to estimate (default is 0.5).
 #' @param outcome_locations Defaults to the vector `1:ncol(y)` but if the
 #'   responses are observed at a different spacing (or appear in a different
@@ -28,7 +28,7 @@
 #'   y4 = rnorm(100), y5 = rnorm(100), y6 = rnorm(100),
 #'   x1 = rnorm(100), x2 = rnorm(100)
 #' )
-#' qr_spec <- smooth_quantile_reg(tau = c(.2, .5, .8), outcome_locations = 1:6)
+#' qr_spec <- smooth_quantile_reg(quantile_levels = c(.2, .5, .8), outcome_locations = 1:6)
 #' ff <- qr_spec %>% fit(cbind(y1, y2, y3, y4, y5, y6) ~ ., data = tib)
 #' p <- predict(ff, new_data = tib)
 #'
@@ -37,7 +37,7 @@
 #' fd <- x[length(x) - 20]
 #' XY <- smoothqr::lagmat(y[1:(length(y) - 20)], c(-20:20))
 #' XY <- tibble::as_tibble(XY)
-#' qr_spec <- smooth_quantile_reg(tau = c(.2, .5, .8), outcome_locations = 20:1)
+#' qr_spec <- smooth_quantile_reg(quantile_levels = c(.2, .5, .8), outcome_locations = 20:1)
 #' tt <- qr_spec %>% fit_xy(x = XY[, 21:41], y = XY[, 1:20])
 #'
 #' library(tidyr)
@@ -54,7 +54,7 @@
 #'     x = x[length(x) - 20] + ahead / 100 * 2 * pi,
 #'     ahead = NULL
 #'   ) %>%
-#'   pivot_wider(names_from = tau, values_from = q)
+#'   pivot_wider(names_from = quantile_levels, values_from = values)
 #' plot(x, y, pch = 16, xlim = c(pi, 2 * pi), col = "lightgrey")
 #' curve(sin(x), add = TRUE)
 #' abline(v = fd, lty = 2)
@@ -77,23 +77,24 @@ smooth_quantile_reg <- function(
     mode = "regression",
     engine = "smoothqr",
     outcome_locations = NULL,
-    tau = 0.5,
+    quantile_levels = 0.5,
     degree = 3L) {
   # Check for correct mode
-  if (mode != "regression") rlang::abort("`mode` must be 'regression'")
-  if (engine != "smoothqr") rlang::abort("`engine` must be 'smoothqr'")
+  if (mode != "regression") cli_abort("`mode` must be 'regression'")
+  if (engine != "smoothqr") cli_abort("`engine` must be 'smoothqr'")
 
-  arg_is_probabilities(tau)
+  arg_is_probabilities(quantile_levels)
   arg_is_pos_int(degree)
   arg_is_scalar(degree)
   arg_is_numeric(outcome_locations, allow_null = TRUE)
-  if (is.unsorted(tau)) {
-    rlang::warn("Sorting tau to increasing order.")
-    tau <- sort(tau)
+  if (is.unsorted(quantile_levels)) {
+    rlang::warn("Sorting `quantile_levels` to increasing order.")
+    quantile_levels <- sort(quantile_levels)
   }
 
   args <- list(
-    tau = rlang::enquo(tau), degree = rlang::enquo(degree),
+    quantile_levels = rlang::enquo(quantile_levels),
+    degree = rlang::enquo(degree),
     outcome_locations = rlang::enquo(outcome_locations)
   )
 
@@ -123,7 +124,7 @@ make_smooth_quantile_reg <- function() {
   parsnip::set_model_arg(
     model = "smooth_quantile_reg",
     eng = "smoothqr",
-    parsnip = "tau",
+    parsnip = "quantile_levels",
     original = "tau",
     func = list(pkg = "smoothqr", fun = "smooth_qr"),
     has_submodel = FALSE
diff --git a/R/pivot_quantiles.R b/R/pivot_quantiles.R
index a156bcf90..e632748df 100644
--- a/R/pivot_quantiles.R
+++ b/R/pivot_quantiles.R
@@ -16,7 +16,7 @@ nested_quantiles <- function(x) {
   distributional:::dist_apply(x, .f = function(z) {
     tibble::as_tibble(vec_data(z)) %>%
       dplyr::mutate(dplyr::across(tidyselect::everything(), as.double)) %>%
-      list_of()
+      vctrs::list_of()
   })
 }
 
@@ -130,14 +130,14 @@ pivot_quantiles_wider <- function(.data, ...) {
       .data <- .data %>%
         tidyr::unnest(tidyselect::all_of(col)) %>%
         tidyr::pivot_wider(
-          names_from = "tau", values_from = "q",
+          names_from = "quantile_levels", values_from = "values",
           names_prefix = paste0(col, "_")
         )
     }
   } else {
     .data <- .data %>%
       tidyr::unnest(tidyselect::all_of(cols)) %>%
-      tidyr::pivot_wider(names_from = "tau", values_from = "q")
+      tidyr::pivot_wider(names_from = "quantile_levels", values_from = "values")
   }
   .data
 }
diff --git a/man/arx_args_list.Rd b/man/arx_args_list.Rd
index b4bc26f55..e5d2391c8 100644
--- a/man/arx_args_list.Rd
+++ b/man/arx_args_list.Rd
@@ -10,11 +10,12 @@ arx_args_list(
   n_training = Inf,
   forecast_date = NULL,
   target_date = NULL,
-  levels = c(0.05, 0.95),
+  quantile_levels = c(0.05, 0.95),
   symmetrize = TRUE,
   nonneg = TRUE,
   quantile_by_key = character(0L),
-  nafill_buffer = Inf
+  nafill_buffer = Inf,
+  ...
 )
 }
 \arguments{
@@ -35,7 +36,7 @@ The default \code{NULL} will attempt to determine this automatically.}
 \item{target_date}{Date. The date for which the forecast is intended.
 The default \code{NULL} will attempt to determine this automatically.}
 
-\item{levels}{Vector or \code{NULL}. A vector of probabilities to produce
+\item{quantile_levels}{Vector or \code{NULL}. A vector of probabilities to produce
 prediction intervals. These are created by computing the quantiles of
 training residuals. A \code{NULL} value will result in point forecasts only.}
 
@@ -63,6 +64,8 @@ we require at least \code{min(lags)} rows of recent data per \code{geo_value} to
 create a prediction. For this reason, setting \code{nafill_buffer < min(lags)}
 will be treated as \emph{additional} allowed recent data rather than the
 total amount of recent data to examine.}
+
+\item{...}{Space to handle future expansions (unused).}
 }
 \value{
 A list containing updated parameter choices with class \code{arx_flist}.
@@ -73,5 +76,5 @@ Constructs a list of arguments for \code{\link[=arx_forecaster]{arx_forecaster()
 \examples{
 arx_args_list()
 arx_args_list(symmetrize = FALSE)
-arx_args_list(levels = c(.1, .3, .7, .9), n_training = 120)
+arx_args_list(quantile_levels = c(.1, .3, .7, .9), n_training = 120)
 }
diff --git a/man/arx_class_args_list.Rd b/man/arx_class_args_list.Rd
index 2e96f0341..fa7a407f0 100644
--- a/man/arx_class_args_list.Rd
+++ b/man/arx_class_args_list.Rd
@@ -16,7 +16,8 @@ arx_class_args_list(
   method = c("rel_change", "linear_reg", "smooth_spline", "trend_filter"),
   log_scale = FALSE,
   additional_gr_args = list(),
-  nafill_buffer = Inf
+  nafill_buffer = Inf,
+  ...
 )
 }
 \arguments{
@@ -82,6 +83,8 @@ we require at least \code{min(lags)} rows of recent data per \code{geo_value} to
 create a prediction. For this reason, setting \code{nafill_buffer < min(lags)}
 will be treated as \emph{additional} allowed recent data rather than the
 total amount of recent data to examine.}
+
+\item{...}{Space to handle future expansions (unused).}
 }
 \value{
 A list containing updated parameter choices with class \code{arx_clist}.
diff --git a/man/arx_fcast_epi_workflow.Rd b/man/arx_fcast_epi_workflow.Rd
index 7a6b66305..8c76bcdd7 100644
--- a/man/arx_fcast_epi_workflow.Rd
+++ b/man/arx_fcast_epi_workflow.Rd
@@ -49,7 +49,7 @@ arx_fcast_epi_workflow(
 arx_fcast_epi_workflow(jhu, "death_rate",
   c("case_rate", "death_rate"),
   trainer = quantile_reg(),
-  args_list = arx_args_list(levels = 1:9 / 10)
+  args_list = arx_args_list(quantile_levels = 1:9 / 10)
 )
 }
 \seealso{
diff --git a/man/arx_forecaster.Rd b/man/arx_forecaster.Rd
index e121f272c..7a042c65c 100644
--- a/man/arx_forecaster.Rd
+++ b/man/arx_forecaster.Rd
@@ -49,7 +49,7 @@ out <- arx_forecaster(
 out <- arx_forecaster(jhu, "death_rate",
   c("case_rate", "death_rate"),
   trainer = quantile_reg(),
-  args_list = arx_args_list(levels = 1:9 / 10)
+  args_list = arx_args_list(quantile_levels = 1:9 / 10)
 )
 }
 \seealso{
diff --git a/man/dist_quantiles.Rd b/man/dist_quantiles.Rd
index 739bae5a8..57d2f3b3b 100644
--- a/man/dist_quantiles.Rd
+++ b/man/dist_quantiles.Rd
@@ -4,12 +4,12 @@
 \alias{dist_quantiles}
 \title{A distribution parameterized by a set of quantiles}
 \usage{
-dist_quantiles(x, tau)
+dist_quantiles(values, quantile_levels)
 }
 \arguments{
-\item{x}{A vector of values}
+\item{values}{A vector of values}
 
-\item{tau}{A vector of probabilities corresponding to \code{x}}
+\item{quantile_levels}{A vector of probabilities corresponding to \code{values}}
 }
 \description{
 A distribution parameterized by a set of quantiles
@@ -20,5 +20,8 @@ quantile(dstn, p = c(.1, .25, .5, .9))
 median(dstn)
 
 # it's a bit annoying to inspect the data
-vctrs::vec_data(vctrs::vec_data(dstn[1])[[1]])
+distributional::parameters(dstn[1])
+nested_quantiles(dstn[1])[[1]]
+
+dist_quantiles(1:4, 1:4 / 5)
 }
diff --git a/man/extract_argument.Rd b/man/extract_argument.Rd
index 3a83c4dd4..69c610c98 100644
--- a/man/extract_argument.Rd
+++ b/man/extract_argument.Rd
@@ -24,7 +24,7 @@ Extract an argument made to a frosting layer or recipe step
 \examples{
 f <- frosting() \%>\%
   layer_predict() \%>\%
-  layer_residual_quantiles(probs = c(0.0275, 0.975), symmetrize = FALSE) \%>\%
+  layer_residual_quantiles(quantile_levels = c(0.0275, 0.975), symmetrize = FALSE) \%>\%
   layer_naomit(.pred)
 
 extract_argument(f, "layer_residual_quantiles", "symmetrize")
diff --git a/man/extrapolate_quantiles.Rd b/man/extrapolate_quantiles.Rd
index cc6cb2c3c..619b2aa07 100644
--- a/man/extrapolate_quantiles.Rd
+++ b/man/extrapolate_quantiles.Rd
@@ -4,12 +4,12 @@
 \alias{extrapolate_quantiles}
 \title{Summarize a distribution with a set of quantiles}
 \usage{
-extrapolate_quantiles(x, p, ...)
+extrapolate_quantiles(x, probs, ...)
 }
 \arguments{
 \item{x}{a \code{distribution} vector}
 
-\item{p}{a vector of probabilities at which to calculate quantiles}
+\item{probs}{a vector of probabilities at which to calculate quantiles}
 
 \item{...}{additional arguments passed on to the \code{quantile} method}
 }
@@ -22,16 +22,16 @@ Summarize a distribution with a set of quantiles
 \examples{
 library(distributional)
 dstn <- dist_normal(c(10, 2), c(5, 10))
-extrapolate_quantiles(dstn, p = c(.25, 0.5, .75))
+extrapolate_quantiles(dstn, probs = c(.25, 0.5, .75))
 
 dstn <- dist_quantiles(list(1:4, 8:11), list(c(.2, .4, .6, .8)))
 # because this distribution is already quantiles, any extra quantiles are
 # appended
-extrapolate_quantiles(dstn, p = c(.25, 0.5, .75))
+extrapolate_quantiles(dstn, probs = c(.25, 0.5, .75))
 
 dstn <- c(
   dist_normal(c(10, 2), c(5, 10)),
   dist_quantiles(list(1:4, 8:11), list(c(.2, .4, .6, .8)))
 )
-extrapolate_quantiles(dstn, p = c(.25, 0.5, .75))
+extrapolate_quantiles(dstn, probs = c(.25, 0.5, .75))
 }
diff --git a/man/flatline_args_list.Rd b/man/flatline_args_list.Rd
index dcae448f1..059dfa038 100644
--- a/man/flatline_args_list.Rd
+++ b/man/flatline_args_list.Rd
@@ -9,11 +9,12 @@ flatline_args_list(
   n_training = Inf,
   forecast_date = NULL,
   target_date = NULL,
-  levels = c(0.05, 0.95),
+  quantile_levels = c(0.05, 0.95),
   symmetrize = TRUE,
   nonneg = TRUE,
   quantile_by_key = character(0L),
-  nafill_buffer = Inf
+  nafill_buffer = Inf,
+  ...
 )
 }
 \arguments{
@@ -34,7 +35,7 @@ The default \code{NULL} will attempt to determine this automatically.}
 \item{target_date}{Date. The date for which the forecast is intended.
 The default \code{NULL} will attempt to determine this automatically.}
 
-\item{levels}{Vector or \code{NULL}. A vector of probabilities to produce
+\item{quantile_levels}{Vector or \code{NULL}. A vector of probabilities to produce
 prediction intervals. These are created by computing the quantiles of
 training residuals. A \code{NULL} value will result in point forecasts only.}
 
@@ -62,6 +63,8 @@ we require at least \code{min(lags)} rows of recent data per \code{geo_value} to
 create a prediction. For this reason, setting \code{nafill_buffer < min(lags)}
 will be treated as \emph{additional} allowed recent data rather than the
 total amount of recent data to examine.}
+
+\item{...}{Space to handle future expansions (unused).}
 }
 \value{
 A list containing updated parameter choices with class \code{flatline_alist}.
@@ -72,5 +75,5 @@ Constructs a list of arguments for \code{\link[=flatline_forecaster]{flatline_fo
 \examples{
 flatline_args_list()
 flatline_args_list(symmetrize = FALSE)
-flatline_args_list(levels = c(.1, .3, .7, .9), n_training = 120)
+flatline_args_list(quantile_levels = c(.1, .3, .7, .9), n_training = 120)
 }
diff --git a/man/layer_point_from_distn.Rd b/man/layer_point_from_distn.Rd
index cc2dcf2fe..7ad69a332 100644
--- a/man/layer_point_from_distn.Rd
+++ b/man/layer_point_from_distn.Rd
@@ -42,7 +42,7 @@ r <- epi_recipe(jhu) \%>\%
   step_epi_ahead(death_rate, ahead = 7) \%>\%
   step_epi_naomit()
 
-wf <- epi_workflow(r, quantile_reg(tau = c(.25, .5, .75))) \%>\% fit(jhu)
+wf <- epi_workflow(r, quantile_reg(quantile_levels = c(.25, .5, .75))) \%>\% fit(jhu)
 
 latest <- get_test_data(recipe = r, x = jhu)
 
diff --git a/man/layer_quantile_distn.Rd b/man/layer_quantile_distn.Rd
index f53e78356..167282760 100644
--- a/man/layer_quantile_distn.Rd
+++ b/man/layer_quantile_distn.Rd
@@ -7,7 +7,7 @@
 layer_quantile_distn(
   frosting,
   ...,
-  levels = c(0.25, 0.75),
+  quantile_levels = c(0.25, 0.75),
   truncate = c(-Inf, Inf),
   name = ".pred_distn",
   id = rand_id("quantile_distn")
@@ -18,7 +18,7 @@ layer_quantile_distn(
 
 \item{...}{Unused, include for consistency with other layers.}
 
-\item{levels}{a vector of probabilities (quantiles) to extract}
+\item{quantile_levels}{a vector of probabilities to extract}
 
 \item{truncate}{Do we truncate the distribution to an interval}
 
@@ -34,7 +34,7 @@ quantiles will be added to the predictions.
 This function calculates quantiles when the prediction was \emph{distributional}.
 Currently, the only distributional engine is \code{quantile_reg()}.
 If this engine is used, then this layer will grab out estimated (or extrapolated)
-quantiles at the requested levels.
+quantiles at the requested quantile values.
 }
 \examples{
 jhu <- case_death_rate_subset \%>\%
@@ -45,7 +45,8 @@ r <- epi_recipe(jhu) \%>\%
   step_epi_ahead(death_rate, ahead = 7) \%>\%
   step_epi_naomit()
 
-wf <- epi_workflow(r, quantile_reg(tau = c(.25, .5, .75))) \%>\% fit(jhu)
+wf <- epi_workflow(r, quantile_reg(quantile_levels = c(.25, .5, .75))) \%>\%
+  fit(jhu)
 
 latest <- get_test_data(recipe = r, x = jhu)
 
diff --git a/man/layer_residual_quantiles.Rd b/man/layer_residual_quantiles.Rd
index 412dbc86e..bf0e05be1 100644
--- a/man/layer_residual_quantiles.Rd
+++ b/man/layer_residual_quantiles.Rd
@@ -7,7 +7,7 @@
 layer_residual_quantiles(
   frosting,
   ...,
-  probs = c(0.05, 0.95),
+  quantile_levels = c(0.05, 0.95),
   symmetrize = TRUE,
   by_key = character(0L),
   name = ".pred_distn",
@@ -19,7 +19,7 @@ layer_residual_quantiles(
 
 \item{...}{Unused, include for consistency with other layers.}
 
-\item{probs}{numeric vector of probabilities with values in (0,1)
+\item{quantile_levels}{numeric vector of probabilities with values in (0,1)
 referring to the desired quantile.}
 
 \item{symmetrize}{logical. If \code{TRUE} then interval will be symmetric.}
@@ -53,7 +53,7 @@ latest <- get_test_data(recipe = r, x = jhu)
 
 f <- frosting() \%>\%
   layer_predict() \%>\%
-  layer_residual_quantiles(probs = c(0.0275, 0.975), symmetrize = FALSE) \%>\%
+  layer_residual_quantiles(quantile_levels = c(0.0275, 0.975), symmetrize = FALSE) \%>\%
   layer_naomit(.pred)
 wf1 <- wf \%>\% add_frosting(f)
 
@@ -61,7 +61,7 @@ p <- predict(wf1, latest)
 
 f2 <- frosting() \%>\%
   layer_predict() \%>\%
-  layer_residual_quantiles(probs = c(0.3, 0.7), by_key = "geo_value") \%>\%
+  layer_residual_quantiles(quantile_levels = c(0.3, 0.7), by_key = "geo_value") \%>\%
   layer_naomit(.pred)
 wf2 <- wf \%>\% add_frosting(f2)
 
diff --git a/man/quantile_reg.Rd b/man/quantile_reg.Rd
index dce711455..8e576ac84 100644
--- a/man/quantile_reg.Rd
+++ b/man/quantile_reg.Rd
@@ -4,7 +4,7 @@
 \alias{quantile_reg}
 \title{Quantile regression}
 \usage{
-quantile_reg(mode = "regression", engine = "rq", tau = 0.5)
+quantile_reg(mode = "regression", engine = "rq", quantile_levels = 0.5)
 }
 \arguments{
 \item{mode}{A single character string for the type of model.
@@ -13,7 +13,7 @@ The only possible value for this model is "regression".}
 \item{engine}{Character string naming the fitting function. Currently, only
 "rq" is supported.}
 
-\item{tau}{A scalar or vector of values in (0, 1) to determine which
+\item{quantile_levels}{A scalar or vector of values in (0, 1) to determine which
 quantiles to estimate (default is 0.5).}
 }
 \description{
@@ -23,7 +23,7 @@ only supported engine is "rq" which uses \code{\link[quantreg:rq]{quantreg::rq()
 }
 \examples{
 tib <- data.frame(y = rnorm(100), x1 = rnorm(100), x2 = rnorm(100))
-rq_spec <- quantile_reg(tau = c(.2, .8)) \%>\% set_engine("rq")
+rq_spec <- quantile_reg(quantile_levels = c(.2, .8)) \%>\% set_engine("rq")
 ff <- rq_spec \%>\% fit(y ~ ., data = tib)
 predict(ff, new_data = tib)
 }
diff --git a/man/smooth_quantile_reg.Rd b/man/smooth_quantile_reg.Rd
index b938541f1..42a951759 100644
--- a/man/smooth_quantile_reg.Rd
+++ b/man/smooth_quantile_reg.Rd
@@ -8,7 +8,7 @@ smooth_quantile_reg(
   mode = "regression",
   engine = "smoothqr",
   outcome_locations = NULL,
-  tau = 0.5,
+  quantile_levels = 0.5,
   degree = 3L
 )
 }
@@ -24,7 +24,7 @@ responses are observed at a different spacing (or appear in a different
 order), that information should be used here. This
 argument will be mapped to the \code{ahead} argument of \code{\link[smoothqr:smooth_qr]{smoothqr::smooth_qr()}}.}
 
-\item{tau}{A scalar or vector of values in (0, 1) to determine which
+\item{quantile_levels}{A scalar or vector of values in (0, 1) to determine which
 quantiles to estimate (default is 0.5).}
 
 \item{degree}{the number of polynomials used for response smoothing. Must
@@ -41,7 +41,7 @@ tib <- data.frame(
   y4 = rnorm(100), y5 = rnorm(100), y6 = rnorm(100),
   x1 = rnorm(100), x2 = rnorm(100)
 )
-qr_spec <- smooth_quantile_reg(tau = c(.2, .5, .8), outcome_locations = 1:6)
+qr_spec <- smooth_quantile_reg(quantile_levels = c(.2, .5, .8), outcome_locations = 1:6)
 ff <- qr_spec \%>\% fit(cbind(y1, y2, y3, y4, y5, y6) ~ ., data = tib)
 p <- predict(ff, new_data = tib)
 
@@ -50,7 +50,7 @@ y <- sin(x) + rnorm(length(x), sd = .1)
 fd <- x[length(x) - 20]
 XY <- smoothqr::lagmat(y[1:(length(y) - 20)], c(-20:20))
 XY <- tibble::as_tibble(XY)
-qr_spec <- smooth_quantile_reg(tau = c(.2, .5, .8), outcome_locations = 20:1)
+qr_spec <- smooth_quantile_reg(quantile_levels = c(.2, .5, .8), outcome_locations = 20:1)
 tt <- qr_spec \%>\% fit_xy(x = XY[, 21:41], y = XY[, 1:20])
 
 library(tidyr)
@@ -67,7 +67,7 @@ pl <- pl \%>\%
     x = x[length(x) - 20] + ahead / 100 * 2 * pi,
     ahead = NULL
   ) \%>\%
-  pivot_wider(names_from = tau, values_from = q)
+  pivot_wider(names_from = quantile_levels, values_from = values)
 plot(x, y, pch = 16, xlim = c(pi, 2 * pi), col = "lightgrey")
 curve(sin(x), add = TRUE)
 abline(v = fd, lty = 2)
diff --git a/tests/testthat/test-arx_args_list.R b/tests/testthat/test-arx_args_list.R
index dcd7a1cfe..138a75e87 100644
--- a/tests/testthat/test-arx_args_list.R
+++ b/tests/testthat/test-arx_args_list.R
@@ -13,9 +13,9 @@ test_that("arx_args checks inputs", {
   expect_error(arx_args_list(symmetrize = 4))
   expect_error(arx_args_list(nonneg = 4))
 
-  expect_error(arx_args_list(levels = -.1))
-  expect_error(arx_args_list(levels = 1.1))
-  expect_type(arx_args_list(levels = NULL), "list")
+  expect_error(arx_args_list(quantile_levels = -.1))
+  expect_error(arx_args_list(quantile_levels = 1.1))
+  expect_type(arx_args_list(quantile_levels = NULL), "list")
 
   expect_error(arx_args_list(target_date = "2022-01-01"))
   expect_identical(
@@ -25,8 +25,8 @@ test_that("arx_args checks inputs", {
 })
 
 test_that("arx forecaster disambiguates quantiles", {
-  alist <- eval(formals(arx_args_list)$levels)
-  tlist <- eval(formals(quantile_reg)$tau)
+  alist <- eval(formals(arx_args_list)$quantile_levels)
+  tlist <- eval(formals(quantile_reg)$quantile_levels)
   expect_identical( # both default
     compare_quantile_args(alist, tlist),
     sort(c(alist, tlist))
@@ -36,7 +36,7 @@ test_that("arx forecaster disambiguates quantiles", {
     compare_quantile_args(alist, tlist),
     sort(unique(alist))
   )
-  alist <- eval(formals(arx_args_list)$levels)
+  alist <- eval(formals(arx_args_list)$quantile_levels)
   tlist <- c(.05, .95, tlist)
   expect_identical( # alist is default, should give tlist
     compare_quantile_args(alist, tlist),
diff --git a/tests/testthat/test-dist_quantiles.R b/tests/testthat/test-dist_quantiles.R
index 07d1530d2..4fc5587d4 100644
--- a/tests/testthat/test-dist_quantiles.R
+++ b/tests/testthat/test-dist_quantiles.R
@@ -20,7 +20,7 @@ test_that("tail functions give reasonable output", {
 })
 
 test_that("single dist_quantiles works, quantiles are accessible", {
-  z <- new_quantiles(q = 1:5, tau = c(.2, .4, .5, .6, .8))
+  z <- new_quantiles(values = 1:5, quantile_levels = c(.2, .4, .5, .6, .8))
   expect_s3_class(z, "dist_quantiles")
   expect_equal(median(z), 3)
   expect_equal(quantile(z, c(.2, .4, .5, .6, .8)), 1:5)
@@ -30,7 +30,7 @@ test_that("single dist_quantiles works, quantiles are accessible", {
   expect_equal(quantile(z, c(.3, .7), middle = "cubic"), Q(c(.3, .7)))
   expect_identical(
     extrapolate_quantiles(z, c(.3, .7), middle = "linear"),
-    new_quantiles(q = c(1, 1.5, 2, 3, 4, 4.5, 5), tau = 2:8 / 10)
+    new_quantiles(values = c(1, 1.5, 2, 3, 4, 4.5, 5), quantile_levels = 2:8 / 10)
   )
 })
 
diff --git a/tests/testthat/test-extract_argument.R b/tests/testthat/test-extract_argument.R
index 974a50888..0654304ba 100644
--- a/tests/testthat/test-extract_argument.R
+++ b/tests/testthat/test-extract_argument.R
@@ -1,26 +1,26 @@
 test_that("layer argument extractor works", {
   f <- frosting() %>%
     layer_predict() %>%
-    layer_residual_quantiles(probs = c(0.0275, 0.975), symmetrize = FALSE) %>%
+    layer_residual_quantiles(quantile_levels = c(0.0275, 0.975), symmetrize = FALSE) %>%
     layer_naomit(.pred)
 
   expect_error(extract_argument(f$layers[[1]], "uhoh", "bubble"))
   expect_error(extract_argument(f$layers[[1]], "layer_predict", "bubble"))
   expect_identical(
-    extract_argument(f$layers[[2]], "layer_residual_quantiles", "probs"),
+    extract_argument(f$layers[[2]], "layer_residual_quantiles", "quantile_levels"),
     c(0.0275, 0.9750)
   )
 
-  expect_error(extract_argument(f, "layer_thresh", "probs"))
+  expect_error(extract_argument(f, "layer_thresh", "quantile_levels"))
   expect_identical(
-    extract_argument(f, "layer_residual_quantiles", "probs"),
+    extract_argument(f, "layer_residual_quantiles", "quantile_levels"),
     c(0.0275, 0.9750)
   )
 
   wf <- epi_workflow(postprocessor = f)
-  expect_error(extract_argument(epi_workflow(), "layer_residual_quantiles", "probs"))
+  expect_error(extract_argument(epi_workflow(), "layer_residual_quantiles", "quantile_levels"))
   expect_identical(
-    extract_argument(wf, "layer_residual_quantiles", "probs"),
+    extract_argument(wf, "layer_residual_quantiles", "quantile_levels"),
     c(0.0275, 0.9750)
   )
 
@@ -46,7 +46,7 @@ test_that("recipe argument extractor works", {
   expect_identical(extract_argument(r$steps[[2]], "step_epi_ahead", "ahead"), 7)
 
 
-  expect_error(extract_argument(r, "step_lightly", "probs"))
+  expect_error(extract_argument(r, "step_lightly", "quantile_levels"))
   expect_identical(
     extract_argument(r, "step_epi_lag", "lag"),
     list(c(0, 7, 14), c(0, 7, 14))
diff --git a/tests/testthat/test-layer_residual_quantiles.R b/tests/testthat/test-layer_residual_quantiles.R
index 967eee1a5..a2c7bad4e 100644
--- a/tests/testthat/test-layer_residual_quantiles.R
+++ b/tests/testthat/test-layer_residual_quantiles.R
@@ -14,7 +14,7 @@ test_that("Returns expected number or rows and columns", {
   f <- frosting() %>%
     layer_predict() %>%
     layer_naomit(.pred) %>%
-    layer_residual_quantiles(probs = c(0.0275, 0.8, 0.95), symmetrize = FALSE)
+    layer_residual_quantiles(quantile_levels = c(0.0275, 0.8, 0.95), symmetrize = FALSE)
 
   wf1 <- wf %>% add_frosting(f)
 
@@ -28,5 +28,5 @@ test_that("Returns expected number or rows and columns", {
   unnested <- nested %>% tidyr::unnest(.quantiles)
 
   expect_equal(nrow(unnested), 9L)
-  expect_equal(unique(unnested$tau), c(.0275, .8, .95))
+  expect_equal(unique(unnested$quantile_levels), c(.0275, .8, .95))
 })
diff --git a/tests/testthat/test-layer_threshold_preds.R b/tests/testthat/test-layer_threshold_preds.R
index 80b6a42a9..9df7e64ab 100644
--- a/tests/testthat/test-layer_threshold_preds.R
+++ b/tests/testthat/test-layer_threshold_preds.R
@@ -43,7 +43,7 @@ test_that("Specified pred_lower and pred_upper work as intended", {
 test_that("thresholds additional columns", {
   f <- frosting() %>%
     layer_predict() %>%
-    layer_residual_quantiles(probs = c(.1, .9)) %>%
+    layer_residual_quantiles(quantile_levels = c(.1, .9)) %>%
     layer_threshold(.pred, .pred_distn, lower = 0.180, upper = 0.31) %>%
     layer_naomit(.pred)
 
@@ -58,6 +58,6 @@ test_that("thresholds additional columns", {
   p <- p %>%
     dplyr::mutate(.quantiles = nested_quantiles(.pred_distn)) %>%
     tidyr::unnest(.quantiles)
-  expect_equal(round(p$q, digits = 3), c(0.180, 0.31, 0.180, .18, 0.310, .31))
-  expect_equal(p$tau, rep(c(.1, .9), times = 3))
+  expect_equal(round(p$values, digits = 3), c(0.180, 0.31, 0.180, .18, 0.310, .31))
+  expect_equal(p$quantile_levels, rep(c(.1, .9), times = 3))
 })
diff --git a/tests/testthat/test-pivot_quantiles.R b/tests/testthat/test-pivot_quantiles.R
index 9928c5e09..908a75795 100644
--- a/tests/testthat/test-pivot_quantiles.R
+++ b/tests/testthat/test-pivot_quantiles.R
@@ -38,7 +38,7 @@ test_that("quantile pivotting longer behaves", {
   tib$d1 <- d1
   expect_length(pivot_quantiles_longer(tib, d1), 5L)
   expect_identical(nrow(pivot_quantiles_longer(tib, d1)), 7L)
-  expect_identical(pivot_quantiles_longer(tib, d1)$q, as.double(c(1:3, 2:5)))
+  expect_identical(pivot_quantiles_longer(tib, d1)$values, as.double(c(1:3, 2:5)))
 
   d1 <- c(dist_quantiles(1:3, 1:3 / 4), dist_quantiles(2:4, 2:4 / 4))
   tib$d1 <- d1
@@ -62,7 +62,7 @@ test_that("quantile pivotting longer behaves", {
     6L
   )
   expect_identical(
-    pivot_quantiles_longer(tib, d1, d3, .ignore_length_check = TRUE)$d1_q,
+    pivot_quantiles_longer(tib, d1, d3, .ignore_length_check = TRUE)$d1_values,
     as.double(rep(c(1:3, 2:4), each = 4))
   )
 })
diff --git a/vignettes/epipredict.Rmd b/vignettes/epipredict.Rmd
index 17a604504..8a8e20e10 100644
--- a/vignettes/epipredict.Rmd
+++ b/vignettes/epipredict.Rmd
@@ -157,7 +157,7 @@ Another property of the basic model is the predictive interval. We describe this
 ```{r differential-levels}
 out_q <- arx_forecaster(jhu, "death_rate", c("case_rate", "death_rate"),
   args_list = arx_args_list(
-    levels = c(.01, .025, seq(.05, .95, by = .05), .975, .99)
+    quantile_levels = c(.01, .025, seq(.05, .95, by = .05), .975, .99)
   )
 )
 ```
@@ -168,7 +168,8 @@ The column `.pred_dstn` in the `predictions` object is actually a "distribution"
 head(quantile(out_q$predictions$.pred_distn, p = .4))
 ```
 
-or extract the entire distribution into a "long" `epi_df` with `tau` being the probability and `q` being the value associated to that quantile.
+or extract the entire distribution into a "long" `epi_df` with `quantile_levels`
+being the probability and `values` being the value associated to that quantile.
 
 ```{r q2}
 out_q$predictions %>%
@@ -182,7 +183,7 @@ Additional simple adjustments to the basic forecaster can be made using the func
 ```{r, eval = FALSE}
 arx_args_list(
   lags = c(0L, 7L, 14L), ahead = 7L, n_training = Inf,
-  forecast_date = NULL, target_date = NULL, levels = c(0.05, 0.95),
+  forecast_date = NULL, target_date = NULL, quantile_levels = c(0.05, 0.95),
   symmetrize = TRUE, nonneg = TRUE, quantile_by_key = character(0L),
   nafill_buffer = Inf
 )
@@ -213,7 +214,7 @@ out_gb <- arx_forecaster(
 Or quantile regression, using our custom forecasting engine `quantile_reg()`:
 
 ```{r quantreg, warning = FALSE}
-out_gb <- arx_forecaster(
+out_qr <- arx_forecaster(
   jhu, "death_rate", c("case_rate", "death_rate"),
   quantile_reg()
 )
@@ -340,7 +341,7 @@ intervals at 0. The code to do this (inside the forecaster) is
 f <- frosting() %>%
   layer_predict() %>%
   layer_residual_quantiles(
-    probs = c(.01, .025, seq(.05, .95, by = .05), .975, .99),
+    quantile_levels = c(.01, .025, seq(.05, .95, by = .05), .975, .99),
     symmetrize = TRUE
   ) %>%
   layer_add_forecast_date() %>%
diff --git a/vignettes/preprocessing-and-models.Rmd b/vignettes/preprocessing-and-models.Rmd
index f85f35f71..60291ffd1 100644
--- a/vignettes/preprocessing-and-models.Rmd
+++ b/vignettes/preprocessing-and-models.Rmd
@@ -354,7 +354,7 @@ f <- frosting() %>%
     df_pop_col = "pop"
   )
 
-wf <- epi_workflow(r, quantile_reg(tau = c(.05, .5, .95))) %>%
+wf <- epi_workflow(r, quantile_reg(quantile_levels = c(.05, .5, .95))) %>%
   fit(jhu) %>%
   add_frosting(f)
 
@@ -373,7 +373,7 @@ p %>%
   select(geo_value, target_date, .pred_scaled, .pred_distn_scaled) %>%
   mutate(.pred_distn_scaled = nested_quantiles(.pred_distn_scaled)) %>%
   unnest(.pred_distn_scaled) %>%
-  pivot_wider(names_from = tau, values_from = q)
+  pivot_wider(names_from = quantile_levels, values_from = values)
 ```
 
 Last but not least, let's take a look at the regression fit and check the