Issue 188: Benchmark infection feedback performance (#189)

_targets_eval_postprocessing.R

@@ -1464,6 +1464,22 @@ benchmarks <- list(
       wwinference_version = eval_config$wwinference_version,
       overwrite_benchmark = eval_config$overwrite_benchmark
     )
+  ),
+  tar_target(
+    name = plot_benchmark_by_loc,
+    command = plot_benchmarks(
+      grouping_var = "location",
+      benchmark_scope = "all_forecasts",
+      benchmark_dir = benchmark_config$benchmark_dir
+    )
+  ),
+  tar_target(
+    name = plot_benchmark_by_forecast_date,
+    command = plot_benchmarks(
+      grouping_var = "forecast_date",
+      benchmark_scope = "all_forecasts",
+      benchmark_dir = benchmark_config$benchmark_dir
+    )
   )
 )
 

_targets_subset_benchmarking.R

@@ -82,7 +82,7 @@ combined_targets <- list(
 # Benchmarking----------------------------------------------------------
 benchmarks <- list(
   tar_target(
-    name = write_benchmark_table_subset_run,
+    name = benchmark_table_subset_run,
     command = benchmark_performance(
       ww_scores = ww_scores,
       hosp_scores = hosp_scores,
@@ -91,6 +91,24 @@ benchmarks <- list(
       wwinference_version = benchmark_config$wwinference_version,
       overwrite_benchmark = benchmark_config$overwrite_benchmark
     )
+  ),
+  tar_target(
+    name = plot_benchmark_by_loc,
+    command = plot_benchmarks(
+      grouping_var = "location",
+      benchmark_scope = "subset_forecasts",
+      benchmark_dir = benchmark_config$benchmark_dir,
+      scores_list = benchmark_table_subset_run
+    )
+  ),
+  tar_target(
+    name = plot_benchmark_by_forecast_date,
+    command = plot_benchmarks(
+      grouping_var = "forecast_date",
+      benchmark_scope = "subset_forecasts",
+      benchmark_dir = benchmark_config$benchmark_dir,
+      scores_list = benchmark_table_subset_run
+    )
   )
 )
 

command_line_eval.R

@@ -5,7 +5,6 @@
 # in parallel on Azure batch
 
 library(argparser, quietly = TRUE)
-library(cfaforecastrenewalww)
 library(wweval)
 library(ggplot2)
 # some functions from plots.R complain about aes() function not existing if we don't load ggplot2

input/params.toml

@@ -41,9 +41,10 @@ offset_ref_initial_exp_growth_rate_prior_sd = 0.025
 
 autoreg_p_hosp_a = 1 # shape1 parameter of autoreg term on IHR(t) trend
 autoreg_p_hosp_b = 100 # shape2 parameter of autoreg term on IHR(t) trend
-eta_sd_sd = 0.01
-infection_feedback_prior_logmean = 6.37408 # log(mode) + q^2 mode = 500, q = 0.4
-infection_feedback_prior_logsd = 0.4
+eta_sd_sd = 0.0097 # this is the mean, sd is 0.0097
+eta_sd_mean = 0.0278
+infection_feedback_prior_logmean = 4.498 # mode=100 log(mode) + q^2 mode = 100, q = 1
+infection_feedback_prior_logsd = 0.636 # roughly 16 to 600 in 95% ci
 
 [hospital_admission_observation_process]
 # Hospitalization parameters (informative priors)

output/benchmarking/latest_subset_forecasts_by_forecast_date.tsv

@@ -1,7 +1,7 @@
 forecast_date	crps_hosp	crps_ww	bias_hosp	bias_ww	ae_hosp	ae_ww	wweval_commit_hash	wwinference_version	time_stamp
-2023-10-16	5.145227727175137	4.405852136963368	0.48245657894736843	0.45490394736842105	6.950701593150161	6.041583425908238	d23b28a	v0.1.0	2024-10-24T17:48:21Z
-2023-11-13	4.66700940996615	4.019961671884518	-0.3123618421052632	-0.22545263157894738	6.145553133172695	5.324888683286438	d23b28a	v0.1.0	2024-10-24T17:48:21Z
-2023-12-11	8.859086840566519	8.837148573389927	-0.5840236842105263	-0.511771052631579	11.50992521124024	11.344388151022924	d23b28a	v0.1.0	2024-10-24T17:48:21Z
-2024-01-08	17.39089396390878	15.609039964313409	0.25387105263157894	0.3958921052631579	22.182123575526884	20.51478614874109	d23b28a	v0.1.0	2024-10-24T17:48:21Z
-2024-02-05	3.324198357881111	5.079010312322833	0.26249342105263157	0.3663105263157895	4.572098587350526	7.344139936498404	d23b28a	v0.1.0	2024-10-24T17:48:21Z
-2024-03-04	3.279328388379696	4.706364986307469	0.49960394736842106	0.6323973684210527	4.718323317619102	6.606204242479576	d23b28a	v0.1.0	2024-10-24T17:48:21Z
+2023-10-16	5.874998928099601	4.842755313016674	0.4281315789473684	0.40362236842105265	8.202210281292677	6.916940992684773	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+2023-11-13	5.514058896168103	4.2595719435319985	-0.4940644736842105	-0.3491868421052632	7.648424162053213	5.968291847862045	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+2023-12-11	5.367768349490639	5.48631763373991	-0.34426973684210527	-0.20431052631578947	8.157250251202527	7.768904286971368	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+2024-01-08	35.16983309770466	35.53006450254746	0.5046828947368421	0.5802315789473684	44.28701944409863	45.39614861289537	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+2024-02-05	2.773028914102106	3.711888477140153	-0.07066578947368421	0.15726973684210527	3.12398987017981	5.094208909464286	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+2024-03-04	2.751103566063376	4.452629873235709	0.3015644736842105	0.4553171052631579	3.7821969793270878	6.373191104611754	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z

output/benchmarking/latest_subset_forecasts_by_location.tsv

@@ -1,7 +1,7 @@
 crps_hosp	crps_ww	bias_hosp	bias_ww	ae_hosp	ae_ww	location	wweval_commit_hash	wwinference_version	time_stamp
-1.0465392954054191	1.0437418123363529	-0.052574561403508774	-0.05236184210526316	1.430541441794062	1.4249102564297214	AK	d23b28a	v0.1.0	2024-10-24T17:48:21Z
-11.624616908337584	11.001390757513024	0.23765679824561403	0.30562280701754385	15.2169931884645	14.921735155883097	MA	d23b28a	v0.1.0	2024-10-24T17:48:21Z
-1.9743211076933298	1.7647838214371094	0.10642543859649123	0.11702741228070176	2.6690277913463327	2.360288987441267	NH	d23b28a	v0.1.0	2024-10-24T17:48:21Z
-17.371192144075163	17.18502224109214	0.16711842105263158	0.24489583333333334	22.505788770781994	22.54254318981423	NJ	d23b28a	v0.1.0	2024-10-24T17:48:21Z
-3.5381177843863316	4.552876071939314	0.0430734649122807	0.31171600877192984	4.9099199893294525	6.397181233712242	WA	d23b28a	v0.1.0	2024-10-24T17:48:21Z
-7.110957447979565	7.109562940863587	0.10033991228070176	0.18538004385964912	9.346454236343268	9.529331764656112	all	d23b28a	v0.1.0	2024-10-24T17:48:21Z
+0.9048110671853405	0.8654847567756644	-0.07057236842105263	-0.07581030701754386	1.2465896546056678	1.2173595112401607	AK	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+14.762354885967529	12.681244008639368	0.18578947368421053	0.30892434210526315	19.250396069656812	16.835413776804067	MA	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+2.100055219513157	1.6857065406611689	-0.036726973684210525	0.038880482456140356	2.984162257895195	2.31003750944109	NH	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+26.360153847097493	25.59449289371047	0.1790548245614035	0.2609703947368421	33.681565511153316	33.56127504975967	NJ	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+3.7482847732602225	7.7424282528899155	0.013604166666666669	0.3361546052631579	5.504862330150637	10.673985614829668	WA	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+9.575131958604748	9.713871290535318	0.05422982456140351	0.17382390350877194	12.533515164692325	12.919614292414932	all	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z

output/benchmarking/subset_forecasts_by_forecast_date.tsv

@@ -1,6 +1,18 @@
 forecast_date	crps_hosp	crps_ww	bias_hosp	bias_ww	ae_hosp	ae_ww	wweval_commit_hash	wwinference_version	time_stamp
-2023-10-16	5.145227727175137	4.405852136963368	0.48245657894736843	0.45490394736842105	6.950701593150161	6.041583425908238	d23b28a	v0.1.0	2024-10-24T17:48:21Z
-2023-11-13	4.66700940996615	4.019961671884518	-0.3123618421052632	-0.22545263157894738	6.145553133172695	5.324888683286438	d23b28a	v0.1.0	2024-10-24T17:48:21Z
+2023-10-16	5.874998928099601	4.842755313016674	0.4281315789473684	0.40362236842105265	8.202210281292677	6.916940992684773	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+2023-11-13	5.514058896168103	4.2595719435319985	-0.4940644736842105	-0.3491868421052632	7.648424162053213	5.968291847862045	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+2023-12-11	5.367768349490639	5.48631763373991	-0.34426973684210527	-0.20431052631578947	8.157250251202527	7.768904286971368	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+2024-01-08	35.16983309770466	35.53006450254746	0.5046828947368421	0.5802315789473684	44.28701944409863	45.39614861289537	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+2024-02-05	2.773028914102106	3.711888477140153	-0.07066578947368421	0.15726973684210527	3.12398987017981	5.094208909464286	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+2024-03-04	2.751103566063376	4.452629873235709	0.3015644736842105	0.4553171052631579	3.7821969793270878	6.373191104611754	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+2023-10-16	9.243950678640694	6.806804188258501	0.5824578947368421	0.4910565789473684	12.075811587773888	9.080835846088256	51ddee5	227-inf-feedback	2024-10-25T17:16:51Z
+2023-11-13	6.798783532681564	4.377860426544979	-0.572275	-0.32701842105263157	8.806163742433748	5.917972294575069	51ddee5	227-inf-feedback	2024-10-25T17:16:51Z
+2023-12-11	5.578345847009017	5.404499318529154	-0.22773947368421052	-0.12964736842105262	7.0379402652994925	6.487141426491198	51ddee5	227-inf-feedback	2024-10-25T17:16:51Z
+2024-01-08	73.02250974305282	64.10560452659152	0.5030894736842105	0.5657434210526315	90.10392730630215	79.63998299207216	51ddee5	227-inf-feedback	2024-10-25T17:16:51Z
+2024-02-05	3.2543113736957956	4.331650777414377	-0.16031315789473685	0.14943684210526317	4.889963029507771	5.687676001142328	51ddee5	227-inf-feedback	2024-10-25T17:16:51Z
+2024-03-04	2.4436722259876116	4.390295509081186	0.2651552631578947	0.5207526315789474	3.274852617671909	6.112645792027634	51ddee5	227-inf-feedback	2024-10-25T17:16:51Z
+2023-10-16	5.145227727175137	4.405852136963368	0.4824565789473684	0.45490394736842105	6.950701593150161	6.041583425908238	d23b28a	v0.1.0	2024-10-24T17:48:21Z
+2023-11-13	4.66700940996615	4.019961671884518	-0.3123618421052632	-0.22545263157894735	6.145553133172695	5.324888683286438	d23b28a	v0.1.0	2024-10-24T17:48:21Z
 2023-12-11	8.859086840566519	8.837148573389927	-0.5840236842105263	-0.511771052631579	11.50992521124024	11.344388151022924	d23b28a	v0.1.0	2024-10-24T17:48:21Z
 2024-01-08	17.39089396390878	15.609039964313409	0.25387105263157894	0.3958921052631579	22.182123575526884	20.51478614874109	d23b28a	v0.1.0	2024-10-24T17:48:21Z
 2024-02-05	3.324198357881111	5.079010312322833	0.26249342105263157	0.3663105263157895	4.572098587350526	7.344139936498404	d23b28a	v0.1.0	2024-10-24T17:48:21Z

output/benchmarking/subset_forecasts_by_location.tsv

@@ -1,7 +1,19 @@
 crps_hosp	crps_ww	bias_hosp	bias_ww	ae_hosp	ae_ww	location	wweval_commit_hash	wwinference_version	time_stamp
-1.0465392954054191	1.0437418123363529	-0.052574561403508774	-0.05236184210526316	1.430541441794062	1.4249102564297214	AK	d23b28a	v0.1.0	2024-10-24T17:48:21Z
+0.9048110671853405	0.8654847567756644	-0.07057236842105263	-0.07581030701754386	1.2465896546056678	1.2173595112401607	AK	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+14.762354885967529	12.681244008639368	0.18578947368421053	0.30892434210526315	19.250396069656812	16.835413776804067	MA	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+2.100055219513157	1.6857065406611689	-0.036726973684210525	0.038880482456140356	2.984162257895195	2.31003750944109	NH	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+26.360153847097493	25.59449289371047	0.1790548245614035	0.2609703947368421	33.681565511153316	33.56127504975967	NJ	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+3.7482847732602225	7.7424282528899155	0.013604166666666669	0.3361546052631579	5.504862330150637	10.673985614829668	WA	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+9.575131958604748	9.713871290535318	0.05422982456140351	0.17382390350877194	12.533515164692325	12.919614292414932	all	11a9898	227-inf-feedback-posterior-mod-eta	2024-10-30T19:46:57Z
+1.0376521512900907	1.0118669366748054	-0.03699122807017543	-0.00765679824561403	1.4245526802625763	1.4264265157766278	AK	51ddee5	227-inf-feedback	2024-10-25T17:16:51Z
+24.626800317489625	19.49484775911998	0.2060657894736842	0.4117872807017544	30.689131156635455	24.78340014717931	MA	51ddee5	227-inf-feedback	2024-10-25T17:16:51Z
+2.9839338399865656	1.8872799196790224	-0.01010087719298246	0.04486293859649123	4.248055546108953	2.2765310519617756	NH	51ddee5	227-inf-feedback	2024-10-25T17:16:51Z
+50.50136383077716	41.39684785751621	0.1817642543859649	0.2782510964912281	62.26077007486285	51.60308220739496	NJ	51ddee5	227-inf-feedback	2024-10-25T17:16:51Z
+4.468227694679483	10.723086482359754	-0.01542543859649122	0.33135855263157893	6.534705999620968	14.01577203801786	WA	51ddee5	227-inf-feedback	2024-10-25T17:16:51Z
+16.723595566844583	14.902785791069954	0.0650625	0.2117206140350877	21.03144309149816	18.821042392066108	all	51ddee5	227-inf-feedback	2024-10-25T17:16:51Z
+1.0465392954054191	1.0437418123363529	-0.05257456140350877	-0.05236184210526316	1.430541441794062	1.4249102564297214	AK	d23b28a	v0.1.0	2024-10-24T17:48:21Z
 11.624616908337584	11.001390757513024	0.23765679824561403	0.30562280701754385	15.2169931884645	14.921735155883097	MA	d23b28a	v0.1.0	2024-10-24T17:48:21Z
-1.9743211076933298	1.7647838214371094	0.10642543859649123	0.11702741228070176	2.6690277913463327	2.360288987441267	NH	d23b28a	v0.1.0	2024-10-24T17:48:21Z
-17.371192144075163	17.18502224109214	0.16711842105263158	0.24489583333333334	22.505788770781994	22.54254318981423	NJ	d23b28a	v0.1.0	2024-10-24T17:48:21Z
+1.9743211076933296	1.7647838214371094	0.10642543859649124	0.11702741228070176	2.6690277913463327	2.360288987441267	NH	d23b28a	v0.1.0	2024-10-24T17:48:21Z
+17.371192144075163	17.18502224109214	0.16711842105263158	0.24489583333333337	22.505788770781994	22.54254318981423	NJ	d23b28a	v0.1.0	2024-10-24T17:48:21Z
 3.5381177843863316	4.552876071939314	0.0430734649122807	0.31171600877192984	4.9099199893294525	6.397181233712242	WA	d23b28a	v0.1.0	2024-10-24T17:48:21Z
 7.110957447979565	7.109562940863587	0.10033991228070176	0.18538004385964912	9.346454236343268	9.529331764656112	all	d23b28a	v0.1.0	2024-10-24T17:48:21Z

scratch/get_empirical_posterior_vals.R

@@ -0,0 +1,72 @@
+# Quick estimate of posterior parameters
+
+benchmark_config <- yaml::read_yaml(file.path(
+  "input", "config",
+  "eval", "benchmark_config.yaml"
+))
+
+vars <- c("eta_sd", "inf_feedback")
+eta_sd_draws <- tibble::tibble()
+inf_feedback_draws <- tibble::tibble()
+
+for (i in seq_along(benchmark_config$forecast_date_hosp)) {
+  this_location <- benchmark_config$location_hosp[i]
+  this_forecast_date <- benchmark_config$forecast_date_hosp[i]
+  this_scenario <- "no_wastewater"
+  for (j in seq_along(vars)) {
+    fp_var <- wweval::get_filepath(benchmark_config$output_dir,
+      scenario = this_scenario,
+      forecast_date = this_forecast_date,
+      model_type = "hosp",
+      location = this_location,
+      output_type = vars[j],
+      file_extension = "tsv"
+    )
+
+    these_var_draws <- readr::read_tsv(fp_var)
+    var_draws <- these_var_draws |>
+      dplyr::mutate(
+        location = this_location,
+        forecast_date = this_forecast_date
+      )
+    if (vars[j] == "eta_sd") {
+      eta_sd_draws <- dplyr::bind_rows(eta_sd_draws, var_draws)
+    }
+    if (vars[j] == "inf_feedback") {
+      inf_feedback_draws <- dplyr::bind_rows(
+        inf_feedback_draws,
+        var_draws
+      )
+    }
+  } # end loop over vars
+} # end loop over forecast date-locations
+
+# Get empirical mean, sd, logmean, and logsd------------------------------
+## eta_sd---------------------------------------------------------------
+mean_eta_sd <- mean(eta_sd_draws$eta_sd)
+sd_eta_sd <- sd(eta_sd_draws$eta_sd)
+
+message("Empirical mean of RW step size across 5 locations: ", mean_eta_sd)
+message("Empirical sd of RW step size across 5 locations: ", sd_eta_sd)
+
+## inf_feedback----------------------------------------------------------
+logmean_inf_feedback <- mean(log(inf_feedback_draws$infection_feedback))
+logsd_inf_feedback <- sd(log(inf_feedback_draws$infection_feedback))
+
+
+message(
+  "Empirical logmean of infection feedback across 5 locations: ",
+  logmean_inf_feedback
+)
+message(
+  "Empirical logsd of infection feedback across 5 locations: ",
+  logsd_inf_feedback
+)
+
+posterior_params <- list(
+  mean_eta_sd = mean_eta_sd,
+  sd_eta_sd = sd_eta_sd,
+  logmean_inf_feedback = logmean_inf_feedback,
+  logsd_inf_feedback = logsd_inf_feedback
+)
+yaml::write_yaml(posterior_params, "output/posterior_params.yaml")

src/setup_eval.R

@@ -32,5 +32,7 @@ write_eval_config(
   overwrite_summary_table = FALSE, # Set as TRUE if trying to get a baseline
   # score for all locations one forecast date
   overwrite_benchmark = TRUE, # Set as TRUE if want to save outputs of
-  wwinference_version = "v0.1.0"
+
+  # benchmarking in directory
+  wwinference_version = "227-inf-feedback-posterior-mod-eta"
 )

src/setup_subset_benchmarking.R

@@ -29,5 +29,5 @@ write_eval_config(
   overwrite_benchmark = TRUE, # Set as TRUE if want to save outputs of
   # benchmarking in directory,
   name_of_config = "benchmark_config",
-  wwinference_version = "v0.1.0"
+  wwinference_version = "227-inf-feedback-posterior-mod-eta"
 )

src/write_eval_config.R

@@ -6,9 +6,10 @@
 #' @param scenatios the scenarios (which will pertain to site ids) to
 #' run the model on
 #' @param config_dir the directory where we want to save the config file
+#' @param scenario_dir the directory where the files defining scenarios
+#' (default `.tsv` format) are located
 #' @param benchmark_dir the directory where to save the benchmarked performance
 #' for this run
-#' @param ms_fig_dir the directory to save the manuscript figures in
 #' @param eval_date the data of the evaluation dataset, in ISO YYYY-MM-DD format
 #' @param overwrite_summary_table Boolean indicating whether or not to overwrite
 #' internal summary table
@@ -191,6 +192,7 @@ write_eval_config <- function(locations, forecast_dates,
     adapt_delta = adapt_delta,
     max_treedepth = max_treedepth,
     seed = seed,
+    name_of_config = name_of_config,
     # Input delay distributions
     generation_interval = generation_interval,
     infection_feedback_pmf = generation_interval,

wweval/DESCRIPTION

@@ -90,5 +90,5 @@ Config/Needs/check: rcmdcheck, testthat
 RoxygenNote: 7.3.2
 Remotes:
     stan-dev/cmdstanr,
-    wwinference=CDCgov/ww-inference-model@v0.1.0
+    wwinference=CDCgov/ww-inference-model@227-inf-feedback-mod
 LazyData: true

wweval/NAMESPACE

@@ -87,6 +87,7 @@ export(nhsn_soda_query)
 export(order_horizons)
 export(order_periods)
 export(order_phases)
+export(plot_benchmarks)
 export(plot_components)
 export(plot_quantiles)
 export(pull_nhsn)