dev/main comparison

.Rbuildignore

@@ -1,2 +1,4 @@
+^renv$
+^renv\.lock$
 ^.*\.Rproj$
 ^\.Rproj\.user$

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)
@@ -33,6 +34,7 @@ export(add_weights_related)
 export(assert)
 export(create_dir_not_exist)
 export(create_params_list)
+export(create_target_lookup)
 export(cv_revision_forecast)
 export(data_filteration)
 export(data_preprocessing)
@@ -68,7 +70,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/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) {
@@ -94,17 +96,11 @@ revision_forecast <- function(train_data, test_data, taus,
 
   test_data_list <- list()
 
-  if (train_models) {
-    sqrt_max_raw <- sqrt(max(train_data$value_7dav, na.rm=TRUE))
-    train_result <- add_sqrtscale(train_data, sqrt_max_raw)
-    train_data <- train_result$data
-    kept_bins <- train_result$kept_bins
-    #for (col in kept_bins) {
-    #  proportion <- sum(train_data[[col]]) / nrow(train_data)
-    #  cat(sprintf("Sum of %s: %.4f\n", col, proportion))
-    #}
-    train_data <- train_data[, c(basic_cols, params_list, extra_cols, kept_bins, response)] %>% drop_na()
-  }
+  sqrt_max_raw <- sqrt(max(train_data$value_7dav, na.rm=TRUE))
+  train_result <- add_sqrtscale(train_data, sqrt_max_raw)
+  train_data <- train_result$data
+  kept_bins <- train_result$kept_bins
+  train_data <- train_data[, c(basic_cols, params_list, extra_cols, kept_bins, response)] %>% drop_na()
 
   # pre-process the test data with max_raw
   if (make_predictions) {
@@ -317,6 +313,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 +324,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)
 
@@ -340,8 +338,8 @@ DelphiRF <- function(df, testing_start_date, taus=TAUS,
 
       # Detect weekly spacing
       if (length(lag_diffs) == 1 && lag_diffs == 7) {
+        if (temporal_resol != "weekly") message("Auto-detected weekly temporal resolution from lag spacing.")
         temporal_resol <- "weekly"
-        message("Auto-detected weekly temporal resolution from lag spacing.")
       }
     }
   }
@@ -388,13 +386,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)
@@ -210,11 +219,17 @@ exponentiate_preds <- function(test_data, taus) {
 get_model <- function(model_path, train_data, covariates, response, tau,
                       sqrt_max_raw, kept_bins,
                       lambda, gamma, lp_solver, train_models) {
-  if (train_models || !file.exists(model_path)) {
-    if (!train_models && !file.exists(model_path)) {
-      warning(str_interp("user requested use of cached model but file {model_path}"),
-        " does not exist; training new model")
+  if (!train_models && !file.exists(model_path)) {
+    candidates <- list.files(dirname(model_path), pattern = "\\.rds$", full.names = TRUE)
+    if (length(candidates) > 0L) {
+      alt_path <- sort(candidates, decreasing = TRUE)[1L]
+      message(sprintf("Exact model not found; loading most recent cached model: %s", basename(alt_path)))
+      return(readRDS(alt_path))
     }
+    warning(str_interp("user requested use of cached model but file ${model_path} does not exist and no alternatives found; training new model"))
+    train_models <- TRUE
+  }
+  if (train_models || !file.exists(model_path)) {
     # Quantile regression
     vec_7dav <- train_data[["value_7dav_diff"]]
     vec_slope <- train_data[["value_slope_diff"]]
@@ -332,27 +347,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
 }

R/preprocessing.R

@@ -36,6 +36,21 @@ fill_rows <- function(df, refd_col, lag_col, min_refd, max_refd, ref_lag) {
   return (df_new)
 }
 
+normalize_weekly_observations <- function(df, refd_col, lag_col) {
+  if (nrow(df) == 0) return(df)
+  original_reference_date <- as.Date(df[[refd_col]])
+  original_report_date <- original_reference_date + df[[lag_col]]
+  epiweek_end <- function(x) x + ((6L - as.POSIXlt(x)$wday) %% 7L)
+  df[[refd_col]] <- epiweek_end(original_reference_date)
+  df$report_date <- epiweek_end(original_report_date)
+  df[[lag_col]] <- as.numeric(df$report_date - df[[refd_col]])
+  ordering <- order(df[[refd_col]], df$report_date,
+                    original_report_date, original_reference_date)
+  df <- df[ordering, , drop = FALSE]
+  key <- paste(df[[refd_col]], df$report_date, sep = "\r")
+  df[!duplicated(key, fromLast = TRUE), , drop = FALSE]
+}
+
 #' Fill missing updates in a time series dataset with lagged values
 #' Get pivot table, filling NANs. If there is no update on issue date D but
 #' previous reports exist for issue date D_p < D, all the dates between
@@ -56,14 +71,19 @@ 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) {
     return(data.frame())
   }
   df$report_date <- df[[refd_col]] + df[[lag_col]]
 
+  if (temporal_resol == "weekly") {
+    df <- normalize_weekly_observations(df, refd_col, lag_col)
+  }
+
   # Generate a sequence of all possible dates
   if (temporal_resol == "daily"){
     all_reference_dates <- seq(min(df[[refd_col]]), max(df[[refd_col]]), by = "day")
@@ -73,21 +93,19 @@ fill_missing_updates <- function(df, value_col, refd_col, lag_col, temporal_reso
     all_reference_dates <- seq(min(df[[refd_col]]), max(df[[refd_col]]), by = "7 days")
     all_report_dates <- seq(min(df[["report_date"]]), max(df[["report_date"]]), by = "7 days")
 
-    # Check if all reference dates in df are within the generated sequence
-    if (!all(df[[refd_col]] %in% all_reference_dates)) {
-      stop("The reference dates do not regularly have a gap of 7 days. Some reference dates will be ignored. Please check your input data.")
-    }
-
-    # Check if all reference dates in df are within the generated sequence
-    if (!all(df$report_date %in% all_report_dates)) {
-      stop("The report dates do not regularly have a gap of 7 days. Some report dates will be ignored. Please check your input data.")
-    }
-
     gap <- 7
   } else {
     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,
@@ -118,10 +136,10 @@ fill_missing_updates <- function(df, value_col, refd_col, lag_col, temporal_reso
   backfill_df <- pivot_df %>%
     tidyr::pivot_longer(-report_date, values_to = "value_raw", names_to = refd_col) %>%
     mutate(
-      reference_date := as.Date(.data[[refd_col]]),
+      reference_date = as.Date(.data[[refd_col]]),
       lag = as.numeric(report_date - reference_date)
     ) %>%
-    filter(lag>=0)
+    filter(lag >= 0)
   return (as.data.frame(backfill_df))
 }
 
@@ -177,6 +195,3 @@ get_weekofmonth <- function(date) {
 
   return(week_number)
 }
-
-
-

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