refactor residual quantiles to incorporate dist_quantiles.

Author: dajmcdon

NAMESPACE

@@ -33,8 +33,12 @@ S3method(refresh_blueprint,default_epi_recipe_blueprint)
 S3method(run_mold,default_epi_recipe_blueprint)
 S3method(slather,layer_naomit)
 S3method(slather,layer_predict)
-S3method(slather,layer_residual_quantile)
+S3method(slather,layer_residual_quantiles)
 S3method(slather,layer_threshold)
+S3method(snap,default)
+S3method(snap,dist_default)
+S3method(snap,dist_quantiles)
+S3method(snap,distribution)
 S3method(vec_ptype_abbr,dist_quantiles)
 S3method(vec_ptype_full,dist_quantiles)
 export("%>%")
@@ -70,7 +74,7 @@ export(knnarx_forecaster)
 export(layer)
 export(layer_naomit)
 export(layer_predict)
-export(layer_residual_quantile)
+export(layer_residual_quantiles)
 export(layer_threshold)
 export(nested_quantiles)
 export(new_default_epi_recipe_blueprint)

R/compat-purrr.R

@@ -2,7 +2,7 @@
 
 
 map <- function(.x, .f, ...) {
-  .f <- rlang::as_function(.f, env = global_env())
+  .f <- rlang::as_function(.f, env = rlang::global_env())
   lapply(.x, .f, ...)
 }
 walk <- function(.x, .f, ...) {
@@ -23,14 +23,14 @@ map_chr <- function(.x, .f, ...) {
   .rlang_purrr_map_mold(.x, .f, character(1), ...)
 }
 .rlang_purrr_map_mold <- function(.x, .f, .mold, ...) {
-  .f <- rlang::as_function(.f, env = global_env())
+  .f <- rlang::as_function(.f, env = rlang::global_env())
   out <- vapply(.x, .f, .mold, ..., USE.NAMES = FALSE)
   names(out) <- names(.x)
   out
 }
 
 map2 <- function(.x, .y, .f, ...) {
-  .f <- rlang::as_function(.f, env = global_env())
+  .f <- rlang::as_function(.f, env = rlang::global_env())
   out <- mapply(.f, .x, .y, MoreArgs = list(...), SIMPLIFY = FALSE)
   if (length(out) == length(.x)) {
     rlang::set_names(out, names(.x))

R/dist_quantiles.R

@@ -10,7 +10,8 @@ new_quantiles <- function(q = double(), tau = double()) {
     q <- q[o]
     tau <- tau[o]
   }
-  if (is.unsorted(q)) rlang::abort("q[order(tau)] produces unsorted quantiles.")
+  if (is.unsorted(q, na.rm = TRUE))
+    rlang::abort("q[order(tau)] produces unsorted quantiles.")
 
   new_rcrd(list(q = q, tau = tau),
            class = c("dist_quantiles", "dist_default"))
@@ -93,8 +94,7 @@ extrapolate_quantiles <- function(x, p, ...) {
 #' @export
 extrapolate_quantiles.distribution <- function(x, p, ...) {
   arg_is_probabilities(p)
-  p <- distributional:::arg_listable(p, .ptype = double())
-  dstn <- distributional:::dist_apply(x, extrapolate_quantiles, p = p, ...)
+  dstn <- lapply(vec_data(x), extrapolate_quantiles, p = p, ...)
   distributional:::wrap_dist(dstn)
 }
 

R/layer_residual_quantile.R renamed to R/layer_residual_quantiles.R

@@ -2,46 +2,47 @@
 #'
 #' @param frosting a `frosting` postprocessor
 #' @param ... Unused, include for consistency with other layers.
-#' @param probs numeric vector of probabilities with values in (0,1) referring to the desired quantile.
-#' @param symmetrize logical. If `TRUE` then interval will be symmetrical.
+#' @param probs 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 .flag a logical to determine if the layer is added. Passed on to
 #'   `add_layer()`. Default `TRUE`.
 #' @param id a random id string
 #'
-#' @return an updated `frosting` postprocessor with additional columns of the residual quantiles added to the prediction
+#' @return an updated `frosting` postprocessor with additional columns of the
+#'   residual quantiles added to the prediction
 #' @export
 #' @examples
-#'  jhu <- case_death_rate_subset %>%
+#' jhu <- case_death_rate_subset %>%
 #'   dplyr::filter(time_value > "2021-11-01", geo_value %in% c("ak", "ca", "ny"))
 #'
 #' r <- epi_recipe(jhu) %>%
 #'   step_epi_lag(death_rate, lag = c(0, 7, 14)) %>%
 #'   step_epi_ahead(death_rate, ahead = 7) %>%
-#'   recipes::step_naomit(recipes::all_predictors()) %>%
-#'   recipes::step_naomit(recipes::all_outcomes(), skip = TRUE)
+#'   step_epi_naomit()
 #'
 #' wf <- epi_workflow(r, parsnip::linear_reg()) %>%
 #'  parsnip::fit(jhu)
 #'
 #' latest <- get_test_data(recipe = r, x = jhu)
 #'
-#' f <- epipredict:::frosting() %>%
-#'      layer_predict() %>%
-#'      layer_residual_quantile(probs = c(0.0275, 0.975), symmetrize = FALSE) %>%
-#'      layer_naomit(.pred)
-#' wf1 <- wf %>% epipredict:::add_frosting(f)
+#' f <- frosting() %>%
+#'   layer_predict() %>%
+#'   layer_residual_quantiles(probs = c(0.0275, 0.975), symmetrize = FALSE) %>%
+#'   layer_naomit(.pred)
+#' wf1 <- wf %>% add_frosting(f)
 #'
 #' p <- predict(wf1, latest)
 #' p
-layer_residual_quantile <- function(frosting, ...,
+layer_residual_quantiles <- function(frosting, ...,
                                     probs = c(0.0275, 0.975),
                                     symmetrize = TRUE,
                                     .flag = TRUE,
-                                    id = rand_id("residual_quantile")) {
+                                    id = rand_id("residual_quantiles")) {
   rlang::check_dots_empty()
   add_layer(
     frosting,
-    layer_residual_quantile_new(
+    layer_residual_quantiles_new(
       probs = probs,
       symmetrize = symmetrize,
       id = id
@@ -50,12 +51,12 @@ layer_residual_quantile <- function(frosting, ...,
   )
 }
 
-layer_residual_quantile_new <- function(probs, symmetrize, id) {
-  layer("residual_quantile", probs = probs, symmetrize = symmetrize, id = id)
+layer_residual_quantiles_new <- function(probs, symmetrize, id) {
+  layer("residual_quantiles", probs = probs, symmetrize = symmetrize, id = id)
 }
 
 #' @export
-slather.layer_residual_quantile <-
+slather.layer_residual_quantiles <-
   function(object, components, the_fit, the_recipe, ...) {
     if (is.null(object$probs)) return(components)
 
@@ -64,8 +65,7 @@ slather.layer_residual_quantile <-
     q <- quantile(c(r, s * r), probs = object$probs, na.rm = TRUE)
 
     estimate <- components$predictions$.pred
-    interval <- data.frame(outer(estimate, q, "+"))
-    names(interval)<- probs_to_string(object$probs)
-    components$predictions <- dplyr::bind_cols(components$predictions,interval)
+    dstn <- dist_quantiles(map(estimate, "+", q), object$probs)
+    components$predictions$.quantiles <- dstn
     components
   }

R/layer_threshold_preds.R

@@ -67,11 +67,41 @@ layer_threshold_new <-
     layer("threshold", terms = terms, lower = lower, upper = upper, id = id)
   }
 
-snap <- function(x, lower, upper) {
+
+
+snap <- function(x, lower, upper, ...) {
+  UseMethod("snap")
+}
+
+#' @export
+snap.default <- function(x, lower, upper, ...) {
+  rlang::check_dots_empty()
   arg_is_scalar(lower, upper)
   pmin(pmax(x, lower), upper)
 }
 
+#' @export
+snap.distribution <- function(x, lower, upper, ...) {
+  rlang::check_dots_empty()
+  arg_is_scalar(lower, upper)
+  dstn <- lapply(vec_data(x), snap, lower = lower, upper = upper)
+  distributional:::wrap_dist(dstn)
+}
+
+#' @export
+snap.dist_default <- function(x, lower, upper, ...) {
+  rlang::check_dots_empty()
+  x
+}
+
+#' @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)
+}
+
 #' @export
 slather.layer_threshold <-
   function(object, components, the_fit, the_recipe, ...) {

man/layer_residual_quantile.Rd renamed to man/layer_residual_quantiles.Rd

@@ -19,14 +19,13 @@ test_that("outcome of the two methods are the same", {
     step_epi_lag(death_rate, lag = c(0, 7)) %>%
     step_epi_ahead(death_rate, ahead = 7) %>%
     step_epi_lag(case_rate, lag = c(7)) %>%
-    step_naomit(all_predictors()) %>%
-    step_naomit(all_outcomes())
+    step_epi_naomit()
 
   s <- parsnip::linear_reg()
   f <- frosting() %>%
     layer_predict() %>%
     layer_naomit(.pred) %>%
-    layer_residual_quantile()
+    layer_residual_quantiles()
 
   ef <- epi_workflow(r, s, f)
   ef2 <- epi_workflow(r, s) %>% add_frosting(f)

tests/testthat/test-epi_workflow.R

@@ -56,21 +56,20 @@ test_that("layer_predict is added by default if missing", {
   r <- epi_recipe(jhu) %>%
     step_epi_lag(death_rate, lag = c(0, 7, 14)) %>%
     step_epi_ahead(death_rate, ahead = 7) %>%
-    step_naomit(all_predictors()) %>%
-    step_naomit(all_outcomes(), skip = TRUE)
+    step_epi_naomit()
 
   wf <- epi_workflow(r, parsnip::linear_reg()) %>% fit(jhu)
 
   latest <- get_test_data(recipe = r, x = jhu)
 
   f1 <- frosting() %>%
     layer_naomit(.pred) %>%
-    layer_residual_quantile()
+    layer_residual_quantiles()
 
   f2 <- frosting() %>%
     layer_predict() %>%
     layer_naomit(.pred) %>%
-    layer_residual_quantile()
+    layer_residual_quantiles()
 
   wf1 <- wf %>% add_frosting(f1)
   wf2 <- wf %>% add_frosting(f2)

tests/testthat/test-frosting.R

@@ -4,8 +4,7 @@ jhu <- case_death_rate_subset %>%
 r <- epi_recipe(jhu) %>%
   step_epi_lag(death_rate, lag = c(0, 7, 14)) %>%
   step_epi_ahead(death_rate, ahead = 7) %>%
-  step_naomit(all_predictors()) %>%
-  step_naomit(all_outcomes(), skip = TRUE)
+  step_epi_naomit()
 
 wf <- epi_workflow(r, parsnip::linear_reg()) %>% fit(jhu)
 latest <- get_test_data(recipe = r, x = jhu)
@@ -15,13 +14,19 @@ test_that("Returns expected number or rows and columns", {
   f <- frosting() %>%
     layer_predict() %>%
     layer_naomit(.pred) %>%
-    layer_residual_quantile(probs = c(0.0275, 0.8, 0.95), symmetrize = FALSE)
+    layer_residual_quantiles(probs = c(0.0275, 0.8, 0.95), symmetrize = FALSE)
 
   wf1 <- wf %>% add_frosting(f)
 
   expect_silent(p <- predict(wf1, latest))
-  expect_equal(ncol(p), 6L)
+  expect_equal(ncol(p), 4L)
   expect_s3_class(p, "epi_df")
   expect_equal(nrow(p), 3L)
-  expect_named(p, c("geo_value", "time_value",".pred","q0.0275","q0.8","q0.95"))
+  expect_named(p, c("geo_value", "time_value",".pred",".quantiles"))
+
+  nested <- p %>% dplyr::mutate(.quantiles = nested_quantiles(.quantiles))
+  unnested <- nested %>% tidyr::unnest(.quantiles)
+
+  expect_equal(nrow(unnested), 9L)
+  expect_equal(unique(unnested$tau), c(.0275, .8, .95))
 })

tests/testthat/test-layer_residual_quantile.R renamed to tests/testthat/test-layer_residual_quantiles.R

@@ -3,8 +3,8 @@ jhu <- case_death_rate_subset %>%
 r <- epi_recipe(jhu) %>%
   step_epi_lag(death_rate, lag = c(0, 7, 14)) %>%
   step_epi_ahead(death_rate, ahead = 7) %>%
-  step_naomit(all_predictors()) %>%
-  step_naomit(all_outcomes(), skip = TRUE)
+  step_epi_naomit()
+
 wf <- epi_workflow(r, parsnip::linear_reg()) %>% fit(jhu)
 latest <- jhu %>%
   dplyr::filter(time_value >= max(time_value) - 14)
@@ -46,17 +46,21 @@ test_that("thresholds additional columns", {
 
   f <- frosting() %>%
     layer_predict() %>%
-    layer_residual_quantile(probs = c(.1, .9)) %>%
-    layer_threshold(.pred, dplyr::starts_with("q"), lower = 0.180, upper = 0.31) %>%
+    layer_residual_quantiles(probs = c(.1, .9)) %>%
+    layer_threshold(.pred, .quantiles, lower = 0.180, upper = 0.31) %>%
     layer_naomit(.pred)
+
   wf2 <- wf %>% add_frosting(f)
 
   expect_silent(p <- predict(wf2, latest))
-  expect_equal(ncol(p), 5L)
+  expect_equal(ncol(p), 4L)
   expect_s3_class(p, "epi_df")
   expect_equal(nrow(p), 3L)
   expect_equal(round(p$.pred, digits = 3), c(0.180, 0.180, 0.310))
-  expect_equal(round(p$q0.1, digits = 3), c(0.180, 0.180, 0.310))
-  expect_equal(round(p$q0.9, digits = 3), c(0.310, 0.180, 0.310))
-  expect_named(p, c("geo_value", "time_value", ".pred", "q0.1", "q0.9"))
+  expect_named(p, c("geo_value", "time_value", ".pred", ".quantiles"))
+  p <- p %>%
+    dplyr::mutate(.quantiles = nested_quantiles(.quantiles)) %>%
+    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))
 })