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

@@ -32,6 +32,37 @@ add_dayofweek <- function(df, time_col, suffix, wd = WEEKDAYS_ABBR) {
 }
 
 
+#' Add grouped day-of-week one-hot columns
+#'
+#' Creates one binary column per group in `onehot_weekdays`. Each group is a
+#' character vector of day abbreviations (from `WEEKDAYS_ABBR`); a row is 1 if
+#' the date falls on any day in the group.  Column names are derived from the
+#' list names when present, otherwise by concatenating the day abbreviations.
+#'
+#' @param df A data frame containing the date column.
+#' @param time_col Name of the date column.
+#' @param suffix Column name suffix (e.g. `"_ref"` or `"_issue"`).
+#' @param onehot_weekdays Named or unnamed list of character vectors, each
+#'   specifying a group of days (e.g. `list(Mon=c("Mon"), Weekends=c("Sat","Sun"))`).
+#'
+#' @return `df` with one additional integer column per group.
+#' @export
+add_grouped_dayofweek <- function(df, time_col, suffix, onehot_weekdays) {
+  df <- df %>% mutate({{ time_col }} := as.Date(.data[[time_col]]))
+  dayofweek <- as.numeric(format(df[[time_col]], format = "%u"))
+  group_names <- if (!is.null(names(onehot_weekdays))) {
+    names(onehot_weekdays)
+  } else {
+    vapply(onehot_weekdays, function(grp) paste0(grp, collapse = ""), character(1))
+  }
+  for (ii in seq_along(onehot_weekdays)) {
+    day_indices <- match(onehot_weekdays[[ii]], WEEKDAYS_ABBR)
+    df[[paste0(group_names[ii], suffix)]] <- as.integer(dayofweek %in% day_indices)
+  }
+  df
+}
+
+
 #' Add one-hot encoding for week of the month based on issue date
 #'
 #' This function calculates the week of the month for each date in the specified
@@ -188,8 +219,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(
@@ -243,31 +279,98 @@ add_log_transformed <- function(df, lagged_term_list) {
 #' @param lag_col Column name representing the lag between the reference and issue date.
 #' @param temporal_resol A string indicating the temporal resolution ("daily" or "weekly").
 #'                       Defaults to "daily".
+#' @param onehot_weekdays Named or unnamed list of character vectors defining day groups
+#'   for one-hot encoding.  Each group becomes one binary column per date column.
+#'   Names, when present, are used as column prefixes; otherwise day abbreviations are
+#'   concatenated.  Defaults to `list(Mon=c("Mon"), Weekends=c("Sat","Sun"))`.
 #'
 #' @details
-#' - If `temporal_resol` is "daily", one-hot encoded day-of-week columns are added
+#' - If `temporal_resol` is "daily", one-hot encoded day-of-week group columns are added
 #'   for both `refd_col` (reference date) and `"report_date"`.
 #' - One-hot encoded week-of-month columns are added for `"report_date"` in all cases.
 #'
 #' @return A modified data frame with additional date-related feature columns.
 #'
 #' @export
-add_params_for_dates <- function(df, refd_col, lag_col, temporal_resol="daily") {
+add_params_for_dates <- function(df, refd_col, lag_col, temporal_resol = "daily",
+                                  onehot_weekdays = list(Mon = c("Mon"), Weekends = c("Sat", "Sun"))) {
   df$report_date <- df[[refd_col]] + df[[lag_col]]
-  if (temporal_resol=="daily"){
-    # Add columns for day-of-week effect
-    df <- add_dayofweek(df, refd_col, "_ref", WEEKDAYS_ABBR)
-    df <- add_dayofweek(df, "report_date", "_issue", WEEKDAYS_ABBR)
-    # Add columns for weekends
-    df$Weekends_issue <- as.integer(df$Sat_issue == 1 | df$Sun_issue == 1)
-    df$Weekends_ref <- as.integer(df$Sat_ref == 1 | df$Sun_ref == 1)
+  if (temporal_resol == "daily") {
+    df <- add_grouped_dayofweek(df, refd_col, "_ref", onehot_weekdays)
+    df <- add_grouped_dayofweek(df, "report_date", "_issue", onehot_weekdays)
   }
-  # Add columns for week-of-month effect
   df <- add_weekofmonth(df, "report_date", WEEK_ISSUES)
+  return(as.data.frame(df))
+}
 
-  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 +390,23 @@ 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,
+                               onehot_weekdays = list(Mon = c("Mon"), Weekends = c("Sat", "Sun"))) {
   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.")
@@ -367,7 +479,26 @@ 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)
+  merged_df <- add_params_for_dates(merged_df, "reference_date", "lag", temporal_resol, onehot_weekdays)
+
+  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="",
@@ -49,9 +50,10 @@ revision_forecast <- function(train_data, test_data, taus,
                               indicator="testdata", signal="",
                               geo_level="state", signal_suffix="",
                               training_end_date="",
-                              training_days =365,
+                              training_days=365,
                               train_models = TRUE,
-                              make_predictions=TRUE) {
+                              make_predictions=TRUE,
+                              onehot_weekdays = list(Mon = c("Mon"), Weekends = c("Sat", "Sun"))) {
 
 
 
@@ -77,7 +79,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, onehot_weekdays, extra_params)
   }
 
   if (smoothed_target) {
@@ -317,6 +319,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,
@@ -327,7 +330,8 @@ DelphiRF <- function(df, testing_start_date, taus=TAUS,
                      training_end_date="",
                      training_days=365,
                      train_models = TRUE,
-                     make_predictions = TRUE) {
+                     make_predictions = TRUE,
+                     onehot_weekdays = list(Mon = c("Mon"), Weekends = c("Sat", "Sun"))) {
 
   testing_start_date <- as.Date(testing_start_date)
 
@@ -388,13 +392,13 @@ 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,
                                  signal_suffix, as.character(testing_start_date),
                                  training_days, train_models,
-                                 make_predictions)
+                                 make_predictions, onehot_weekdays)
 
     test_data_list <- append(test_data_list, list(results))
   }

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,27 +341,26 @@ 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,
+                               onehot_weekdays = list(Mon = c("Mon"), Weekends = c("Sat", "Sun")),
+                               extra_params = NULL) {
   params_list <- c(
     WEEK_ISSUES[1],
     Y7DAV,
     paste0("log_value_7dav_lag", lagged_term_list),
     paste0("log_delta_value_7dav_lag", lagged_term_list)
   )
-  # Include log lag adjustments if multiple lags exist
-  if (length(unique(train_data$lag)) > 1){
+  if (length(unique(train_data$lag)) > 1) {
     params_list <- c(params_list, LOG_LAG)
   }
 
-  dayofweek <- c("Mon", "Weekends")
-  extra_params_for_daily <- c(
-    paste0(dayofweek, "_ref"),
-    paste0(dayofweek, "_issue")
-  )
-
-  if (temporal_resol == "daily"){
-    return (c(params_list, extra_params_for_daily))
+  group_names <- if (!is.null(names(onehot_weekdays))) {
+    names(onehot_weekdays)
   } else {
-    return(params_list)
+    vapply(onehot_weekdays, function(grp) paste0(grp, collapse = ""), character(1))
   }
+  extra_params_for_daily <- c(paste0(group_names, "_ref"), paste0(group_names, "_issue"))
+
+  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
 }