Merge pull request #429 from cmu-delphi/distn_quantileNA

fix for dist_quantile all `NA` s

Author: dsweber2

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: epipredict
 Title: Basic epidemiology forecasting methods
-Version: 0.1.5
+Version: 0.1.6
 Authors@R: c(
     person("Daniel J.", "McDonald", , "[email protected]", role = c("aut", "cre")),
     person("Ryan", "Tibshirani", , "[email protected]", role = "aut"),

NEWS.md

@@ -23,6 +23,7 @@ Pre-1.0.0 numbering scheme: 0.x will indicate releases, while 0.0.x will indicat
 ## Bug fixes
 - Shifting no columns results in no error for either `step_epi_ahead` and `step_epi_lag`
 - Quantiles produced by `grf` were sometimes out of order.
+- dist_quantiles can have all `NA` values without causing unrelated errors
 
 # epipredict 0.1
 

R/dist_quantiles.R

@@ -128,10 +128,12 @@ is_dist_quantiles <- function(x) {
 median.dist_quantiles <- function(x, na.rm = FALSE, ..., middle = c("cubic", "linear")) {
   quantile_levels <- field(x, "quantile_levels")
   values <- field(x, "values")
+  # we have exactly that quantile
   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) {
+  # if there's only 1 quantile_level (and it isn't 0.5), or the smallest quantile is larger than 0.5 or the largest smaller than 0.5, or if every value is NA, return NA
+  if (length(quantile_levels) < 2 || min(quantile_levels) > 0.5 || max(quantile_levels) < 0.5 || all(is.na(values))) {
     return(NA)
   }
   if (length(quantile_levels) < 3 || min(quantile_levels) > .25 || max(quantile_levels) < .75) {
@@ -161,6 +163,9 @@ quantile_extrapolate <- function(x, tau_out, middle) {
   tau <- field(x, "quantile_levels")
   qvals <- field(x, "values")
   nas <- is.na(qvals)
+  if (all(nas)) {
+    return(rep(NA, times = length(tau_out)))
+  }
   qvals_out <- rep(NA, length(tau_out))
   qvals <- qvals[!nas]
   tau <- tau[!nas]

man/step_adjust_latency.Rd

@@ -110,3 +110,14 @@ test_that("arithmetic works on quantiles", {
   expect_snapshot(error = TRUE, sum(dstn))
   expect_snapshot(error = TRUE, suppressWarnings(dstn + distributional::dist_normal()))
 })
+
+test_that("quantile.dist_quantile works for NA vectors", {
+  distn <- dist_quantiles(
+    list(c(NA, NA)),
+    list(1:2 / 3)
+  )
+  expect_true(is.na(quantile(distn, p = 0.5)))
+  expect_true(is.na(median(distn)))
+  expect_true(is.na(mean(distn)))
+  expect_equal(format(distn), "quantiles(NA)[2]")
+})

tests/testthat/test-dist_quantiles.R

@@ -385,8 +385,7 @@ canada_archive_faux <- epix_as_of(canada_archive, canada_archive$versions_end) %
 smooth_cases <- function(epi_df) {
   epi_df %>%
     group_by(geo_value) %>%
-    epi_slide_mean("case_rate", .window_size = 7, na.rm = TRUE) %>%
-    rename(cr_7dav = slide_value_case_rate)
+    epi_slide_mean("case_rate", .window_size = 7, na.rm = TRUE, .suffix = "_{.n}dav")
 }
 forecast_dates <- seq.Date(
   from = min(canada_archive$DT$version),
@@ -401,8 +400,8 @@ canada_forecasts <- bind_rows(
     ~ forecast_k_week_ahead(
       canada_archive_faux,
       ahead = .x,
-      outcome = "cr_7dav",
-      predictors = "cr_7dav",
+      outcome = "case_rate_7dav",
+      predictors = "case_rate_7dav",
       forecast_dates = forecast_dates,
       process_data = smooth_cases
     ) %>% mutate(version_aware = FALSE)
@@ -412,8 +411,8 @@ canada_forecasts <- bind_rows(
     ~ forecast_k_week_ahead(
       canada_archive,
       ahead = .x,
-      outcome = "cr_7dav",
-      predictors = "cr_7dav",
+      outcome = "case_rate_7dav",
+      predictors = "case_rate_7dav",
       forecast_dates = forecast_dates,
       process_data = smooth_cases
     ) %>% mutate(version_aware = TRUE)
@@ -435,15 +434,15 @@ case_rate_data <- bind_rows(
     ~ canada_archive %>%
       epix_as_of(.x) %>%
       smooth_cases() %>%
-      mutate(case_rate = cr_7dav, version = .x)
+      mutate(case_rate = case_rate_7dav, version = .x)
   ) %>%
     bind_rows() %>%
     mutate(version_aware = TRUE),
   # Latest data for the version-unaware forecasts
   canada_archive %>%
     epix_as_of(doctor_visits$versions_end) %>%
     smooth_cases() %>%
-    mutate(case_rate = cr_7dav, version_aware = FALSE)
+    mutate(case_rate = case_rate_7dav, version_aware = FALSE)
 ) %>%
   filter(geo_value == geo_choose)