From f98cbafdca1667404f40ac5bfb08f24e6e6af87a Mon Sep 17 00:00:00 2001 From: dsweber2 Date: Thu, 11 Jun 2026 14:04:51 -0500 Subject: [PATCH 1/5] dropping stale packages --- DESCRIPTION | 4 +--- NAMESPACE | 1 - R/model.R | 17 +++++++++++++---- R/utils.R | 12 +++--------- 4 files changed, 17 insertions(+), 17 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index db43f7e..73d6981 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -12,11 +12,9 @@ RoxygenNote: 7.3.2 Suggests: testthat (>= 3.0.0) Config/testthat/edition: 3 -Imports: +Imports: arrow, - covidcast, dplyr, - evalcast, english, jsonlite, lubridate, diff --git a/NAMESPACE b/NAMESPACE index 0c0d5a1..aa5d5b7 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -68,7 +68,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) diff --git a/R/model.R b/R/model.R index 1af536f..5d0e3c7 100644 --- a/R/model.R +++ b/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) diff --git a/R/utils.R b/R/utils.R index 57f88f8..7103110 100644 --- a/R/utils.R +++ b/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 From c9068c9b3b0ef9742dbd76341eb740d398707155 Mon Sep 17 00:00:00 2001 From: dsweber2 Date: Thu, 9 Jul 2026 12:57:40 -0500 Subject: [PATCH 2/5] ref_lag not pinned to 7 days --- R/feature_engineering.R | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/R/feature_engineering.R b/R/feature_engineering.R index c253065..3d45ec8 100644 --- a/R/feature_engineering.R +++ b/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( From b6dfc4dc616abee11f0bc4a0170b097ffd19c99b Mon Sep 17 00:00:00 2001 From: dsweber2 Date: Mon, 27 Jul 2026 16:58:13 -0500 Subject: [PATCH 3/5] initial exogenous --- R/feature_engineering.R | 86 +++++++++++++++++++++++++++++++++++++++-- R/forecast.R | 6 ++- R/model.R | 9 ++--- 3 files changed, 90 insertions(+), 11 deletions(-) diff --git a/R/feature_engineering.R b/R/feature_engineering.R index 3d45ec8..445c1fa 100644 --- a/R/feature_engineering.R +++ b/R/feature_engineering.R @@ -273,6 +273,66 @@ 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) { + filled_df <- fill_missing_updates(df, "value", "reference_date", "lag", temporal_resol) + 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, @@ -292,14 +352,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.") @@ -374,6 +442,18 @@ 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)) { + for (nm in names(aux_triangles)) { + aux_processed <- process_aux_triangle( + aux_triangles[[nm]], nm, lagged_term_list, temporal_resol, smoothed + ) + merged_df <- dplyr::left_join( + merged_df, aux_processed, + by = c("reference_date", "report_date", "lag") + ) + } + } + merged_df <- merged_df %>% filter(.data$lag < ref_lag) diff --git a/R/forecast.R b/R/forecast.R index 5ab1b21..f54127a 100644 --- a/R/forecast.R +++ b/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, diff --git a/R/model.R b/R/model.R index 5d0e3c7..70006e3 100644 --- a/R/model.R +++ b/R/model.R @@ -341,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, @@ -359,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 } From b900add70d7227395a81c6236b3ba0dc08a03c94 Mon Sep 17 00:00:00 2001 From: dsweber2 Date: Tue, 28 Jul 2026 10:47:48 -0500 Subject: [PATCH 4/5] deal with timing issues --- DESCRIPTION | 2 +- NAMESPACE | 1 + R/feature_engineering.R | 11 +++-- R/preprocessing.R | 11 ++++- tests/testthat/test-feature_engineering.R | 51 +++++++++++++++++++++++ 5 files changed, 71 insertions(+), 5 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 73d6981..f929e62 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -8,7 +8,7 @@ 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 diff --git a/NAMESPACE b/NAMESPACE index aa5d5b7..310924f 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/feature_engineering.R b/R/feature_engineering.R index 445c1fa..03710e6 100644 --- a/R/feature_engineering.R +++ b/R/feature_engineering.R @@ -294,8 +294,10 @@ add_params_for_dates <- function(df, refd_col, lag_col, temporal_resol="daily") #' 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) { - filled_df <- fill_missing_updates(df, "value", "reference_date", "lag", temporal_resol) +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 (!smoothed && temporal_resol == "daily") { filled_df <- add_7davs(filled_df, "value_raw", "reference_date", "lag") } else { @@ -443,14 +445,17 @@ data_preprocessing <- function(df, value_col, refd_col, lag_col, ref_lag, 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 + aux_triangles[[nm]], nm, lagged_term_list, temporal_resol, smoothed, + max_report_override = primary_max_report ) 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")]]), ] } } diff --git a/R/preprocessing.R b/R/preprocessing.R index 4371cec..20a10e2 100644 --- a/R/preprocessing.R +++ b/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, diff --git a/tests/testthat/test-feature_engineering.R b/tests/testthat/test-feature_engineering.R index c4a577f..13b81c4 100644 --- a/tests/testthat/test-feature_engineering.R +++ b/tests/testthat/test-feature_engineering.R @@ -341,3 +341,54 @@ test_that("Testing add weighted related features", { expect_true(all(expected_columns %in% colnames(result))) }) + + +make_daily_tri <- function(ref_dates, lags) { + do.call(rbind, lapply(ref_dates, function(rd) { + data.frame( + reference_date = rd, + lag = lags, + value = as.numeric(as.Date(rd) - as.Date("2022-12-31")) + ) + })) +} + +test_that("data_preprocessing drops rows where aux has no coverage (earlier min)", { + ref_all <- seq(as.Date("2023-01-01"), as.Date("2023-01-15"), by = "day") + ref_late <- seq(as.Date("2023-01-08"), as.Date("2023-01-15"), by = "day") + + primary <- make_daily_tri(ref_all, 1:3) + aux <- make_daily_tri(ref_late, 1:3) + + result <- data_preprocessing( + primary, + value_col = "value", refd_col = "reference_date", lag_col = "lag", + ref_lag = 5L, temporal_resol = "daily", + aux_triangles = list(beds = aux) + ) + + expect_true(all(result$reference_date >= as.Date("2023-01-08"))) + expect_true("beds_value_raw" %in% colnames(result)) + expect_false(anyNA(result[["beds_value_raw"]])) +}) + +test_that("data_preprocessing forward-fills aux to primary max report_date (later max)", { + ref_dates <- seq(as.Date("2023-01-01"), as.Date("2023-01-10"), by = "day") + primary <- make_daily_tri(ref_dates, 1:5) + aux <- make_daily_tri(ref_dates, 1:3) + # primary max report_date = 2023-01-10 + 5 = 2023-01-15 + # aux max report_date = 2023-01-10 + 3 = 2023-01-13 + + result <- data_preprocessing( + primary, + value_col = "value", refd_col = "reference_date", lag_col = "lag", + ref_lag = 6L, temporal_resol = "daily", + aux_triangles = list(beds = aux) + ) + + # rows with lag 4 and 5 (report_dates past aux max) should be present + expect_true(any(result$lag == 4L)) + expect_true(any(result$lag == 5L)) + expect_true("beds_value_raw" %in% colnames(result)) + expect_false(anyNA(result[["beds_value_raw"]])) +}) From 787882839448f04e4d2ebef4767c2ccd1f60db29 Mon Sep 17 00:00:00 2001 From: dsweber2 Date: Tue, 28 Jul 2026 11:34:11 -0500 Subject: [PATCH 5/5] guard against either being empty --- R/feature_engineering.R | 10 +++++++ tests/testthat/test-feature_engineering.R | 36 +++++++++++++++++++++++ 2 files changed, 46 insertions(+) diff --git a/R/feature_engineering.R b/R/feature_engineering.R index 03710e6..ed35057 100644 --- a/R/feature_engineering.R +++ b/R/feature_engineering.R @@ -298,6 +298,12 @@ process_aux_triangle <- function(df, name, lagged_term_list, temporal_resol, smo 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 { @@ -451,6 +457,10 @@ data_preprocessing <- function(df, value_col, refd_col, lag_col, ref_lag, 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") diff --git a/tests/testthat/test-feature_engineering.R b/tests/testthat/test-feature_engineering.R index 13b81c4..fc79e77 100644 --- a/tests/testthat/test-feature_engineering.R +++ b/tests/testthat/test-feature_engineering.R @@ -392,3 +392,39 @@ test_that("data_preprocessing forward-fills aux to primary max report_date (late expect_true("beds_value_raw" %in% colnames(result)) expect_false(anyNA(result[["beds_value_raw"]])) }) + +test_that("data_preprocessing returns 0 rows when aux has no data for the geo", { + ref_dates <- seq(as.Date("2023-01-01"), as.Date("2023-01-15"), by = "day") + primary <- make_daily_tri(ref_dates, 1:4) + empty_aux <- data.frame( + reference_date = as.Date(character()), + report_date = as.Date(character()), + lag = integer(), + value = numeric() + ) + + result <- data_preprocessing( + primary, + value_col = "value", refd_col = "reference_date", lag_col = "lag", + ref_lag = 5L, temporal_resol = "daily", + aux_triangles = list(beds = empty_aux) + ) + + expect_equal(nrow(result), 0L) +}) + +test_that("process_aux_triangle returns zero-row df with correct columns when aux is empty", { + empty_aux <- data.frame( + reference_date = as.Date(character()), + report_date = as.Date(character()), + lag = integer(), + value = numeric() + ) + lagged_term_list <- c(1L, 7L) + + result <- process_aux_triangle(empty_aux, "beds", lagged_term_list, "daily", TRUE) + + expect_equal(nrow(result), 0L) + expect_true(all(aux_feature_names("beds", lagged_term_list) %in% colnames(result))) + expect_true(all(c("reference_date", "report_date", "lag") %in% colnames(result))) +})