dev/main comparison

DESCRIPTION

@@ -8,15 +8,13 @@ Description: Provide real-time revision forecasts.
 License: MIT + file LICENSE
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 Suggests: 
     testthat (>= 3.0.0)
 Config/testthat/edition: 3
-Imports: 
+Imports:
     arrow,
-    covidcast,
     dplyr,
-    evalcast,
     english,
     jsonlite,
     lubridate,

NAMESPACE

@@ -21,6 +21,7 @@ export(WEEK_ISSUES)
 export(Y7DAV)
 export(YITL)
 export(add_7davs)
+export(aux_feature_names)
 export(add_dayofweek)
 export(add_lagged_terms)
 export(add_log_transformed)
@@ -68,7 +69,6 @@ importFrom(dplyr,slice_max)
 importFrom(dplyr,starts_with)
 importFrom(dplyr,ungroup)
 importFrom(english,english)
-importFrom(evalcast,weighted_interval_score)
 importFrom(jsonlite,read_json)
 importFrom(lubridate,days_in_month)
 importFrom(lubridate,make_date)

R/feature_engineering.R

@@ -188,8 +188,13 @@ add_lagged_terms <- function(df, value_col, refd_col, lag_col, lagged_term_list=
 #' @importFrom dplyr rename
 #' @export
 add_targets <- function(df, value_col, refd_col, lag_col, ref_lag, temporal_resol) {
-  # Add target
-  target_df <- df[df[[lag_col]]==ref_lag, c(refd_col, "report_date", value_col, "value_7dav")]
+  available_lags <- sort(unique(df[[lag_col]]))
+  effective_ref_lag <- min(available_lags[available_lags >= ref_lag])
+  if (is.infinite(effective_ref_lag))
+    stop(sprintf("ref_lag %d exceeds all available lags (max %d)", ref_lag, max(available_lags)))
+  if (effective_ref_lag != ref_lag)
+    message(sprintf("ref_lag %d not available; using next lag %d", ref_lag, effective_ref_lag))
+  target_df <- df[df[[lag_col]] == effective_ref_lag, c(refd_col, "report_date", value_col, "value_7dav")]
   # Rename columns for clarity
   target_df <- target_df %>%
     dplyr::rename(
@@ -268,6 +273,74 @@ add_params_for_dates <- function(df, refd_col, lag_col, temporal_resol="daily")
   return (as.data.frame(df))
 }
 
+#' Process an auxiliary reporting triangle into prefixed feature columns
+#'
+#' Runs a single auxiliary data frame (one geo, one signal) through the
+#' fill → 7-day-average → lagged-terms → log-transform pipeline and
+#' renames all value-derived columns with a `{name}_` prefix so they can
+#' be safely joined to the primary preprocessed data frame.
+#'
+#' @param df Data frame with columns `reference_date`, `report_date`, `lag`,
+#'   and `value` (the signal values).
+#' @param name Character scalar used as the column prefix (e.g. `"nssp"`).
+#' @param lagged_term_list Numeric vector of lag values (same as used for the
+#'   primary signal).
+#' @param temporal_resol `"daily"` or `"weekly"`.
+#' @param smoothed Logical; if `FALSE` and `temporal_resol == "daily"`, a 7-day
+#'   moving average is computed. Otherwise `value_7dav` is set equal to
+#'   `value_raw`.
+#'
+#' @return Data frame with columns `reference_date`, `report_date`, `lag`, and
+#'   all value/log columns prefixed with `{name}_`.  Useful predictors are
+#'   `{name}_log_value_7dav_lag{N}` and `{name}_log_delta_value_7dav_lag{N}`;
+#'   see [aux_feature_names()].
+process_aux_triangle <- function(df, name, lagged_term_list, temporal_resol, smoothed,
+                                max_report_override = NULL) {
+  filled_df <- fill_missing_updates(df, "value", "reference_date", "lag", temporal_resol,
+                                    max_report_override = max_report_override)
+  if (nrow(filled_df) == 0) {
+    expected_cols <- c("reference_date", "report_date", "lag",
+                       aux_feature_names(name, lagged_term_list))
+    return(setNames(data.frame(matrix(ncol = length(expected_cols), nrow = 0)),
+                    expected_cols))
+  }
+  if (!smoothed && temporal_resol == "daily") {
+    filled_df <- add_7davs(filled_df, "value_raw", "reference_date", "lag")
+  } else {
+    filled_df$value_7dav <- filled_df$value_raw
+  }
+  filled_df <- add_lagged_terms(
+    filled_df, "value_7dav", "reference_date", "lag", lagged_term_list, temporal_resol
+  )
+  filled_df <- add_log_transformed(filled_df, lagged_term_list)
+
+  value_cols <- grep("^(value_|log_)", colnames(filled_df), value = TRUE)
+  colnames(filled_df)[colnames(filled_df) %in% value_cols] <- paste0(name, "_", value_cols)
+
+  filled_df[, c("reference_date", "report_date", "lag", paste0(name, "_", value_cols))]
+}
+
+
+#' Return the feature column names produced by an auxiliary triangle
+#'
+#' Gives the column names that [process_aux_triangle()] adds for a given
+#' auxiliary signal, matching the log-value and log-delta features used by
+#' the primary signal in [create_params_list()].
+#'
+#' @param name Character scalar; the aux triangle name (must match what was
+#'   passed to [data_preprocessing()]).
+#' @param lagged_term_list Numeric vector of lag values.
+#'
+#' @return Character vector of feature column names.
+#' @export
+aux_feature_names <- function(name, lagged_term_list) {
+  c(
+    paste0(name, "_log_value_7dav_lag", lagged_term_list),
+    paste0(name, "_log_delta_value_7dav_lag", lagged_term_list)
+  )
+}
+
+
 #' Data Preprocessing Function
 #'
 #' This function processes input data by handling missing values, computing lagged terms,
@@ -287,14 +360,22 @@ add_params_for_dates <- function(df, refd_col, lag_col, temporal_resol="daily")
 #' @param value_type Character indicating the type of values ('count' or 'fraction').
 #' @param temporal_resol Character specifying temporal resolution ('daily' or 'weekly').
 #' @param smoothed Logical indicating whether smoothing should be applied.
-#'
-#' @importFrom dplyr full_join distinct
+#' @param aux_triangles Named list of auxiliary reporting-triangle data frames.
+#'   Each element must have columns `reference_date`, `report_date`, `lag`, and
+#'   `value` (already filtered to the same single geo as `df`).  Each is run
+#'   through the same fill → lag → log pipeline as the primary signal and
+#'   joined to the result; columns are prefixed with the list element name.
+#'   Use [aux_feature_names()] to obtain the resulting predictor column names
+#'   for [create_params_list()].
+#'
+#' @importFrom dplyr full_join left_join distinct
 #' @importFrom english english
 #'
 #' @export
 data_preprocessing <- function(df, value_col, refd_col, lag_col, ref_lag,
                                suffixes=c(""), lagged_term_list = NULL, value_type="count",
-                               temporal_resol="daily", smoothed=FALSE) {
+                               temporal_resol="daily", smoothed=FALSE,
+                               aux_triangles = NULL) {
   if (value_type == "count") {
     if (length(value_col) > 1) warning("Multiple value column names provided; only the first one will be used.")
     if (length(unique(suffixes)) > 1) warning("Multiple suffixes provided; only the first one will be used.")
@@ -369,6 +450,25 @@ data_preprocessing <- function(df, value_col, refd_col, lag_col, ref_lag,
   merged_df$inv_log_lag <- 1/(merged_df$lag + 1)
   merged_df <- add_params_for_dates(merged_df, "reference_date", "lag", temporal_resol)
 
+  if (!is.null(aux_triangles)) {
+    primary_max_report <- max(merged_df$report_date)
+    for (nm in names(aux_triangles)) {
+      aux_processed <- process_aux_triangle(
+        aux_triangles[[nm]], nm, lagged_term_list, temporal_resol, smoothed,
+        max_report_override = primary_max_report
+      )
+      if (nrow(aux_processed) == 0L) {
+        merged_df <- merged_df[0L, ]
+        break
+      }
+      merged_df <- dplyr::left_join(
+        merged_df, aux_processed,
+        by = c("reference_date", "report_date", "lag")
+      )
+      merged_df <- merged_df[!is.na(merged_df[[paste0(nm, "_value_raw")]]), ]
+    }
+  }
+
   merged_df <- merged_df %>%
     filter(.data$lag < ref_lag)
 

R/forecast.R

@@ -41,6 +41,7 @@ revision_forecast <- function(train_data, test_data, taus,
                               smoothed_target=TRUE,
                               lagged_term_list=NULL,
                               params_list=NULL,
+                              extra_params=NULL,
                               temporal_resol="daily",
                               lambda = 0.1, gamma = 0.1,
                               lp_solver=LP_SOLVER, test_lag_group="",
@@ -77,7 +78,7 @@ revision_forecast <- function(train_data, test_data, taus,
   }
 
   if (is.null(params_list)) {
-    params_list <- create_params_list(train_data, lagged_term_list, temporal_resol)
+    params_list <- create_params_list(train_data, lagged_term_list, temporal_resol, extra_params)
   }
 
   if (smoothed_target) {
@@ -317,6 +318,7 @@ DelphiRF <- function(df, testing_start_date, taus=TAUS,
                      smoothed_target=TRUE,
                      lagged_term_list=NULL,
                      params_list=NULL,
+                     extra_params=NULL,
                      lambda=LAMBDA, gamma=GAMMA, lag_pad=LAG_PAD,
                      temporal_resol="daily",
                      lp_solver=LP_SOLVER,
@@ -388,7 +390,7 @@ DelphiRF <- function(df, testing_start_date, taus=TAUS,
 
     results <- revision_forecast(train_data, test_data, taus,
                                  smoothed_target, lagged_term_list,
-                                 params_list, temporal_resol,
+                                 params_list, extra_params, temporal_resol,
                                  l, g, lp_solver, test_lag_group,
                                  geo, value_type, model_save_dir,
                                  indicator, signal, geo_level,

R/model.R

@@ -142,17 +142,26 @@ get_prediction <- function(test_data, taus, covariates, response, obj,
   return (as.data.frame(test_data))
 }
 
+#' Weighted interval score for a single observation
+#'
+#' Inlined from the evalcast package.
+#'
+#' @param taus Numeric vector of quantile levels.
+#' @param residuals Numeric vector of (quantile_prediction - actual) values.
+#' @param point_pred Unused; kept for interface compatibility.
+#' @keywords internal
+weighted_interval_score <- function(taus, residuals, point_pred) {
+  alpha <- 2 * pmin(taus, 1 - taus)
+  mean(alpha * (abs(residuals) + (residuals) * (2 * (taus >= 0.5) - 1)))
+}
+
 #' Evaluation of the test results based on WIS score
-#' The WIS score calculation is based on the weighted_interval_score function
-#' from the `evalcast` package from Delphi
 #'
 #' @param test_data dataframe with a column containing the prediction results of
 #'    each requested quantile. Each row represents an update with certain
 #'    (reference_date, report_date, location) combination.
 #' @template taus-template
 #'
-#' @importFrom evalcast weighted_interval_score
-#'
 #' @export
 evaluate <- function(test_data, taus, response) {
   n_row <- nrow(test_data)
@@ -332,7 +341,7 @@ generate_filename <- function(indicator, signal,
 #'
 #' @importFrom dplyr mutate select
 #'
-create_params_list <- function(train_data, lagged_term_list, temporal_resol) {
+create_params_list <- function(train_data, lagged_term_list, temporal_resol, extra_params = NULL) {
   params_list <- c(
     WEEK_ISSUES[1],
     Y7DAV,
@@ -350,9 +359,6 @@ create_params_list <- function(train_data, lagged_term_list, temporal_resol) {
     paste0(dayofweek, "_issue")
   )
 
-  if (temporal_resol == "daily"){
-    return (c(params_list, extra_params_for_daily))
-  } else {
-    return(params_list)
-  }
+  base_params <- if (temporal_resol == "daily") c(params_list, extra_params_for_daily) else params_list
+  if (!is.null(extra_params)) c(base_params, extra_params) else base_params
 }

R/preprocessing.R

@@ -56,7 +56,8 @@ fill_rows <- function(df, refd_col, lag_col, min_refd, max_refd, ref_lag) {
 #' @importFrom tidyr fill pivot_wider pivot_longer replace_na expand_grid
 #' @importFrom dplyr %>% select left_join mutate arrange distinct filter everything
 #' @export
-fill_missing_updates <- function(df, value_col, refd_col, lag_col, temporal_resol="daily") {
+fill_missing_updates <- function(df, value_col, refd_col, lag_col, temporal_resol = "daily",
+                              max_report_override = NULL) {
   df <- df %>% distinct()  # Remove duplicates if any
 
   if (nrow(df) == 0) {
@@ -88,6 +89,14 @@ fill_missing_updates <- function(df, value_col, refd_col, lag_col, temporal_reso
     stop("Invalid temporal_resol. Choose either 'daily' or 'weekly'.")
   }
 
+  if (!is.null(max_report_override)) {
+    max_report_override <- as.Date(max_report_override)
+    if (max_report_override > max(all_report_dates)) {
+      extra <- seq(max(all_report_dates) + gap, max_report_override, by = gap)
+      all_report_dates <- c(all_report_dates, extra)
+    }
+  }
+
   # Create a complete grid of all combinations of reference_date and report_date
   complete_grid <- tidyr::expand_grid(
     !!refd_col := all_reference_dates,

R/utils.R

@@ -200,19 +200,13 @@ training_days_check <- function(report_date, training_days) {
 
 #' Subset list of counties to those included in the 200 most populous in the US
 #'
+#' Requires the covidcast package, which is not installed by default.
+#'
 #' @importFrom dplyr select %>% arrange desc pull
 #' @importFrom rlang .data
 #' @importFrom utils head
 get_populous_counties <- function() {
-  return(
-    covidcast::county_census %>%
-      dplyr::select(pop = .data$POPESTIMATE2019, fips = .data$FIPS) %>%
-      # Drop megacounties (states)
-      filter(!endsWith(.data$fips, "000")) %>%
-      arrange(desc(.data$pop)) %>%
-      pull(.data$fips) %>%
-      head(n=200)
-  )
+  stop("get_populous_counties() requires the covidcast package. Install it with renv::install(\"cmu-delphi/covidcast/R-packages/covidcast\").")
 }
 
 #' Write a message to the console with the current time