Updates to running, and generating figures

.gitignore

@@ -11,4 +11,5 @@ data/simulated/simulated_catchment_populations.csv
 data/simulated/simulated_facility_distances.csv
 data/simulated/simulated_hospital_counts.csv
 data/simulated/simulated_ww_data.csv
-R/bug_testing
+R/bug_testing
+figs/spatial_paper/

R/WA_realdata/1_nonspatial_realdata.R

@@ -1,12 +1,10 @@
 #Install pacman if missing
 if(!require(pacman)) install.packages("pacman")
 
-#If you have installed a NON-MAIN branch version of the package, then re-install the main branch
-prior <- packageDescription("wwinference")$GithubRef
-if(packageDescription("wwinference")$GithubRef != "HEAD"){
-  devtools::install_github("CDCgov/ww-inference-model")
+#If you have installed the NON-SPATIAL version of the package, then re-install the spatial branch
+if(packageDescription("wwinference")$GithubRef != "spatial-main"){
+  devtools::install_github("CDCgov/ww-inference-model@spatial-main")
 }
-current <- packageDescription("wwinference")$GithubRef
 
 #Load packages
 pacman::p_load(

R/WA_realdata/2_spatial_realdata.R

@@ -2,11 +2,9 @@
 if(!require(pacman)) install.packages("pacman")
 
 #If you have installed the NON-SPATIAL version of the package, then re-install the spatial branch
-prior <- packageDescription("wwinference")$GithubRef
 if(packageDescription("wwinference")$GithubRef != "spatial-main"){
   devtools::install_github("CDCgov/ww-inference-model@spatial-main")
 }
-current <- packageDescription("wwinference")$GithubRef
 
 #Load packages
 pacman::p_load(

R/functions/WA_spatial_run.R

@@ -63,7 +63,7 @@ WA_spatial_run <- function(
   these_runs <- 1:length(random_dates)
 
   if(!is.na(repeat_subset)){
-    these_runs <- as.numeric(unlist(strsplit(repeat_subset, ";")))
+    these_runs <- as.numeric(unlist(strsplit(as.character(repeat_subset), ";")))
   }
 
   # Loop through runs ----------------------------------------------------
@@ -617,6 +617,13 @@ WA_spatial_run <- function(
     group_by(site, forecast_or_fit, metric) %>%
     mutate(no_correlation_value = value[model == "No correlation"])
 
+  hosp_score_gathered <- model_score_raw_all %>%
+    gather(bias:se_mean,
+           key = "metric",
+           value = "value") %>%
+    group_by(forecast_or_fit, metric) %>%
+    mutate(no_correlation_value = value[model == "No correlation"])
+
   ww_score_gathered_format <- ww_score_gathered %>%
     group_by(site, metric, forecast_or_fit, model) %>%
     mutate(difference = 100 * (abs(no_correlation_value) - abs(value))/abs(value)) %>%
@@ -638,6 +645,23 @@ WA_spatial_run <- function(
       metric_model = paste(metric, model, sep = " - ")
     )
 
+  hosp_score_gathered_format <- hosp_score_gathered %>%
+    group_by(site, metric, forecast_or_fit, model) %>%
+    mutate(difference = 100 * (abs(no_correlation_value) - abs(value))/abs(value)) %>%
+    filter(model != "No correlation") %>%
+    mutate(metric_full = case_when(
+      metric == "mad" ~ "Mean absolute deviation",
+      metric == "bias" ~ "Bias",
+      metric == "dss" ~ "Dawid-Sebastiani score",
+      metric == "crps" ~ "Ranked Probability score",
+      metric == "ae_median" ~ "Absolute error (median)",
+      metric == "se_mean" ~ "Standard error (mean)",
+      metric == "log_score" ~ "Log score",
+      metric == "overprediction" ~ "Overprediction",
+      metric == "underprediction" ~ "Underprediction",
+      metric == "dispersion" ~ "Dispersion"
+    ))
+
   #Generate overall plot
   metric_better_plot <- ggplot(
     data = ww_score_gathered_format %>%
@@ -667,6 +691,54 @@ WA_spatial_run <- function(
     theme_bw() +
     facet_wrap(~model, ncol = 1)
 
+  metric_better_plot_hosp <- ggplot(
+    data = hosp_score_gathered_format %>%
+      filter(forecast_or_fit == "Forecast"),
+    mapping = aes(
+      x = model,
+      y = metric_full,
+      fill = difference
+    )
+  ) +
+    geom_tile() +
+    labs("% spatial improvement") +
+    colorspace::scale_fill_continuous_divergingx(
+      palette = 'RdBu', 
+      mid = 0.0, 
+      p1 = 0.5,
+      p2 = 0.75,
+      # l3 = 0,
+      # p3 = .5,
+      p4 = 1,
+      rev = T
+    ) +
+    labs(x = "Site",
+         y = "",
+         fill = "% metric improvement\nfrom wastewater data") +
+    theme_bw() +
+    facet_wrap(~model, ncol = 1)
+
+
+
+  ggplot(
+    data = hosp_score_gathered_format %>%
+      filter(forecast_or_fit == "Forecast" &
+               metric_full != "Bias"),
+    mapping = aes(
+      y = metric_full,
+      x = difference,
+      color = model
+    )
+  ) +
+    geom_point() +
+    labs("% spatial improvement") +
+    labs(y = "",
+         color = "",
+         x = "% metric improvement\nfrom wastewater data") +
+    theme_bw() +
+    theme(axis.text = element_text(angle = 0)) +
+    geom_vline(xintercept = 0, linetype = "dashed")
+
   # Output summary plots and statistics -------------------------------------
 
   #Output

R/functions/process_WA_data.R

@@ -82,7 +82,7 @@ process_WA_data <- function(
     ) %>%
     mutate(
       log_genome_copies_per_ml = log(simulated_raw_concentration),
-      log_lod = 5,
+      log_lod = 5,      #This is currently arbitrary - should be updated with correct value
       site = lab
     ) %>%
     left_join(

R/scripts/1_assess_nonspatial_state.R

@@ -1,12 +1,10 @@
 #Install pacman if missing
 if(!require(pacman)) install.packages("pacman")
 
-#If you have installed a NON-MAIN branch version of the package, then re-install the main branch
-prior <- packageDescription("wwinference")$GithubRef
-if(packageDescription("wwinference")$GithubRef != "HEAD"){
-  devtools::install_github("CDCgov/ww-inference-model")
+#If you have installed the NON-SPATIAL version of the package, then re-install the spatial branch
+if(packageDescription("wwinference")$GithubRef != "spatial-main"){
+  devtools::install_github("CDCgov/ww-inference-model@spatial-main")
 }
-current <- packageDescription("wwinference")$GithubRef
 
 #Load packages
 pacman::p_load(
@@ -41,8 +39,8 @@ compiled_model <- compile_model_upd(update_files = T)
 
 #Specify fit options
 fit_this <- get_mcmc_options(
-  iter_warmup = 250,
-  iter_sampling = 500,
+  iter_warmup = 250,           #Keep >100, 
+  iter_sampling = 500,         #Keep >100,
   seed = 1
 )
 
@@ -59,7 +57,9 @@ WA_nonspatial_run(
   #Population by catchment
   pop_data,
   #If listing individual sites, combine with ; (i.e. 2;3;4)
-  sites = "all",                         
+  sites = "all",
+  #Repeat subset
+  repeat_subset = "9",
   #WA population
   WA_population = 7786000,    
   #Number of days you are going to predict to
@@ -69,10 +69,10 @@ WA_nonspatial_run(
   #How many shuffles of data we want to repeat this on
   runs = 10,  
   #Savename modifier for file output
-  savename = "nonspatial_state",    
+  savename = "updated_nonspatial_state",    
   #Model fit options - specified above
   fit_options = fit_this,
-  #The compiled model
+  #The compiled model;
   compiled_model = compiled_model
 )
 

R/scripts/2_assess_nonspatial_substate.R

@@ -1,12 +1,10 @@
 #Install pacman if missing
 if(!require(pacman)) install.packages("pacman")
 
-#If you have installed a NON-MAIN branch version of the package, then re-install the main branch
-prior <- packageDescription("wwinference")$GithubRef
-if(packageDescription("wwinference")$GithubRef != "HEAD"){
-  devtools::install_github("CDCgov/ww-inference-model")
+#If you have installed the NON-SPATIAL version of the package, then re-install the spatial branch
+if(packageDescription("wwinference")$GithubRef != "spatial-main"){
+  devtools::install_github("CDCgov/ww-inference-model@spatial-main")
 }
-current <- packageDescription("wwinference")$GithubRef
 
 #Load packages
 pacman::p_load(

R/scripts/3_assess_spatial_state.R

@@ -2,11 +2,9 @@
 if(!require(pacman)) install.packages("pacman")
 
 #If you have installed the NON-SPATIAL version of the package, then re-install the spatial branch
-prior <- packageDescription("wwinference")$GithubRef
 if(packageDescription("wwinference")$GithubRef != "spatial-main"){
   devtools::install_github("CDCgov/ww-inference-model@spatial-main")
 }
-current <- packageDescription("wwinference")$GithubRef
 
 #Load packages
 pacman::p_load(
@@ -62,7 +60,9 @@ WA_spatial_run(
   #Population by catchment
   pop_data,
   #If listing individual sites, combine with ; (i.e. 2;3;4)
-  sites = "all",                         
+  sites = "all",               
+  #Repeat subset
+  repeat_subset = "9",
   #WA population
   WA_population = 7786000,    
   #Number of days you are going to predict to
@@ -72,7 +72,7 @@ WA_spatial_run(
   #How many shuffles of data we want to repeat this on
   runs = 10,  
   #Save name modifier for file output
-  savename = "spatial_state",    
+  savename = "spatial_parameter_update",    
   #Model fit options - specified above
   fit_options = fit_this,
   #The compiled model

output/model_runs/summary/nonspatial_realdata/key.csv

@@ -0,0 +1,11 @@
+"sites","forecast_horizon","calibration_time","runs","savename","random_dates"
+"all",28,90,10,"nonspatial_realdata",2023-01-08
+"all",28,90,10,"nonspatial_realdata",2023-04-17
+"all",28,90,10,"nonspatial_realdata",2023-03-10
+"all",28,90,10,"nonspatial_realdata",2023-04-12
+"all",28,90,10,"nonspatial_realdata",2022-12-14
+"all",28,90,10,"nonspatial_realdata",2022-11-15
+"all",28,90,10,"nonspatial_realdata",2022-12-22
+"all",28,90,10,"nonspatial_realdata",2023-01-25
+"all",28,90,10,"nonspatial_realdata",2022-11-22
+"all",28,90,10,"nonspatial_realdata",2023-02-15

output/model_runs/summary/spatial_parameter_update/key.csv

@@ -0,0 +1,11 @@
+"sites","forecast_horizon","calibration_time","runs","savename","spatial","random_dates"
+"all",28,90,10,"spatial_parameter_update",TRUE,2023-01-08
+"all",28,90,10,"spatial_parameter_update",TRUE,2023-04-17
+"all",28,90,10,"spatial_parameter_update",TRUE,2023-03-10
+"all",28,90,10,"spatial_parameter_update",TRUE,2023-04-12
+"all",28,90,10,"spatial_parameter_update",TRUE,2022-12-14
+"all",28,90,10,"spatial_parameter_update",TRUE,2022-11-15
+"all",28,90,10,"spatial_parameter_update",TRUE,2022-12-22
+"all",28,90,10,"spatial_parameter_update",TRUE,2023-01-25
+"all",28,90,10,"spatial_parameter_update",TRUE,2022-11-22
+"all",28,90,10,"spatial_parameter_update",TRUE,2023-02-15

output/model_runs/summary/spatial_state/key.csv

@@ -1,4 +1,4 @@
-"sites","forecast_horizon","calibration_time","repeats","savename","spatial","random_dates"
+"sites","forecast_horizon","calibration_time","runs","savename","spatial","random_dates"
 "all",28,90,10,"spatial_state",TRUE,2023-01-08
 "all",28,90,10,"spatial_state",TRUE,2023-04-17
 "all",28,90,10,"spatial_state",TRUE,2023-03-10