From b172b67dcc538dcd0db737f86cc58845e97807a2 Mon Sep 17 00:00:00 2001
From: "Daniel J. McDonald" <dajmcdon@gmail.com>
Date: Mon, 7 Jul 2025 09:15:02 -0700
Subject: [PATCH 1/8] cp epipredict impl

---
 R/quant-weighted_interval.R | 109 ++++++++++++++++++++++++++++++++++++
 1 file changed, 109 insertions(+)
 create mode 100644 R/quant-weighted_interval.R

diff --git a/R/quant-weighted_interval.R b/R/quant-weighted_interval.R
new file mode 100644
index 00000000..47187fe5
--- /dev/null
+++ b/R/quant-weighted_interval.R
@@ -0,0 +1,109 @@
+#' Compute weighted interval score
+#'
+#' Weighted interval score (WIS), a well-known quantile-based
+#' approximation of the commonly-used continuous ranked probability score
+#' (CRPS). WIS is a proper score, and can be thought of as a distributional
+#' generalization of absolute error. For example, see [Bracher et
+#' al. (2020)](https://arxiv.org/abs/2005.12881) for discussion in the context
+#' of COVID-19 forecasting.
+#'
+#' @param x A vector of class `quantile_pred`.
+#' @param actual double. Actual value(s)
+#' @param quantile_levels probabilities. If specified, the score will be
+#'   computed at this set of levels. Otherwise, those present in `x` will be
+#'   used.
+#' @param na_handling character. Determines missing values are handled.
+#'   For `"impute"`, missing values will be
+#'   calculated if possible using the available quantiles. For `"drop"`,
+#'   explicitly missing values are ignored in the calculation of the score, but
+#'   implicitly missing values are imputed if possible.
+#'   For `"propogate"`, the resulting score will be `NA` if any missing values
+#'   exist. Finally, if
+#'   `quantile_levels` is specified, `"fail"` will result in
+#'   the score being `NA` when any required quantile levels (implicit or explicit)
+#'   do not have corresponding values.
+#' @param ... not used
+#'
+#' @return a vector of nonnegative scores.
+#'
+#' @export
+#' @examples
+#' quantile_levels <- c(.2, .4, .6, .8)
+#' predq1 <- 1:4 #
+#' predq2 <- 8:11
+#' dstn <- quantile_pred(rbind(predq1, predq2), quantile_levels)
+#' actual <- c(3.3, 7.1)
+#' weighted_interval_score(dstn, actual)
+#' weighted_interval_score(dstn, actual, c(.25, .5, .75))
+#'
+#' # Missing value behaviours
+#' dstn <- quantile_pred(matrix(c(1, 2, NA, 4), nrow = 1), 1:4 / 5)
+#' weighted_interval_score(dstn, 2.5)
+#' weighted_interval_score(dstn, 2.5, 1:9 / 10)
+#' weighted_interval_score(dstn, 2.5, 1:9 / 10, na_handling = "drop")
+#' weighted_interval_score(dstn, 2.5, na_handling = "propagate")
+#' weighted_interval_score(
+#'   quantile_pred(matrix(1:4, nrow = 1), 1:4 / 5),
+#'   actual = 2.5,
+#'   quantile_levels = 1:9 / 10,
+#'   na_handling = "fail"
+#' )
+#'
+#'
+#' # Using some actual forecasts --------
+#' library(dplyr)
+#' training <- covid_case_death_rates %>%
+#'   filter(time_value >= "2021-10-01", time_value <= "2021-12-01")
+#' preds <- flatline_forecaster(
+#'   training, "death_rate",
+#'   flatline_args_list(quantile_levels = c(.01, .025, 1:19 / 20, .975, .99))
+#' )$predictions
+#' actuals <- covid_case_death_rates %>%
+#'   filter(time_value == as.Date("2021-12-01") + 7) %>%
+#'   select(geo_value, time_value, actual = death_rate)
+#' preds <- left_join(preds, actuals,
+#'   by = c("target_date" = "time_value", "geo_value")
+#' ) %>%
+#'   mutate(wis = weighted_interval_score(.pred_distn, actual))
+#' preds
+weighted_interval_score <- function(
+    x,
+    actual,
+    quantile_levels = NULL,
+    na_handling = c("impute", "drop", "propagate", "fail"),
+    ...) {
+  UseMethod("weighted_interval_score")
+}
+
+
+#' @export
+weighted_interval_score.quantile_pred <- function(
+    x, actual,
+    quantile_levels = NULL,
+    na_handling = c("impute", "drop", "propagate", "fail"),
+    ...) {
+  rlang::check_dots_empty()
+  n <- vctrs::vec_size(x)
+  if (length(actual) == 1L) actual <- rep(actual, n)
+  assert_numeric(actual, finite = TRUE, len = n)
+  assert_numeric(quantile_levels, lower = 0, upper = 1, null.ok = TRUE)
+  na_handling <- rlang::arg_match(na_handling)
+  old_quantile_levels <- x %@% "quantile_levels"
+  if (na_handling == "fail") {
+    if (is.null(quantile_levels)) {
+      cli_abort('`na_handling = "fail"` requires `quantile_levels` to be specified.')
+    }
+    if (!all(quantile_levels %in% old_quantile_levels)) {
+      return(rep(NA_real_, n))
+    }
+  }
+  tau <- quantile_levels %||% old_quantile_levels
+  x <- extrapolate_quantiles(x, tau, replace_na = (na_handling == "impute"))
+  x <- as.matrix(x)[, attr(x, "quantile_levels") %in% tau, drop = FALSE]
+  na_rm <- (na_handling == "drop")
+  map2_dbl(vctrs::vec_chop(x), actual, ~ wis_one_quantile(.x, tau, .y, na_rm))
+}
+
+wis_one_quantile <- function(q, tau, actual, na_rm) {
+  2 * mean(pmax(tau * (actual - q), (1 - tau) * (q - actual)), na.rm = na_rm)
+}

From da8d2301e4ad5437f0bdf9adceb8473ab7cccd58 Mon Sep 17 00:00:00 2001
From: "Daniel J. McDonald" <dajmcdon@gmail.com>
Date: Mon, 7 Jul 2025 10:15:31 -0700
Subject: [PATCH 2/8] add quantile_pred validator

---
 R/validation.R | 40 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 40 insertions(+)

diff --git a/R/validation.R b/R/validation.R
index 547e072a..d75be0de 100644
--- a/R/validation.R
+++ b/R/validation.R
@@ -443,3 +443,43 @@ validate_case_weights <- function(case_weights, size, call = caller_env()) {
 
   invisible(NULL)
 }
+
+validate_numeric_truth_quantile_estimate <- function(
+    truth,
+    estimate,
+    call = caller_env()
+) {
+  if (!is.numeric(truth)) {
+    cli::cli_abort(
+      "{.arg truth} should be a numeric vector,
+      not {.obj_type_friendly {truth}}.",
+      call = call
+    )
+  }
+
+  if (!inherits(estimate, "quantile_pred")) {
+    cli::cli_abort(
+      "{.arg estimate} should be a {.cls quantile_pred} object,
+      not {.obj_type_friendly {estimate}}.",
+      call = call
+    )
+  }
+
+  if (is.matrix(truth)) {
+    cli::cli_abort(
+      "{.arg truth} should be a numeric vector, not a numeric matrix.",
+      call = call
+    )
+  }
+
+  n_truth <- length(truth)
+  n_estimate <- vctrs::vec_size(estimate)
+
+  if (n_truth != n_estimate) {
+    cli::cli_abort(
+      "{.arg truth} ({n_truth}) and
+      {.arg estimate} ({n_estimate}) must be the same length.",
+      call = call
+    )
+  }
+}

From 8150d3974d616627c2b8c43c3a5cece19a93df3e Mon Sep 17 00:00:00 2001
From: "Daniel J. McDonald" <dajmcdon@gmail.com>
Date: Mon, 7 Jul 2025 10:57:21 -0700
Subject: [PATCH 3/8] add a check for quantile metrics

---
 R/check-metric.R | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/R/check-metric.R b/R/check-metric.R
index eebdfb46..8b2042a4 100644
--- a/R/check-metric.R
+++ b/R/check-metric.R
@@ -15,6 +15,7 @@
 #'   - For `check_ordered_prob_metric()`, an ordered factor.
 #'   - For `check_dynamic_survival_metric()`, a Surv object.
 #'   - For `check_static_survival_metric()`, a Surv object.
+#'   - For `check_quantile_metric()`, a numeric vector.
 #'
 #' @param estimate The realized `estimate` result.
 #'   - For `check_numeric_metric()`, a numeric vector.
@@ -25,6 +26,7 @@
 #'     a numeric matrix for multic-class `truth`.
 #'   - For `check_dynamic_survival_metric()`, list-column of data.frames.
 #'   - For `check_static_survival_metric()`, a numeric vector.
+#'   - For `check_quantile_metric()`, a `hardhat::quantile_pred` vector.
 #'
 #' @param case_weights The realized case weights, as a numeric vector. This must
 #'   be the same length as `truth`.
@@ -120,3 +122,15 @@ check_static_survival_metric <- function(
   validate_case_weights(case_weights, size = nrow(truth), call = call)
   validate_surv_truth_numeric_estimate(truth, estimate, call = call)
 }
+
+#' @rdname check_metric
+#' @export
+check_quantile_metric <- function(
+  truth,
+  estimate,
+  case_weights,
+  call = caller_env()
+) {
+  validate_numeric_truth_quantile_estimate(truth, estimate, call = call)
+  validate_case_weights(case_weights, size = nrow(truth), call = call)
+}

From 3cb06ea4476771d25abca730f699903259ccfedc Mon Sep 17 00:00:00 2001
From: "Daniel J. McDonald" <dajmcdon@gmail.com>
Date: Mon, 7 Jul 2025 10:58:03 -0700
Subject: [PATCH 4/8] internal functions for interpolating/extrapolating the
 quantile_pred to additional quantiles

---
 R/interpolate-quantile_pred.R | 149 ++++++++++++++++++++++++++++++++++
 1 file changed, 149 insertions(+)
 create mode 100644 R/interpolate-quantile_pred.R

diff --git a/R/interpolate-quantile_pred.R b/R/interpolate-quantile_pred.R
new file mode 100644
index 00000000..e44f38b3
--- /dev/null
+++ b/R/interpolate-quantile_pred.R
@@ -0,0 +1,149 @@
+# Move these to hardhat? --------------------------------------------------
+
+#' restrict various objects to the interval \[lower, upper\]
+#' @param x the object to restrict
+#' @param lower numeric, the lower bound
+#' @param upper numeric, the upper bound
+#' @param ... unused
+#' @export
+#' @keywords internal
+snap <- function(x, lower, upper, ...) {
+  UseMethod("snap")
+}
+
+#' @export
+snap.default <- function(x, lower, upper, ...) {
+  rlang::check_dots_empty()
+  check_number_decimal(lower)
+  check_number_decimal(upper)
+  if (!is.numeric(x)) {
+    cli::cli_abort(
+      "{.arg x} should be should be numeric,
+      not {.obj_type_friendly {x}}."
+    )
+  }
+  pmin(pmax(x, lower), upper)
+}
+
+#' @export
+snap.quantile_pred <- function(x, lower, upper, ...) {
+  values <- as.matrix(x)
+  quantile_levels <- x %@% "quantile_levels"
+  values <- lapply(vctrs::vec_chop(values), function(.x) snap(.x, lower, upper))
+  hardhat::quantile_pred(do.call(rbind, values), quantile_levels = quantile_levels)
+}
+
+
+impute_quantiles <- function(
+    x, probs = seq(0, 1, 0.25), na.rm = FALSE, lower = -Inf,
+    upper = Inf, middle = c("cubic", "linear"), ...)
+{
+  check_bool(na.rm)
+  if (is.unsorted(probs)) probs <- sort(probs)
+  hardhat::check_quantile_levels(probs)
+  check_number_decimal(lower, allow_na = na.rm)
+  check_number_decimal(upper, allow_na = na.rm)
+  if (lower > upper) {
+    cli::cli_abort("`lower` must be less than `upper`.")
+  }
+  middle <- rlang::arg_match(middle)
+  snap(impute_quantile_internal(x, probs, middle), lower, upper)
+}
+
+impute_quantile_internal <- function(x, tau_out, middle) {
+  tau <- x %@% "quantile_levels"
+  qvals <- as.matrix(x)
+  if (all(tau_out %in% tau) && !anyNA(qvals)) {
+    return(qvals[, match(tau_out, tau), drop = FALSE])
+  }
+  if (length(tau) < 2) {
+    cli::cli_abort(
+      "Quantile interpolation is not possible when fewer than 2 quantiles
+      are avaliable."
+    )
+  }
+  qvals_out <- lapply(
+    vctrs::vec_chop(qvals),
+    function(.x) impute_quantiles_single(.x, tau, tau_out, middle)
+  )
+  qvals_out <- do.call(rbind, qvals_out)
+  qvals_out
+}
+
+impute_quantiles_single <- function(qvals, tau, tau_out, middle) {
+  qvals_out <- rep(NA, length(tau_out))
+  good <- !is.na(qvals)
+  if (!any(good)) {
+    return(qvals_out)
+  }
+  qvals <- qvals[good]
+  tau <- tau[good]
+
+  # in case we only have one point, and it matches something we wanted
+  if (length(good) < 2) {
+    matched_one <- tau_out %in% tau
+    qvals_out[matched_one] <- qvals[matched_one]
+    return(qvals_out)
+  }
+
+  indl <- tau_out < min(tau)
+  indr <- tau_out > max(tau)
+  indm <- !indl & !indr
+
+  if (middle == "cubic") {
+    method <- "cubic"
+    result <- tryCatch(
+      {
+        Q <- stats::splinefun(tau, qvals, method = "hyman")
+        quartiles <- Q(c(.25, .5, .75))
+      },
+      error = function(e) {
+        return(NA)
+      }
+    )
+  }
+  if (middle == "linear" || any(is.na(result))) {
+    method <- "linear"
+    quartiles <- stats::approx(tau, qvals, c(.25, .5, .75))$y
+  }
+  if (any(indm)) {
+    qvals_out[indm] <- switch(
+      method,
+      linear = stats::approx(tau, qvals, tau_out[indm])$y,
+      cubic = Q(tau_out[indm])
+    )
+  }
+  if (any(indl) || any(indr)) {
+    qv <- data.frame(
+      q = c(tau, tau_out[indm]),
+      v = c(qvals, qvals_out[indm])
+    ) %>%
+      dplyr::distinct(q, .keep_all = TRUE) %>%
+      dplyr::arrange(q)
+  }
+  if (any(indl)) {
+    qvals_out[indl] <- tail_extrapolate(tau_out[indl], utils::head(qv, 2))
+  }
+  if (any(indr)) {
+    qvals_out[indr] <- tail_extrapolate(tau_out[indr], utils::tail(qv, 2))
+  }
+  qvals_out
+}
+
+logit <- function(p) {
+  p <- pmax(pmin(p, 1), 0)
+  log(p) - log(1 - p)
+}
+
+# extrapolates linearly on the logistic scale using
+# the two points nearest the tail
+tail_extrapolate <- function(tau_out, qv) {
+  if (nrow(qv) == 1L) {
+    return(rep(qv$v[1], length(tau_out)))
+  }
+  x <- logit(qv$q)
+  x0 <- logit(tau_out)
+  y <- qv$v
+  m <- diff(y) / diff(x)
+  m * (x0 - x[1]) + y[1]
+}

From f97f0a084007f6a6fedc84663f7a33a99acd0929 Mon Sep 17 00:00:00 2001
From: "Daniel J. McDonald" <dajmcdon@gmail.com>
Date: Mon, 7 Jul 2025 10:58:40 -0700
Subject: [PATCH 5/8] potentially viable impl

---
 R/quant-weighted_interval.R | 111 +++++++++++++++++++++---------------
 1 file changed, 65 insertions(+), 46 deletions(-)

diff --git a/R/quant-weighted_interval.R b/R/quant-weighted_interval.R
index 47187fe5..545c00da 100644
--- a/R/quant-weighted_interval.R
+++ b/R/quant-weighted_interval.R
@@ -12,7 +12,7 @@
 #' @param quantile_levels probabilities. If specified, the score will be
 #'   computed at this set of levels. Otherwise, those present in `x` will be
 #'   used.
-#' @param na_handling character. Determines missing values are handled.
+#' @param na_handling character. Determines how missing values are handled.
 #'   For `"impute"`, missing values will be
 #'   calculated if possible using the available quantiles. For `"drop"`,
 #'   explicitly missing values are ignored in the calculation of the score, but
@@ -50,60 +50,79 @@
 #' )
 #'
 #'
-#' # Using some actual forecasts --------
-#' library(dplyr)
-#' training <- covid_case_death_rates %>%
-#'   filter(time_value >= "2021-10-01", time_value <= "2021-12-01")
-#' preds <- flatline_forecaster(
-#'   training, "death_rate",
-#'   flatline_args_list(quantile_levels = c(.01, .025, 1:19 / 20, .975, .99))
-#' )$predictions
-#' actuals <- covid_case_death_rates %>%
-#'   filter(time_value == as.Date("2021-12-01") + 7) %>%
-#'   select(geo_value, time_value, actual = death_rate)
-#' preds <- left_join(preds, actuals,
-#'   by = c("target_date" = "time_value", "geo_value")
-#' ) %>%
-#'   mutate(wis = weighted_interval_score(.pred_distn, actual))
-#' preds
-weighted_interval_score <- function(
-    x,
-    actual,
-    quantile_levels = NULL,
-    na_handling = c("impute", "drop", "propagate", "fail"),
-    ...) {
-  UseMethod("weighted_interval_score")
+wis <- function(data,...){
+  UseMethod("wis")
 }
-
+wis <- new_numeric_metric(
+  mae,
+  direction = "minimize"
+)
 
 #' @export
-weighted_interval_score.quantile_pred <- function(
-    x, actual,
+#' @rdname wis
+wis_vec <- function(
+    truth,
+    estimate,
     quantile_levels = NULL,
-    na_handling = c("impute", "drop", "propagate", "fail"),
-    ...) {
-  rlang::check_dots_empty()
-  n <- vctrs::vec_size(x)
-  if (length(actual) == 1L) actual <- rep(actual, n)
-  assert_numeric(actual, finite = TRUE, len = n)
-  assert_numeric(quantile_levels, lower = 0, upper = 1, null.ok = TRUE)
-  na_handling <- rlang::arg_match(na_handling)
-  old_quantile_levels <- x %@% "quantile_levels"
-  if (na_handling == "fail") {
-    if (is.null(quantile_levels)) {
-      cli_abort('`na_handling = "fail"` requires `quantile_levels` to be specified.')
-    }
-    if (!all(quantile_levels %in% old_quantile_levels)) {
-      return(rep(NA_real_, n))
+    na_rm = TRUE,
+    na_impute = c("none", "explicit", "any"),
+    case_weights = NULL,
+    ...
+) {
+  check_quantile_metric(truth, estimate, case_weights)
+  estimate_quantile_levels <- hardhat::extract_quantile_levels(estimate)
+  na_impute <- rlang::arg_match(na_impute)
+
+  # If we requested particular levels, but we don't impute, and the levels
+  # aren't all there, then return NA
+  if (!is.null(quantile_levels)) {
+    hardhat::check_quantile_levels(quantile_levels)
+    if (na_impute == "none" && !(all(quantile_levels %in% estimate_quantile_levels))) {
+      return(NA_real_)
     }
   }
-  tau <- quantile_levels %||% old_quantile_levels
-  x <- extrapolate_quantiles(x, tau, replace_na = (na_handling == "impute"))
-  x <- as.matrix(x)[, attr(x, "quantile_levels") %in% tau, drop = FALSE]
+  quantile_levels <- quantile_levels %||% estimate_quantile_levels
+
+  if (na_rm) {
+    result <- yardstick_remove_missing(truth, estimate, case_weights)
+
+    truth <- result$truth
+    estimate <- result$estimate
+    case_weights <- result$case_weights
+  } else if (yardstick_any_missing(truth, estimate, case_weights)) {
+    return(NA_real_)
+  }
+
+  wis_impl(
+    truth = truth,
+    estimate = estimate,
+    quantile_levels = quantile_levels,
+    na_impute = na_impute,
+    case_weights = case_weights
+  )
+}
+
+wis_impl <- function(
+    truth,
+    estimate,
+    quantile_levels = NULL,
+    na_impute = c("none", "explicit", "any"),
+    case_weights = NULL,
+    ...
+) {
+  estimate <- impute_quantiles(estimate, quantile_levels, replace_na = (na_handling == "impute"))
+  estimate <- as.matrix(estimate)[, attr(x, "quantile_levels") %in% quantile_levels, drop = FALSE]
   na_rm <- (na_handling == "drop")
-  map2_dbl(vctrs::vec_chop(x), actual, ~ wis_one_quantile(.x, tau, .y, na_rm))
+  errors <- as.vector(mapply(
+    FUN = function(.x, .y) wis_one_quantile(.x, quantile_levels, .y, na_rm),
+    vctrs::vec_chop(estimate),
+    truth
+  ), "double")
+
+  yardstick_mean(vec_score, case_weights = case_weights)
 }
 
+
 wis_one_quantile <- function(q, tau, actual, na_rm) {
   2 * mean(pmax(tau * (actual - q), (1 - tau) * (q - actual)), na.rm = na_rm)
 }

From ed519db67e32100fbd95c09a6d2491669510ee78 Mon Sep 17 00:00:00 2001
From: "Daniel J. McDonald" <dajmcdon@gmail.com>
Date: Wed, 3 Dec 2025 11:18:06 -0800
Subject: [PATCH 6/8] use quant methods in dev hardhat

---
 DESCRIPTION                   |   5 +-
 NAMESPACE                     |   3 +
 R/interpolate-quantile_pred.R | 149 ----------------------------------
 man/check_metric.Rd           |   5 ++
 man/wis.Rd                    |  75 +++++++++++++++++
 5 files changed, 86 insertions(+), 151 deletions(-)
 delete mode 100644 R/interpolate-quantile_pred.R
 create mode 100644 man/wis.Rd

diff --git a/DESCRIPTION b/DESCRIPTION
index a2dbbcb0..7f0ad675 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -23,7 +23,7 @@ Imports:
     cli,
     dplyr (>= 1.1.0),
     generics (>= 0.1.2),
-    hardhat (>= 1.3.0),
+    hardhat (>= 1.4.2.9000),
     lifecycle (>= 1.0.3),
     rlang (>= 1.1.4),
     tibble,
@@ -49,7 +49,7 @@ Config/usethis/last-upkeep: 2025-04-24
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 Collate: 
     'aaa-metrics.R'
     'import-standalone-types-check.R'
@@ -116,6 +116,7 @@ Collate:
     'prob-roc_aunp.R'
     'prob-roc_aunu.R'
     'prob-roc_curve.R'
+    'quant-weighted_interval.R'
     'reexports.R'
     'surv-brier_survival.R'
     'surv-brier_survival_integrated.R'
diff --git a/NAMESPACE b/NAMESPACE
index acddadbc..df108f58 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -149,6 +149,7 @@ export(check_dynamic_survival_metric)
 export(check_numeric_metric)
 export(check_ordered_prob_metric)
 export(check_prob_metric)
+export(check_quantile_metric)
 export(check_static_survival_metric)
 export(class_metric_summarizer)
 export(classification_cost)
@@ -262,6 +263,8 @@ export(specificity_vec)
 export(static_survival_metric_summarizer)
 export(tidy)
 export(validate_estimator)
+export(wis)
+export(wis_vec)
 export(yardstick_any_missing)
 export(yardstick_remove_missing)
 import(rlang)
diff --git a/R/interpolate-quantile_pred.R b/R/interpolate-quantile_pred.R
deleted file mode 100644
index e44f38b3..00000000
--- a/R/interpolate-quantile_pred.R
+++ /dev/null
@@ -1,149 +0,0 @@
-# Move these to hardhat? --------------------------------------------------
-
-#' restrict various objects to the interval \[lower, upper\]
-#' @param x the object to restrict
-#' @param lower numeric, the lower bound
-#' @param upper numeric, the upper bound
-#' @param ... unused
-#' @export
-#' @keywords internal
-snap <- function(x, lower, upper, ...) {
-  UseMethod("snap")
-}
-
-#' @export
-snap.default <- function(x, lower, upper, ...) {
-  rlang::check_dots_empty()
-  check_number_decimal(lower)
-  check_number_decimal(upper)
-  if (!is.numeric(x)) {
-    cli::cli_abort(
-      "{.arg x} should be should be numeric,
-      not {.obj_type_friendly {x}}."
-    )
-  }
-  pmin(pmax(x, lower), upper)
-}
-
-#' @export
-snap.quantile_pred <- function(x, lower, upper, ...) {
-  values <- as.matrix(x)
-  quantile_levels <- x %@% "quantile_levels"
-  values <- lapply(vctrs::vec_chop(values), function(.x) snap(.x, lower, upper))
-  hardhat::quantile_pred(do.call(rbind, values), quantile_levels = quantile_levels)
-}
-
-
-impute_quantiles <- function(
-    x, probs = seq(0, 1, 0.25), na.rm = FALSE, lower = -Inf,
-    upper = Inf, middle = c("cubic", "linear"), ...)
-{
-  check_bool(na.rm)
-  if (is.unsorted(probs)) probs <- sort(probs)
-  hardhat::check_quantile_levels(probs)
-  check_number_decimal(lower, allow_na = na.rm)
-  check_number_decimal(upper, allow_na = na.rm)
-  if (lower > upper) {
-    cli::cli_abort("`lower` must be less than `upper`.")
-  }
-  middle <- rlang::arg_match(middle)
-  snap(impute_quantile_internal(x, probs, middle), lower, upper)
-}
-
-impute_quantile_internal <- function(x, tau_out, middle) {
-  tau <- x %@% "quantile_levels"
-  qvals <- as.matrix(x)
-  if (all(tau_out %in% tau) && !anyNA(qvals)) {
-    return(qvals[, match(tau_out, tau), drop = FALSE])
-  }
-  if (length(tau) < 2) {
-    cli::cli_abort(
-      "Quantile interpolation is not possible when fewer than 2 quantiles
-      are avaliable."
-    )
-  }
-  qvals_out <- lapply(
-    vctrs::vec_chop(qvals),
-    function(.x) impute_quantiles_single(.x, tau, tau_out, middle)
-  )
-  qvals_out <- do.call(rbind, qvals_out)
-  qvals_out
-}
-
-impute_quantiles_single <- function(qvals, tau, tau_out, middle) {
-  qvals_out <- rep(NA, length(tau_out))
-  good <- !is.na(qvals)
-  if (!any(good)) {
-    return(qvals_out)
-  }
-  qvals <- qvals[good]
-  tau <- tau[good]
-
-  # in case we only have one point, and it matches something we wanted
-  if (length(good) < 2) {
-    matched_one <- tau_out %in% tau
-    qvals_out[matched_one] <- qvals[matched_one]
-    return(qvals_out)
-  }
-
-  indl <- tau_out < min(tau)
-  indr <- tau_out > max(tau)
-  indm <- !indl & !indr
-
-  if (middle == "cubic") {
-    method <- "cubic"
-    result <- tryCatch(
-      {
-        Q <- stats::splinefun(tau, qvals, method = "hyman")
-        quartiles <- Q(c(.25, .5, .75))
-      },
-      error = function(e) {
-        return(NA)
-      }
-    )
-  }
-  if (middle == "linear" || any(is.na(result))) {
-    method <- "linear"
-    quartiles <- stats::approx(tau, qvals, c(.25, .5, .75))$y
-  }
-  if (any(indm)) {
-    qvals_out[indm] <- switch(
-      method,
-      linear = stats::approx(tau, qvals, tau_out[indm])$y,
-      cubic = Q(tau_out[indm])
-    )
-  }
-  if (any(indl) || any(indr)) {
-    qv <- data.frame(
-      q = c(tau, tau_out[indm]),
-      v = c(qvals, qvals_out[indm])
-    ) %>%
-      dplyr::distinct(q, .keep_all = TRUE) %>%
-      dplyr::arrange(q)
-  }
-  if (any(indl)) {
-    qvals_out[indl] <- tail_extrapolate(tau_out[indl], utils::head(qv, 2))
-  }
-  if (any(indr)) {
-    qvals_out[indr] <- tail_extrapolate(tau_out[indr], utils::tail(qv, 2))
-  }
-  qvals_out
-}
-
-logit <- function(p) {
-  p <- pmax(pmin(p, 1), 0)
-  log(p) - log(1 - p)
-}
-
-# extrapolates linearly on the logistic scale using
-# the two points nearest the tail
-tail_extrapolate <- function(tau_out, qv) {
-  if (nrow(qv) == 1L) {
-    return(rep(qv$v[1], length(tau_out)))
-  }
-  x <- logit(qv$q)
-  x0 <- logit(tau_out)
-  y <- qv$v
-  m <- diff(y) / diff(x)
-  m * (x0 - x[1]) + y[1]
-}
diff --git a/man/check_metric.Rd b/man/check_metric.Rd
index 9281ae66..2c574037 100644
--- a/man/check_metric.Rd
+++ b/man/check_metric.Rd
@@ -8,6 +8,7 @@
 \alias{check_ordered_prob_metric}
 \alias{check_dynamic_survival_metric}
 \alias{check_static_survival_metric}
+\alias{check_quantile_metric}
 \title{Developer function for checking inputs in new metrics}
 \usage{
 check_numeric_metric(truth, estimate, case_weights, call = caller_env())
@@ -49,6 +50,8 @@ check_static_survival_metric(
   case_weights,
   call = caller_env()
 )
+
+check_quantile_metric(truth, estimate, case_weights, call = caller_env())
 }
 \arguments{
 \item{truth}{The realized vector of \code{truth}.
@@ -59,6 +62,7 @@ check_static_survival_metric(
 \item For \code{check_ordered_prob_metric()}, an ordered factor.
 \item For \code{check_dynamic_survival_metric()}, a Surv object.
 \item For \code{check_static_survival_metric()}, a Surv object.
+\item For \code{check_quantile_metric()}, a numeric vector.
 }}
 
 \item{estimate}{The realized \code{estimate} result.
@@ -71,6 +75,7 @@ a numeric matrix for multic-class \code{truth}.
 a numeric matrix for multic-class \code{truth}.
 \item For \code{check_dynamic_survival_metric()}, list-column of data.frames.
 \item For \code{check_static_survival_metric()}, a numeric vector.
+\item For \code{check_quantile_metric()}, a \code{hardhat::quantile_pred} vector.
 }}
 
 \item{case_weights}{The realized case weights, as a numeric vector. This must
diff --git a/man/wis.Rd b/man/wis.Rd
new file mode 100644
index 00000000..066e9856
--- /dev/null
+++ b/man/wis.Rd
@@ -0,0 +1,75 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/quant-weighted_interval.R
+\name{wis}
+\alias{wis}
+\alias{wis_vec}
+\title{Compute weighted interval score}
+\usage{
+wis(data, ...)
+
+wis_vec(
+  truth,
+  estimate,
+  quantile_levels = NULL,
+  na_rm = TRUE,
+  na_impute = c("none", "explicit", "any"),
+  case_weights = NULL,
+  ...
+)
+}
+\arguments{
+\item{...}{not used}
+
+\item{quantile_levels}{probabilities. If specified, the score will be
+computed at this set of levels. Otherwise, those present in \code{x} will be
+used.}
+
+\item{x}{A vector of class \code{quantile_pred}.}
+
+\item{actual}{double. Actual value(s)}
+
+\item{na_handling}{character. Determines how missing values are handled.
+For \code{"impute"}, missing values will be
+calculated if possible using the available quantiles. For \code{"drop"},
+explicitly missing values are ignored in the calculation of the score, but
+implicitly missing values are imputed if possible.
+For \code{"propogate"}, the resulting score will be \code{NA} if any missing values
+exist. Finally, if
+\code{quantile_levels} is specified, \code{"fail"} will result in
+the score being \code{NA} when any required quantile levels (implicit or explicit)
+do not have corresponding values.}
+}
+\value{
+a vector of nonnegative scores.
+}
+\description{
+Weighted interval score (WIS), a well-known quantile-based
+approximation of the commonly-used continuous ranked probability score
+(CRPS). WIS is a proper score, and can be thought of as a distributional
+generalization of absolute error. For example, see \href{https://arxiv.org/abs/2005.12881}{Bracher et al. (2020)} for discussion in the context
+of COVID-19 forecasting.
+}
+\examples{
+quantile_levels <- c(.2, .4, .6, .8)
+predq1 <- 1:4 #
+predq2 <- 8:11
+dstn <- quantile_pred(rbind(predq1, predq2), quantile_levels)
+actual <- c(3.3, 7.1)
+weighted_interval_score(dstn, actual)
+weighted_interval_score(dstn, actual, c(.25, .5, .75))
+
+# Missing value behaviours
+dstn <- quantile_pred(matrix(c(1, 2, NA, 4), nrow = 1), 1:4 / 5)
+weighted_interval_score(dstn, 2.5)
+weighted_interval_score(dstn, 2.5, 1:9 / 10)
+weighted_interval_score(dstn, 2.5, 1:9 / 10, na_handling = "drop")
+weighted_interval_score(dstn, 2.5, na_handling = "propagate")
+weighted_interval_score(
+  quantile_pred(matrix(1:4, nrow = 1), 1:4 / 5),
+  actual = 2.5,
+  quantile_levels = 1:9 / 10,
+  na_handling = "fail"
+)
+
+
+}

From df5a8cd1b9d65a18707bbcf3c507aca8b9deb9a4 Mon Sep 17 00:00:00 2001
From: "Daniel J. McDonald" <dajmcdon@gmail.com>
Date: Wed, 3 Dec 2025 11:32:11 -0800
Subject: [PATCH 7/8] rename wis fun

---
 NAMESPACE                      |  4 +-
 R/quant-weighted_interval.R    | 65 +++++++++++++++++------------
 man/weighted_interval_score.Rd | 75 ++++++++++++++++++++++++++++++++++
 3 files changed, 116 insertions(+), 28 deletions(-)
 create mode 100644 man/weighted_interval_score.Rd

diff --git a/NAMESPACE b/NAMESPACE
index df108f58..ad601a76 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -263,8 +263,8 @@ export(specificity_vec)
 export(static_survival_metric_summarizer)
 export(tidy)
 export(validate_estimator)
-export(wis)
-export(wis_vec)
+export(weighted_interval_score)
+export(weighted_interval_score_vec)
 export(yardstick_any_missing)
 export(yardstick_remove_missing)
 import(rlang)
diff --git a/R/quant-weighted_interval.R b/R/quant-weighted_interval.R
index 545c00da..000a9d13 100644
--- a/R/quant-weighted_interval.R
+++ b/R/quant-weighted_interval.R
@@ -50,24 +50,24 @@
 #' )
 #'
 #'
-wis <- function(data,...){
-  UseMethod("wis")
+weighted_interval_score <- function(data, ...) {
+  UseMethod("weighted_interval_score")
 }
-wis <- new_numeric_metric(
+weighted_interval_score <- new_numeric_metric(
   mae,
   direction = "minimize"
 )
 
 #' @export
-#' @rdname wis
-wis_vec <- function(
-    truth,
-    estimate,
-    quantile_levels = NULL,
-    na_rm = TRUE,
-    na_impute = c("none", "explicit", "any"),
-    case_weights = NULL,
-    ...
+#' @rdname weighted_interval_score
+weighted_interval_score_vec <- function(
+  truth,
+  estimate,
+  quantile_levels = NULL,
+  na_rm = TRUE,
+  na_impute = c("none", "explicit", "any"),
+  case_weights = NULL,
+  ...
 ) {
   check_quantile_metric(truth, estimate, case_weights)
   estimate_quantile_levels <- hardhat::extract_quantile_levels(estimate)
@@ -77,7 +77,10 @@ wis_vec <- function(
   # aren't all there, then return NA
   if (!is.null(quantile_levels)) {
     hardhat::check_quantile_levels(quantile_levels)
-    if (na_impute == "none" && !(all(quantile_levels %in% estimate_quantile_levels))) {
+    if (
+      na_impute == "none" &&
+        !(all(quantile_levels %in% estimate_quantile_levels))
+    ) {
       return(NA_real_)
     }
   }
@@ -103,21 +106,31 @@ wis_vec <- function(
 }
 
 wis_impl <- function(
-    truth,
-    estimate,
-    quantile_levels = NULL,
-    na_impute = c("none", "explicit", "any"),
-    case_weights = NULL,
-    ...
+  truth,
+  estimate,
+  quantile_levels = NULL,
+  na_impute = c("none", "explicit", "any"),
+  case_weights = NULL,
+  ...
 ) {
-  estimate <- impute_quantiles(estimate, quantile_levels, replace_na = (na_handling == "impute"))
-  estimate <- as.matrix(estimate)[, attr(x, "quantile_levels") %in% quantile_levels, drop = FALSE]
+  estimate <- impute_quantiles(
+    estimate,
+    quantile_levels,
+    replace_na = (na_handling == "impute")
+  )
+  estimate <- as.matrix(estimate)[,
+    attr(x, "quantile_levels") %in% quantile_levels,
+    drop = FALSE
+  ]
   na_rm <- (na_handling == "drop")
-  errors <- as.vector(mapply(
-    FUN = function(.x, .y) wis_one_quantile(.x, quantile_levels, .y, na_rm),
-    vctrs::vec_chop(estimate),
-    truth
-  ), "double")
+  errors <- as.vector(
+    mapply(
+      FUN = function(.x, .y) wis_one_quantile(.x, quantile_levels, .y, na_rm),
+      vctrs::vec_chop(estimate),
+      truth
+    ),
+    "double"
+  )
 
   yardstick_mean(vec_score, case_weights = case_weights)
 }
diff --git a/man/weighted_interval_score.Rd b/man/weighted_interval_score.Rd
new file mode 100644
index 00000000..7b0619bc
--- /dev/null
+++ b/man/weighted_interval_score.Rd
@@ -0,0 +1,75 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/quant-weighted_interval.R
+\name{weighted_interval_score}
+\alias{weighted_interval_score}
+\alias{weighted_interval_score_vec}
+\title{Compute weighted interval score}
+\usage{
+weighted_interval_score(data, ...)
+
+weighted_interval_score_vec(
+  truth,
+  estimate,
+  quantile_levels = NULL,
+  na_rm = TRUE,
+  na_impute = c("none", "explicit", "any"),
+  case_weights = NULL,
+  ...
+)
+}
+\arguments{
+\item{...}{not used}
+
+\item{quantile_levels}{probabilities. If specified, the score will be
+computed at this set of levels. Otherwise, those present in \code{x} will be
+used.}
+
+\item{x}{A vector of class \code{quantile_pred}.}
+
+\item{actual}{double. Actual value(s)}
+
+\item{na_handling}{character. Determines how missing values are handled.
+For \code{"impute"}, missing values will be
+calculated if possible using the available quantiles. For \code{"drop"},
+explicitly missing values are ignored in the calculation of the score, but
+implicitly missing values are imputed if possible.
+For \code{"propogate"}, the resulting score will be \code{NA} if any missing values
+exist. Finally, if
+\code{quantile_levels} is specified, \code{"fail"} will result in
+the score being \code{NA} when any required quantile levels (implicit or explicit)
+do not have corresponding values.}
+}
+\value{
+a vector of nonnegative scores.
+}
+\description{
+Weighted interval score (WIS), a well-known quantile-based
+approximation of the commonly-used continuous ranked probability score
+(CRPS). WIS is a proper score, and can be thought of as a distributional
+generalization of absolute error. For example, see \href{https://arxiv.org/abs/2005.12881}{Bracher et al. (2020)} for discussion in the context
+of COVID-19 forecasting.
+}
+\examples{
+quantile_levels <- c(.2, .4, .6, .8)
+predq1 <- 1:4 #
+predq2 <- 8:11
+dstn <- quantile_pred(rbind(predq1, predq2), quantile_levels)
+actual <- c(3.3, 7.1)
+weighted_interval_score(dstn, actual)
+weighted_interval_score(dstn, actual, c(.25, .5, .75))
+
+# Missing value behaviours
+dstn <- quantile_pred(matrix(c(1, 2, NA, 4), nrow = 1), 1:4 / 5)
+weighted_interval_score(dstn, 2.5)
+weighted_interval_score(dstn, 2.5, 1:9 / 10)
+weighted_interval_score(dstn, 2.5, 1:9 / 10, na_handling = "drop")
+weighted_interval_score(dstn, 2.5, na_handling = "propagate")
+weighted_interval_score(
+  quantile_pred(matrix(1:4, nrow = 1), 1:4 / 5),
+  actual = 2.5,
+  quantile_levels = 1:9 / 10,
+  na_handling = "fail"
+)
+
+
+}

From 87b5f36eff22f3228d7f2195cbb6c62527d74474 Mon Sep 17 00:00:00 2001
From: "Daniel J. McDonald" <dajmcdon@gmail.com>
Date: Wed, 3 Dec 2025 15:12:25 -0800
Subject: [PATCH 8/8] fix NA handling

---
 R/quant-weighted_interval.R    | 167 ++++++++++++++++++++-------------
 man/weighted_interval_score.Rd |  90 ++++++++++++------
 man/wis.Rd                     |  75 ---------------
 3 files changed, 160 insertions(+), 172 deletions(-)
 delete mode 100644 man/wis.Rd

diff --git a/R/quant-weighted_interval.R b/R/quant-weighted_interval.R
index 000a9d13..f519a73e 100644
--- a/R/quant-weighted_interval.R
+++ b/R/quant-weighted_interval.R
@@ -11,17 +11,33 @@
 #' @param actual double. Actual value(s)
 #' @param quantile_levels probabilities. If specified, the score will be
 #'   computed at this set of levels. Otherwise, those present in `x` will be
-#'   used.
-#' @param na_handling character. Determines how missing values are handled.
-#'   For `"impute"`, missing values will be
-#'   calculated if possible using the available quantiles. For `"drop"`,
-#'   explicitly missing values are ignored in the calculation of the score, but
-#'   implicitly missing values are imputed if possible.
-#'   For `"propogate"`, the resulting score will be `NA` if any missing values
-#'   exist. Finally, if
-#'   `quantile_levels` is specified, `"fail"` will result in
-#'   the score being `NA` when any required quantile levels (implicit or explicit)
-#'   do not have corresponding values.
+#'   used. If `quantile_levels` do not exactly match those available in `x`,
+#'   then some quantiles will have implicit missingness. Handling of these
+#'   is determined by `quantile_estimate_nas`.
+#' @param quantile_estimate_nas character. This argument applies only to `x`.
+#'   It handles imputation of individual `quantile_levels` that are necessary to
+#'   compute a score. Because each element of `x` is a [hardhat::quantile_pred],
+#'   it is possible for these to be missing for particular
+#'   `quantile_levels`. There are a number of different possibilities for such
+#'   missingness. The options are as follows:
+#'   * For `"impute"`, both explicit and implicit missing values will be imputed
+#'   using [hardhat::impute_quantiles()] prior to the calculation of the score.
+#'   So the score will be `NA` only if imputation fails.
+#'   * For `"drop"`, any explicit missing values will be removed
+#'   before calculating the score for a particular prediction. This may be
+#'   reasonable due to the weighting. For example, if the estimate has
+#'   `quantile_levels = c(.25, .5, .75)` but the median is `NA` for a particular
+#'   prediction, it may be reasonable to average the accuracy of `c(.25, .75)`
+#'   for that prediction with others that don't have missingness. This option
+#'   is only works if `quantile_levels = NULL` or is a subset of the
+#'   `quantile_levels` in `x`.
+#'   * For `"propagate"`, any missing value predictions will result in that
+#'   element of `x` having a score of `NA`. If `na_rm = TRUE`, then these will
+#'   be removed before averaging.
+#' @param na_rm logical. If `TRUE`, missing values in `actual` or both implicit and
+#'   explicit (values of `NA` present in `x`), will be ignored (dropped) in the
+#'   calculation of the summary score. If `FALSE` (the default), any `NA`s will
+#'   result in the summary being `NA`.
 #' @param ... not used
 #'
 #' @return a vector of nonnegative scores.
@@ -29,26 +45,37 @@
 #' @export
 #' @examples
 #' quantile_levels <- c(.2, .4, .6, .8)
-#' predq1 <- 1:4 #
-#' predq2 <- 8:11
-#' dstn <- quantile_pred(rbind(predq1, predq2), quantile_levels)
-#' actual <- c(3.3, 7.1)
-#' weighted_interval_score(dstn, actual)
-#' weighted_interval_score(dstn, actual, c(.25, .5, .75))
+#' pred1 <- 1:4
+#' pred2 <- 8:11
+#' preds <- quantile_pred(rbind(pred1, pred2), quantile_levels)
+#' truth <- c(3.3, 7.1)
+#' weighted_interval_score_vec(truth, preds)
+#' weighted_interval_score_vec(truth, preds, quantile_levels = c(.25, .5, .75))
 #'
 #' # Missing value behaviours
-#' dstn <- quantile_pred(matrix(c(1, 2, NA, 4), nrow = 1), 1:4 / 5)
-#' weighted_interval_score(dstn, 2.5)
-#' weighted_interval_score(dstn, 2.5, 1:9 / 10)
-#' weighted_interval_score(dstn, 2.5, 1:9 / 10, na_handling = "drop")
-#' weighted_interval_score(dstn, 2.5, na_handling = "propagate")
-#' weighted_interval_score(
-#'   quantile_pred(matrix(1:4, nrow = 1), 1:4 / 5),
-#'   actual = 2.5,
+#'
+#' preds_na <- quantile_pred(rbind(pred1, c(1, 2, NA, 4)), 1:4 / 5)
+#' truth <- c(2.5, 2.5)
+#' weighted_interval_score_vec(truth, preds_na)
+#' weighted_interval_score_vec(truth, preds_na, quantile_levels = 1:9 / 10)
+#' expect_error(weighted_interval_score_vec(
+#'   truth,
+#'   preds_na,
 #'   quantile_levels = 1:9 / 10,
-#'   na_handling = "fail"
+#'   quantile_estimate_nas = "drop"
+#' ))
+#' weighted_interval_score_vec(
+#'   truth,
+#'   preds_na,
+#'   quantile_levels = c(2, 3) / 5,
+#'   quantile_estimate_nas = "drop"
+#' )
+#' weighted_interval_score_vec(
+#'   truth, preds_na, na_rm = TRUE, quantile_estimate_nas = "propagate"
+#' )
+#' weighted_interval_score_vec(
+#'   truth, preds_na, quantile_estimate_nas = "propagate"
 #' )
-#'
 #'
 weighted_interval_score <- function(data, ...) {
   UseMethod("weighted_interval_score")
@@ -64,78 +91,84 @@ weighted_interval_score_vec <- function(
   truth,
   estimate,
   quantile_levels = NULL,
-  na_rm = TRUE,
-  na_impute = c("none", "explicit", "any"),
+  na_rm = FALSE,
+  quantile_estimate_nas = c("impute", "drop", "propagate"),
   case_weights = NULL,
   ...
 ) {
   check_quantile_metric(truth, estimate, case_weights)
   estimate_quantile_levels <- hardhat::extract_quantile_levels(estimate)
-  na_impute <- rlang::arg_match(na_impute)
-
-  # If we requested particular levels, but we don't impute, and the levels
-  # aren't all there, then return NA
+  quantile_estimate_nas <- rlang::arg_match(quantile_estimate_nas)
   if (!is.null(quantile_levels)) {
     hardhat::check_quantile_levels(quantile_levels)
-    if (
-      na_impute == "none" &&
-        !(all(quantile_levels %in% estimate_quantile_levels))
-    ) {
+    all_levels_estimated <- all(quantile_levels %in% estimate_quantile_levels)
+    if (quantile_estimate_nas == "drop" && !all_levels_estimated) {
+      cli::cli_abort(
+        "When `quantile_levels` is not a subset of those available in `estimate`, 
+      `quantile_estimate_nas` may not be `'drop'`."
+      )
+    }
+    if (!all_levels_estimated && (quantile_estimate_nas == "propagate")) {
+      # We requested particular levels, but the levels aren't all there,
+      # and NAs propagate, so return NA
       return(NA_real_)
     }
   }
+
   quantile_levels <- quantile_levels %||% estimate_quantile_levels
+  if (quantile_estimate_nas %in% c("drop", "propagate")) {
+    levels_estimated <- estimate_quantile_levels %in% quantile_levels
+    estimate <- as.matrix(estimate)[, levels_estimated, drop = FALSE]
+  } else {
+    estimate <- as.matrix(hardhat::impute_quantiles(estimate, quantile_levels))
+  }
+
+  vec_wis <- wis_impl(
+    truth = truth,
+    estimate = estimate,
+    quantile_levels = quantile_levels,
+    rowwise_na_rm = (quantile_estimate_nas == "drop")
+  )
 
   if (na_rm) {
-    result <- yardstick_remove_missing(truth, estimate, case_weights)
+    result <- yardstick_remove_missing(truth, vec_wis, case_weights)
 
     truth <- result$truth
-    estimate <- result$estimate
+    vec_wis <- result$estimate
     case_weights <- result$case_weights
-  } else if (yardstick_any_missing(truth, estimate, case_weights)) {
+  } else if (yardstick_any_missing(truth, vec_wis, case_weights)) {
     return(NA_real_)
   }
 
-  wis_impl(
-    truth = truth,
-    estimate = estimate,
-    quantile_levels = quantile_levels,
-    na_impute = na_impute,
-    case_weights = case_weights
-  )
+  yardstick_mean(vec_wis, case_weights = case_weights)
 }
 
 wis_impl <- function(
   truth,
   estimate,
-  quantile_levels = NULL,
-  na_impute = c("none", "explicit", "any"),
-  case_weights = NULL,
-  ...
+  quantile_levels,
+  rowwise_na_rm = TRUE
 ) {
-  estimate <- impute_quantiles(
-    estimate,
-    quantile_levels,
-    replace_na = (na_handling == "impute")
-  )
-  estimate <- as.matrix(estimate)[,
-    attr(x, "quantile_levels") %in% quantile_levels,
-    drop = FALSE
-  ]
-  na_rm <- (na_handling == "drop")
-  errors <- as.vector(
+  as.vector(
     mapply(
-      FUN = function(.x, .y) wis_one_quantile(.x, quantile_levels, .y, na_rm),
+      FUN = function(.x, .y) {
+        wis_one_quantile(.x, quantile_levels, .y, rowwise_na_rm)
+      },
       vctrs::vec_chop(estimate),
       truth
     ),
     "double"
   )
-
-  yardstick_mean(vec_score, case_weights = case_weights)
 }
 
 
-wis_one_quantile <- function(q, tau, actual, na_rm) {
-  2 * mean(pmax(tau * (actual - q), (1 - tau) * (q - actual)), na.rm = na_rm)
+wis_one_quantile <- function(values, quantile_levels, truth, na_rm) {
+  2 *
+    mean(
+      pmax(
+        quantile_levels * (truth - values),
+        (1 - quantile_levels) * (values - truth)
+      ),
+      na.rm = na_rm
+    )
 }
diff --git a/man/weighted_interval_score.Rd b/man/weighted_interval_score.Rd
index 7b0619bc..5791decc 100644
--- a/man/weighted_interval_score.Rd
+++ b/man/weighted_interval_score.Rd
@@ -11,8 +11,8 @@ weighted_interval_score_vec(
   truth,
   estimate,
   quantile_levels = NULL,
-  na_rm = TRUE,
-  na_impute = c("none", "explicit", "any"),
+  na_rm = FALSE,
+  quantile_estimate_nas = c("impute", "drop", "propagate"),
   case_weights = NULL,
   ...
 )
@@ -22,22 +22,41 @@ weighted_interval_score_vec(
 
 \item{quantile_levels}{probabilities. If specified, the score will be
 computed at this set of levels. Otherwise, those present in \code{x} will be
-used.}
+used. If \code{quantile_levels} do not exactly match those available in \code{x},
+then some quantiles will have implicit missingness. Handling of these
+is determined by \code{quantile_estimate_nas}.}
+
+\item{na_rm}{logical. If \code{TRUE}, missing values in \code{actual} or both implicit and
+explicit (values of \code{NA} present in \code{x}), will be ignored (dropped) in the
+calculation of the summary score. If \code{FALSE} (the default), any \code{NA}s will
+result in the summary being \code{NA}.}
+
+\item{quantile_estimate_nas}{character. This argument applies only to \code{x}.
+It handles imputation of individual \code{quantile_levels} that are necessary to
+compute a score. Because each element of \code{x} is a \link[hardhat:quantile_pred]{hardhat::quantile_pred},
+it is possible for these to be missing for particular
+\code{quantile_levels}. There are a number of different possibilities for such
+missingness. The options are as follows:
+\itemize{
+\item For \code{"impute"}, both explicit and implicit missing values will be imputed
+using \code{\link[hardhat:impute_quantiles]{hardhat::impute_quantiles()}} prior to the calculation of the score.
+So the score will be \code{NA} only if imputation fails.
+\item For \code{"drop"}, any explicit missing values will be removed
+before calculating the score for a particular prediction. This may be
+reasonable due to the weighting. For example, if the estimate has
+\code{quantile_levels = c(.25, .5, .75)} but the median is \code{NA} for a particular
+prediction, it may be reasonable to average the accuracy of \code{c(.25, .75)}
+for that prediction with others that don't have missingness. This option
+is only works if \code{quantile_levels = NULL} or is a subset of the
+\code{quantile_levels} in \code{x}.
+\item For \code{"propagate"}, any missing value predictions will result in that
+element of \code{x} having a score of \code{NA}. If \code{na_rm = TRUE}, then these will
+be removed before averaging.
+}}
 
 \item{x}{A vector of class \code{quantile_pred}.}
 
 \item{actual}{double. Actual value(s)}
-
-\item{na_handling}{character. Determines how missing values are handled.
-For \code{"impute"}, missing values will be
-calculated if possible using the available quantiles. For \code{"drop"},
-explicitly missing values are ignored in the calculation of the score, but
-implicitly missing values are imputed if possible.
-For \code{"propogate"}, the resulting score will be \code{NA} if any missing values
-exist. Finally, if
-\code{quantile_levels} is specified, \code{"fail"} will result in
-the score being \code{NA} when any required quantile levels (implicit or explicit)
-do not have corresponding values.}
 }
 \value{
 a vector of nonnegative scores.
@@ -51,25 +70,36 @@ of COVID-19 forecasting.
 }
 \examples{
 quantile_levels <- c(.2, .4, .6, .8)
-predq1 <- 1:4 #
-predq2 <- 8:11
-dstn <- quantile_pred(rbind(predq1, predq2), quantile_levels)
-actual <- c(3.3, 7.1)
-weighted_interval_score(dstn, actual)
-weighted_interval_score(dstn, actual, c(.25, .5, .75))
+pred1 <- 1:4
+pred2 <- 8:11
+preds <- quantile_pred(rbind(pred1, pred2), quantile_levels)
+truth <- c(3.3, 7.1)
+weighted_interval_score_vec(truth, preds)
+weighted_interval_score_vec(truth, preds, quantile_levels = c(.25, .5, .75))
 
 # Missing value behaviours
-dstn <- quantile_pred(matrix(c(1, 2, NA, 4), nrow = 1), 1:4 / 5)
-weighted_interval_score(dstn, 2.5)
-weighted_interval_score(dstn, 2.5, 1:9 / 10)
-weighted_interval_score(dstn, 2.5, 1:9 / 10, na_handling = "drop")
-weighted_interval_score(dstn, 2.5, na_handling = "propagate")
-weighted_interval_score(
-  quantile_pred(matrix(1:4, nrow = 1), 1:4 / 5),
-  actual = 2.5,
+
+preds_na <- quantile_pred(rbind(pred1, c(1, 2, NA, 4)), 1:4 / 5)
+truth <- c(2.5, 2.5)
+weighted_interval_score_vec(truth, preds_na)
+weighted_interval_score_vec(truth, preds_na, quantile_levels = 1:9 / 10)
+expect_error(weighted_interval_score_vec(
+  truth,
+  preds_na,
   quantile_levels = 1:9 / 10,
-  na_handling = "fail"
+  quantile_estimate_nas = "drop"
+))
+weighted_interval_score_vec(
+  truth,
+  preds_na,
+  quantile_levels = c(2, 3) / 5,
+  quantile_estimate_nas = "drop"
+)
+weighted_interval_score_vec(
+  truth, preds_na, na_rm = TRUE, quantile_estimate_nas = "propagate"
+)
+weighted_interval_score_vec(
+  truth, preds_na, quantile_estimate_nas = "propagate"
 )
-
 
 }
diff --git a/man/wis.Rd b/man/wis.Rd
deleted file mode 100644
index 066e9856..00000000
--- a/man/wis.Rd
+++ /dev/null
@@ -1,75 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/quant-weighted_interval.R
-\name{wis}
-\alias{wis}
-\alias{wis_vec}
-\title{Compute weighted interval score}
-\usage{
-wis(data, ...)
-
-wis_vec(
-  truth,
-  estimate,
-  quantile_levels = NULL,
-  na_rm = TRUE,
-  na_impute = c("none", "explicit", "any"),
-  case_weights = NULL,
-  ...
-)
-}
-\arguments{
-\item{...}{not used}
-
-\item{quantile_levels}{probabilities. If specified, the score will be
-computed at this set of levels. Otherwise, those present in \code{x} will be
-used.}
-
-\item{x}{A vector of class \code{quantile_pred}.}
-
-\item{actual}{double. Actual value(s)}
-
-\item{na_handling}{character. Determines how missing values are handled.
-For \code{"impute"}, missing values will be
-calculated if possible using the available quantiles. For \code{"drop"},
-explicitly missing values are ignored in the calculation of the score, but
-implicitly missing values are imputed if possible.
-For \code{"propogate"}, the resulting score will be \code{NA} if any missing values
-exist. Finally, if
-\code{quantile_levels} is specified, \code{"fail"} will result in
-the score being \code{NA} when any required quantile levels (implicit or explicit)
-do not have corresponding values.}
-}
-\value{
-a vector of nonnegative scores.
-}
-\description{
-Weighted interval score (WIS), a well-known quantile-based
-approximation of the commonly-used continuous ranked probability score
-(CRPS). WIS is a proper score, and can be thought of as a distributional
-generalization of absolute error. For example, see \href{https://arxiv.org/abs/2005.12881}{Bracher et al. (2020)} for discussion in the context
-of COVID-19 forecasting.
-}
-\examples{
-quantile_levels <- c(.2, .4, .6, .8)
-predq1 <- 1:4 #
-predq2 <- 8:11
-dstn <- quantile_pred(rbind(predq1, predq2), quantile_levels)
-actual <- c(3.3, 7.1)
-weighted_interval_score(dstn, actual)
-weighted_interval_score(dstn, actual, c(.25, .5, .75))
-
-# Missing value behaviours
-dstn <- quantile_pred(matrix(c(1, 2, NA, 4), nrow = 1), 1:4 / 5)
-weighted_interval_score(dstn, 2.5)
-weighted_interval_score(dstn, 2.5, 1:9 / 10)
-weighted_interval_score(dstn, 2.5, 1:9 / 10, na_handling = "drop")
-weighted_interval_score(dstn, 2.5, na_handling = "propagate")
-weighted_interval_score(
-  quantile_pred(matrix(1:4, nrow = 1), 1:4 / 5),
-  actual = 2.5,
-  quantile_levels = 1:9 / 10,
-  na_handling = "fail"
-)
-
-
-}