Merge pull request #365 from cmu-delphi/361-dist-quantiles-null

361 dist quantiles null

Author: dajmcdon

NAMESPACE

@@ -124,6 +124,10 @@ S3method(tidy,layer)
 S3method(update,layer)
 S3method(vec_ptype_abbr,dist_quantiles)
 S3method(vec_ptype_full,dist_quantiles)
+S3method(weighted_interval_score,default)
+S3method(weighted_interval_score,dist_default)
+S3method(weighted_interval_score,dist_quantiles)
+S3method(weighted_interval_score,distribution)
 export("%>%")
 export(Add_model)
 export(Remove_model)
@@ -207,6 +211,7 @@ export(update_epi_recipe)
 export(update_frosting)
 export(update_model)
 export(validate_layer)
+export(weighted_interval_score)
 import(distributional)
 import(epiprocess)
 import(parsnip)

NEWS.md

@@ -52,5 +52,7 @@ Pre-1.0.0 numbering scheme: 0.x will indicate releases, while 0.0.x will indicat
   `...` args intended for `predict.model_fit()`
 - `bake.epi_recipe()` will now re-infer the geo and time type in case baking the
   steps has changed the appropriate values
+- produce length 0 `dist_quantiles()`
+- add functionality to calculate weighted interval scores for `dist_quantiles()`
 - Add `step_epi_slide` to produce generic sliding computations over an `epi_df`
 - Add quantile random forests (via `{grf}`) as a parsnip engine

R/dist_quantiles.R

@@ -1,10 +1,21 @@
 #' @importFrom vctrs field vec_cast new_rcrd
-new_quantiles <- function(values = double(), quantile_levels = double()) {
+new_quantiles <- function(values = double(1), quantile_levels = double(1)) {
   arg_is_probabilities(quantile_levels)
 
   vec_cast(values, double())
   vec_cast(quantile_levels, double())
+  values <- unname(values)
+  if (length(values) == 0L) {
+    return(new_rcrd(
+      list(
+        values = rep(NA_real_, length(quantile_levels)),
+        quantile_levels = quantile_levels
+      ),
+      class = c("dist_quantiles", "dist_default")
+    ))
+  }
   stopifnot(length(values) == length(quantile_levels))
+
   stopifnot(!vctrs::vec_duplicate_any(quantile_levels))
   if (is.unsorted(quantile_levels)) {
     o <- vctrs::vec_order(quantile_levels)
@@ -37,30 +48,49 @@ format.dist_quantiles <- function(x, digits = 2, ...) {
 
 #' A distribution parameterized by a set of quantiles
 #'
-#' @param values A vector of values
-#' @param quantile_levels A vector of probabilities corresponding to `values`
+#' @param values A vector (or list of vectors) of values.
+#' @param quantile_levels A vector (or list of vectors) of probabilities
+#'   corresponding to `values`.
+#'
+#' When creating multiple sets of `values`/`quantile_levels` resulting in
+#' different distributions, the sizes must match. See the examples below.
+#'
+#' @return A vector of class `"distribution"`.
 #'
 #' @export
 #'
 #' @examples
-#' dstn <- dist_quantiles(list(1:4, 8:11), list(c(.2, .4, .6, .8)))
+#' dist_quantiles(1:4, 1:4 / 5)
+#' dist_quantiles(list(1:3, 1:4), list(1:3 / 4, 1:4 / 5))
+#' dstn <- dist_quantiles(list(1:4, 8:11), c(.2, .4, .6, .8))
+#' dstn
+#'
 #' quantile(dstn, p = c(.1, .25, .5, .9))
 #' median(dstn)
 #'
 #' # it's a bit annoying to inspect the data
 #' 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)
+  if (!is.list(quantile_levels)) {
+    assert_numeric(quantile_levels, lower = 0, upper = 1, any.missing = FALSE, min.len = 1L)
+    quantile_levels <- list(quantile_levels)
+  }
+  if (!is.list(values)) {
+    if (length(values) == 0L) values <- NA_real_
+    values <- list(values)
+  }
 
   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))
+
+  qntls <- as_list_of(
+    map2(args$values, args$quantile_levels, new_quantiles),
+    .ptype = new_quantiles(NA_real_, 0.5)
+  )
   new_vctr(qntls, class = "distribution")
 }
 
@@ -87,59 +117,6 @@ validate_dist_quantiles <- function(values, quantile_levels) {
 }
 
 
-#' Summarize a distribution with a set of quantiles
-#'
-#' @param x a `distribution` vector
-#' @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`
-#' @export
-#'
-#' @examples
-#' library(distributional)
-#' dstn <- dist_normal(c(10, 2), c(5, 10))
-#' 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, 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, probs = c(.25, 0.5, .75))
-extrapolate_quantiles <- function(x, probs, ...) {
-  UseMethod("extrapolate_quantiles")
-}
-
-#' @export
-#' @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, probs, ...) {
-  values <- quantile(x, probs, ...)
-  new_quantiles(values = values, quantile_levels = probs)
-}
-
-#' @export
-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")
 }
@@ -183,18 +160,20 @@ quantile.dist_quantiles <- function(x, p, ..., middle = c("cubic", "linear")) {
 quantile_extrapolate <- function(x, tau_out, middle) {
   tau <- field(x, "quantile_levels")
   qvals <- field(x, "values")
-  r <- range(tau, na.rm = TRUE)
+  nas <- is.na(qvals)
   qvals_out <- rep(NA, length(tau_out))
+  qvals <- qvals[!nas]
+  tau <- tau[!nas]
 
   # short circuit if we aren't actually extrapolating
   # matches to ~15 decimals
   if (all(tau_out %in% tau)) {
     return(qvals[match(tau_out, tau)])
   }
-  if (length(qvals) < 2) {
-    cli::cli_abort(c(
+  if (length(tau) < 2) {
+    cli::cli_abort(
       "Quantile extrapolation is not possible with fewer than 2 quantiles."
-    ))
+    )
     return(qvals_out)
   }
 

R/extrapolate_quantiles.R

@@ -0,0 +1,65 @@
+#' Summarize a distribution with a set of quantiles
+#'
+#' @param x a `distribution` vector
+#' @param probs a vector of probabilities at which to calculate quantiles
+#' @param replace_na logical. If `x` contains `NA`'s, these are imputed if
+#'   possible (if `TRUE`) or retained (if `FALSE`). This only effects
+#'   elements of class `dist_quantiles`.
+#' @param ... additional arguments passed on to the `quantile` method
+#'
+#' @return a `distribution` vector containing `dist_quantiles`. Any elements
+#'   of `x` which were originally `dist_quantiles` will now have a superset
+#'   of the original `quantile_values` (the union of those and `probs`).
+#' @export
+#'
+#' @examples
+#' library(distributional)
+#' dstn <- dist_normal(c(10, 2), c(5, 10))
+#' 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, 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, probs = c(.25, 0.5, .75))
+extrapolate_quantiles <- function(x, probs, replace_na = TRUE, ...) {
+  UseMethod("extrapolate_quantiles")
+}
+
+#' @export
+#' @importFrom vctrs vec_data
+extrapolate_quantiles.distribution <- function(x, probs, replace_na = TRUE, ...) {
+  rlang::check_dots_empty()
+  arg_is_lgl_scalar(replace_na)
+  arg_is_probabilities(probs)
+  if (is.unsorted(probs)) probs <- sort(probs)
+  dstn <- lapply(vec_data(x), extrapolate_quantiles, probs = probs, replace_na = replace_na)
+  new_vctr(dstn, vars = NULL, class = "distribution")
+}
+
+#' @export
+extrapolate_quantiles.dist_default <- function(x, probs, replace_na = TRUE, ...) {
+  values <- quantile(x, probs, ...)
+  new_quantiles(values = values, quantile_levels = probs)
+}
+
+#' @export
+extrapolate_quantiles.dist_quantiles <- function(x, probs, replace_na = TRUE, ...) {
+  orig_probs <- field(x, "quantile_levels")
+  orig_values <- field(x, "values")
+  new_probs <- c(orig_probs, probs)
+  dups <- duplicated(new_probs)
+  if (!replace_na || !anyNA(orig_values)) {
+    new_values <- c(orig_values, quantile(x, probs, ...))
+  } else {
+    nas <- is.na(orig_values)
+    orig_values[nas] <- quantile(x, orig_probs[nas], ...)
+    new_values <- c(orig_values, quantile(x, probs, ...))
+  }
+  new_quantiles(new_values[!dups], new_probs[!dups])
+}

R/flatline_forecaster.R

@@ -1,8 +1,9 @@
 #' Predict the future with today's value
 #'
 #' This is a simple forecasting model for
-#' [epiprocess::epi_df][epiprocess::as_epi_df] data. It uses the most recent observation as the
-#' forcast for any future date, and produces intervals based on the quantiles
+#' [epiprocess::epi_df][epiprocess::as_epi_df] data. It uses the most recent
+#' observation as the
+#' forecast for any future date, and produces intervals based on the quantiles
 #' of the residuals of such a "flatline" forecast over all available training
 #' data.
 #'