Merge branch 'v0.0.5' into smooth-quant-reg

Author: dajmcdon

NAMESPACE

@@ -140,6 +140,7 @@ export(layer_unnest)
 export(nested_quantiles)
 export(new_default_epi_recipe_blueprint)
 export(new_epi_recipe_blueprint)
+export(pivot_quantiles)
 export(prep)
 export(quantile_reg)
 export(remove_frosting)

R/dist_quantiles.R

@@ -107,7 +107,9 @@ extrapolate_quantiles.dist_quantiles <- function(x, p, ...) {
   new_quantiles(q = c(qvals, q), tau = c(tau, p))
 }
 
-
+is_dist_quantiles <- function(x) {
+  is_distribution(x) && all(stats::family(x) == "quantiles")
+}
 
 
 #' Turn a a vector of quantile distributions into a list-col
@@ -124,8 +126,7 @@ extrapolate_quantiles.dist_quantiles <- function(x, p, ...) {
 #' edf_nested <- edf %>% dplyr::mutate(q = nested_quantiles(q))
 #' edf_nested %>% tidyr::unnest(q)
 nested_quantiles <- function(x) {
-  stopifnot(is_distribution(x),
-            all(stats::family(x) == "quantiles"))
+  stopifnot(is_dist_quantiles(x))
   distributional:::dist_apply(x, .f = function(z) {
     tibble::as_tibble(vec_data(z)) %>%
       dplyr::mutate(dplyr::across(tidyselect::everything(), as.double)) %>%
@@ -134,6 +135,68 @@ nested_quantiles <- function(x) {
 }
 
 
+#' Pivot columns containing `dist_quantile` wider
+#'
+#' Any selected columns that contain `dist_quantiles` will be "widened" with
+#' the "taus" (quantile) serving as names and the values in the data frame.
+#' When pivoting multiple columns, the original column name will be used as
+#' a prefix.
+#'
+#' @param .data A data frame, or a data frame extension such as a tibble or
+#'   epi_df.
+#' @param ... <[`tidy-select`][dplyr::dplyr_tidy_select]> One or more unquoted
+#'   expressions separated by commas. Variable names can be used as if they
+#'   were positions in the data frame, so expressions like `x:y` can
+#'   be used to select a range of variables. Any selected columns should
+#'
+#' @return An object of the same class as `.data`
+#' @export
+#'
+#' @examples
+#' d1 <- c(dist_quantiles(1:3, 1:3 / 4), dist_quantiles(2:4, 1:3 / 4))
+#' d2 <- c(dist_quantiles(2:4, 2:4 / 5), dist_quantiles(3:5, 2:4 / 5))
+#' tib <- tibble::tibble(g = c("a", "b"), d1 = d1, d2 = d2)
+#'
+#' pivot_quantiles(tib, c("d1", "d2"))
+#' pivot_quantiles(tib, tidyselect::starts_with("d"))
+#' pivot_quantiles(tib, d2)
+pivot_quantiles <- function(.data, ...) {
+  expr <- rlang::expr(c(...))
+  cols <- names(tidyselect::eval_select(expr, .data))
+  dqs <- map_lgl(cols, ~ is_dist_quantiles(.data[[.x]]))
+  if (!all(dqs)) {
+    nms <- cols[!dqs]
+    cli::cli_abort(
+      "Variables(s) {.var {nms}} are not `dist_quantiles`. Cannot pivot them."
+    )
+  }
+  .data <- .data %>%
+    dplyr::mutate(dplyr::across(tidyselect::all_of(cols), nested_quantiles))
+  checks <- map_lgl(cols, ~ diff(range(vctrs::list_sizes(.data[[.x]]))) == 0L)
+  if (!all(checks)) {
+    nms <- cols[!checks]
+    cli::cli_abort(
+      c("Quantiles must be the same length and have the same set of taus.",
+        i = "Check failed for variables(s) {.var {nms}}."))
+  }
+  if (length(cols) > 1L) {
+    for (col in cols) {
+      .data <- .data %>%
+        tidyr::unnest(tidyselect::all_of(col)) %>%
+        tidyr::pivot_wider(
+          names_from = "tau", values_from = "q",
+          names_prefix = paste0(col, "_")
+        )
+    }
+  } else {
+    .data <- .data %>%
+      tidyr::unnest(tidyselect::all_of(cols)) %>%
+      tidyr::pivot_wider(names_from = "tau", values_from = "q")
+  }
+  .data
+}
+
+
 
 
 #' @export

R/flatline.R

@@ -73,11 +73,31 @@ predict.flatline <- function(object, newdata, ...) {
   object <- object$.pred
   metadata <- names(object)[names(object) != ".pred"]
   ek <- names(newdata)
-  if (! all(metadata %in% ek)) {
+  if (!all(metadata %in% ek)) {
     cli_stop("`newdata` has different metadata than was used",
              "to fit the flatline forecaster")
   }
 
   dplyr::left_join(newdata, object, by = metadata) %>%
     dplyr::pull(.pred)
 }
+
+#' @export
+print.flatline <- function(x, ...) {
+  keys <- colnames(x$.pred)
+  keys <- paste(keys[!(keys %in% ".pred")], collapse = ", ")
+  nloc <- nrow(x$.pred)
+  nres <- nrow(x$residuals)
+  pmsg <- glue::glue(
+    "Predictions produced by {keys} resulting in {nloc} total forecasts."
+  )
+  rmsg <- glue::glue(
+    "A total of {nres} residuals are available from the training set."
+  )
+  cat("Flatline forecaster\n")
+  cat("\n")
+  cat(pmsg)
+  cat("\n")
+  cat(rmsg)
+  cat("\n\n")
+}

R/flatline_forecaster.R

@@ -37,14 +37,14 @@ flatline_forecaster <- function(
     cli_stop("args_list was not created using `flatline_args_list().")
   }
   keys <- epi_keys(epi_data)
-  ek <- keys[-1]
+  ek <- kill_time_value(keys)
   outcome <- rlang::sym(outcome)
 
 
   r <- epi_recipe(epi_data) %>%
     step_epi_ahead(!!outcome, ahead = args_list$ahead, skip = TRUE) %>%
     recipes::update_role(!!outcome, new_role = "predictor") %>%
-    recipes::add_role(dplyr::all_of(keys), new_role = "predictor") %>%
+    recipes::add_role(tidyselect::all_of(keys), new_role = "predictor") %>%
     step_training_window(n_recent = args_list$n_training)
 
   latest <- get_test_data(epi_recipe(epi_data), epi_data)
@@ -65,7 +65,6 @@ flatline_forecaster <- function(
   eng <- parsnip::linear_reg() %>% parsnip::set_engine("flatline")
 
   wf <- epi_workflow(r, eng, f)
-
   wf <- generics::fit(wf, epi_data)
   preds <- suppressWarnings(predict(wf, new_data = latest)) %>%
     tibble::as_tibble() %>%

R/layer_residual_quantiles.R

@@ -75,29 +75,43 @@ slather.layer_residual_quantiles <-
   function(object, components, the_fit, the_recipe, ...) {
     if (is.null(object$probs)) return(components)
 
-
     s <- ifelse(object$symmetrize, -1, NA)
-    r <- dplyr::bind_cols(
-      r = grab_residuals(the_fit, components),
-      geo_value = components$mold$extras$roles$geo_value,
-      components$mold$extras$roles$key)
+    r <- grab_residuals(the_fit, components)
 
     ## Handle any grouping requests
     if (length(object$by_key) > 0L) {
-      common <- intersect(object$by_key, names(r))
-      excess <- setdiff(object$by_key, names(r))
+      key_cols <- dplyr::bind_cols(
+        geo_value = components$mold$extras$roles$geo_value,
+        components$mold$extras$roles$key
+      )
+      common <- intersect(object$by_key, names(key_cols))
+      excess <- setdiff(object$by_key, names(key_cols))
       if (length(excess) > 0L) {
-        cli_warn("Requested residual grouping key(s) {excess} unavailable ",
-            "in the original data. Grouping by the remainder {common}.")
-
+        rlang::warn(
+          "Requested residual grouping key(s) {excess} are unavailable ",
+          "in the original data. Grouping by the remainder: {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(
+            "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))
       }
-      if (length(common) > 0L)
-        r <- r %>% dplyr::group_by(!!!rlang::syms(common))
     }
 
     r <- r %>%
-      dplyr::summarise(
-        q = list(quantile(c(r, s * r), probs = object$probs, na.rm = TRUE))
+      dplyr::summarize(
+        q = list(quantile(
+          c(.resid, s * .resid),
+          probs = object$probs, na.rm = TRUE
+        ))
       )
 
     estimate <- components$predictions$.pred
@@ -112,13 +126,40 @@ slather.layer_residual_quantiles <-
 grab_residuals <- function(the_fit, components) {
   if (the_fit$spec$mode != "regression")
     rlang::abort("For meaningful residuals, the predictor should be a regression model.")
-  r_generic <- attr(utils::methods(class = class(the_fit)[1]), "info")$generic
-  if ("residuals" %in% r_generic) {
-    r <- residuals(the_fit)
-  } else {
-    yhat <- predict(the_fit, new_data = components$mold$predictors)
-    r <- c(components$mold$outcomes - yhat)[[1]]
+  r_generic <- attr(utils::methods(class = class(the_fit$fit)[1]), "info")$generic
+  if ("residuals" %in% r_generic) { # Try to use the available method.
+    cl <- class(the_fit$fit)[1]
+    r <- residuals(the_fit$fit)
+    if (inherits(r, "data.frame")) {
+      if (".resid" %in% names(r)) { # success
+        return(r)
+      } else { # failure
+        rlang::warn(c(
+          "The `residuals()` method for objects of class {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.",
+          i = "This may result in unexpected behaviour."
+        ))
+      }
+    } else if (is.vector(drop(r))) { # also success
+      return(tibble(.resid = drop(r)))
+    } else { # failure
+      rlang::warn(c(
+        "The `residuals()` method for objects of class {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",
+        i = "difference between predictions and observations.",
+        i = "This may result in unexpected behaviour."
+      ))
+    }
   }
+  # The method failed for one reason or another and a warning was issued
+  # Or there was no method available.
+  yhat <- predict(the_fit, new_data = components$mold$predictors)
+  r <- c(components$mold$outcomes - yhat)[[1]] # this will be a vector
+  r <- tibble(.resid = r)
   r
 }