push scoring to the remote

Author: dsweber2

.gitignore

@@ -18,3 +18,7 @@ reports/report.md
 cache/
 data/
 .vscode/
+.envrc
+.nhsn_covid_cache.parquet
+.nhsn_flu_cache.parquet
+meta/

R/targets/score_targets.R

@@ -2,7 +2,7 @@ get_external_forecasts <- function(disease) {
   locations_crosswalk <- get_population_data() %>%
     select(state_id, state_code) %>%
     filter(state_id != "usa")
-  arrow::read_parquet(paste::paste("data/forecasts/{disease}_hosp_forecasts.parquet")) %>%
+  arrow::read_parquet(glue::glue("data/forecasts/{disease}_hosp_forecasts.parquet")) %>%
     filter(output_type == "quantile") %>%
     select(forecaster, geo_value = location, forecast_date, target_end_date, quantile = output_type_id, value) %>%
     inner_join(locations_crosswalk, by = c("geo_value" = "state_code")) %>%
@@ -49,7 +49,7 @@ score_forecasts <- function(nhsn_latest_data, joined_forecasts_and_ensembles) {
 }
 
 
-render_score_plot <- function(disease) {
+render_score_plot <- function(score_report_rmd, scores, forecast_dates, disease) {
   rmarkdown::render(
     score_report_rmd,
     params = list(

R/utils.R

@@ -364,6 +364,30 @@ update_site <- function(sync_to_s3 = TRUE) {
     prod_reports_index <- which(grepl("## Production Reports", report_md_content)) + 1
     report_md_content <- append(report_md_content, report_link, after = prod_reports_index)
   }
+  # add scoring notebooks if they exist
+  score_files <- dir_ls(reports_dir, regexp = ".*_backtesting_2024_2025_on_.*.html")
+  if (length(score_files) > 0) {
+    # a tibble of all score files, along with their generation date and disease
+    score_table <- tibble(
+      filename = score_files,
+      dates = str_match_all(filename, "[0-9]{4}-..-..")
+    ) %>%
+      unnest_wider(dates, names_sep = "_") %>%
+      rename(generation_date = dates_1) %>%
+      mutate(
+        generation_date = ymd(generation_date),
+        disease = str_match(filename, "flu|covid")
+      )
+    used_files <- score_table %>%
+      group_by(disease) %>%
+      slice_max(generation_date)
+    # iterating over the diseases
+    for (row_num in seq_along(used_files$filename)) {
+      scoring_index <- which(grepl("### Scoring this season", report_md_content)) + 1
+      score_link <- sprintf("- [%s Scoring, Rendered %s](%s)", str_to_title(used_files$disease[[row_num]]), used_files$generation_date[[row_num]], used_files$filename[[row_num]])
+      report_md_content <- append(report_md_content, score_link, after = scoring_index)
+    }
+  }
 
   # Write the updated content to report.md
   report_md_path <- path(reports_dir, "report.md")

reports/template.md

@@ -5,6 +5,9 @@
 ## Production Reports
 
 
+### Scoring this season
+
+
 ## Exploration Reports
 
 - [NHSN 2024-2025 Data Analysis](new_data.html)

scripts/covid_hosp_prod.R

@@ -14,6 +14,8 @@ forecast_generation_dates <- Sys.Date()
 # Usually, the forecast_date is the same as the generation date, but you can
 # override this. It should be a Wednesday.
 forecast_dates <- round_date(forecast_generation_dates, "weeks", week_start = 3)
+
+
 # forecast_generation_date needs to follow suit, but it's more complicated
 # because sometimes we forecast on Thursday.
 # forecast_generation_dates <- c(as.Date(c("2024-11-20", "2024-11-27", "2024-12-04", "2024-12-11", "2024-12-18", "2024-12-26", "2025-01-02")), seq.Date(as.Date("2025-01-08"), Sys.Date(), by = 7L))
@@ -412,15 +414,7 @@ if (backtest_mode) {
     tar_target(
       external_forecasts,
       command = {
-        locations_crosswalk <- get_population_data() %>%
-          select(state_id, state_code) %>%
-          filter(state_id != "usa")
-        arrow::read_parquet("data/forecasts/covid_hosp_forecasts.parquet") %>%
-          filter(output_type == "quantile") %>%
-          select(forecaster, geo_value = location, forecast_date, target_end_date, quantile = output_type_id, value) %>%
-          inner_join(locations_crosswalk, by = c("geo_value" = "state_code")) %>%
-          mutate(geo_value = state_id) %>%
-          select(forecaster, geo_value, forecast_date, target_end_date, quantile, value)
+        get_external_forecasts("covid")
       }
     ),
     tar_combine(
@@ -433,58 +427,13 @@ if (backtest_mode) {
     tar_target(
       name = scores,
       command = {
-        truth_data <-
-          nhsn_latest_data %>%
-          select(geo_value, target_end_date = time_value, oracle_value = value) %>%
-          left_join(
-            get_population_data() %>%
-              select(state_id, state_code),
-            by = c("geo_value" = "state_id")
-          ) %>%
-          drop_na() %>%
-          rename(location = state_code) %>%
-          select(-geo_value)
-        forecasts_formatted <-
-          joined_forecasts_and_ensembles %>%
-          format_scoring_utils(disease = "covid")
-        scores <- forecasts_formatted %>%
-          filter(location != "US") %>%
-          hubEvals::score_model_out(
-            truth_data,
-            metrics = c("wis", "ae_median", "interval_coverage_50", "interval_coverage_90"),
-            summarize = FALSE,
-            by = c("model_id")
-          )
-        scores %>%
-          left_join(
-            get_population_data() %>%
-              select(state_id, state_code),
-            by = c("location" = "state_code")
-          ) %>%
-          rename(
-            forecaster = model_id,
-            forecast_date = reference_date,
-            ahead = horizon,
-            geo_value = state_id
-          ) %>%
-          select(-location)
+        score_forecasts(nhsn_latest_data, joined_forecasts_and_ensembles)
       }
     ),
     tar_target(
       name = score_plot,
       command = {
-        rmarkdown::render(
-          score_report_rmd,
-          params = list(
-            scores = scores,
-            forecast_dates = forecast_dates,
-            disease = "covid"
-          ),
-          output_file = here::here(
-            "reports",
-            sprintf("covid_backtesting_2024_2025_on_%s.html", as.Date(Sys.Date()))
-          )
-        )
+        render_score_plot(score_report_rmd, scores, forecast_dates, "covid")
       },
       cue = tar_cue("always")
     )

scripts/flu_hosp_prod.R

@@ -467,15 +467,7 @@ if (backtest_mode) {
     tar_target(
       external_forecasts,
       command = {
-        locations_crosswalk <- get_population_data() %>%
-          select(state_id, state_code) %>%
-          filter(state_id != "usa")
-        arrow::read_parquet("data/forecasts/flu_hosp_forecasts.parquet") %>%
-          filter(output_type == "quantile") %>%
-          select(forecaster, geo_value = location, forecast_date, target_end_date, quantile = output_type_id, value) %>%
-          inner_join(locations_crosswalk, by = c("geo_value" = "state_code")) %>%
-          mutate(geo_value = state_id) %>%
-          select(forecaster, geo_value, forecast_date, target_end_date, quantile, value)
+        get_external_forecasts("flu")
       }
     ),
     tar_combine(
@@ -488,31 +480,18 @@ if (backtest_mode) {
     tar_target(
       name = scores,
       command = {
-        truth_data <- nhsn_latest_data %>%
-          select(geo_value, target_end_date = time_value, true_value = value) %>%
-          mutate(target_end_date = target_end_date + 3)
-        joined_forecasts_and_ensembles %>%
-          select(-source) %>%
-          rename("model" = "forecaster", "prediction" = "value") %>%
-          evaluate_predictions(forecasts = ., truth_data = truth_data) %>%
-          rename("forecaster" = "model")
+        nhsn_latest_end_of_week <-
+          nhsn_latest_data %>%
+          mutate(
+            time_value = ceiling_date(time_value, unit = "week", week_start = 6)
+          )
+        score_forecasts(nhsn_latest_end_of_week, joined_forecasts_and_ensembles)
       }
     ),
     tar_target(
       name = score_plot,
       command = {
-        rmarkdown::render(
-          score_report_rmd,
-          params = list(
-            scores = scores,
-            forecast_dates = forecast_dates,
-            disease = "flu"
-          ),
-          output_file = here::here(
-            "reports",
-            sprintf("flu_backtesting_2024_2025_on_%s.html", as.Date(Sys.Date()))
-          )
-        )
+        render_score_plot(score_report_rmd, scores, forecast_dates, "flu")
       }
     )
   )