Revert "Testing slightly different viral kinetics curve"

This reverts commit 7ccaf24.
jameshay218 · May 20, 2021 · 8ad1e19 · 8ad1e19
1 parent 7ccaf24
commit 8ad1e19
Show file tree

Hide file tree

Showing 31 changed files with 55 additions and 89 deletions.
diff --git a/1.nursing_homes/1.nursing_home_analyses_combined.R b/1.nursing_homes/1.nursing_home_analyses_combined.R
@@ -27,7 +27,7 @@ setwd(main_wd)
 
 ## IMPORTANT - change this flag to TRUE if running the MCMC for the first time
 rerun_mcmc_seeirr <- FALSE
-rerun_mcmc_ct <- TRUE
+rerun_mcmc_ct <- FALSE
 
 ## Clean data and set up model parameters. Also some data plots.
 source("1.nursing_homes/1.nursing_home_analyses_headers.R")

diff --git a/5.bwh/5.fit_ma_gp.R b/5.bwh/5.fit_ma_gp.R
@@ -27,10 +27,10 @@ devtools::load_all(paste0(HOME_WD,"/virosolver"))
 index <- 1994
 set.seed(index)
 n_samp <- 1000
-runname <- "ma_gp_test"
+runname <- "ma_gp"
 #runname <- "ma_gp_free"
 run_version <- "gp" ##gp, seir or exp##
-rerun_mcmc <- TRUE
+rerun_mcmc <- FALSE
 
 ## CHANGE TO MAIN WD
 ## Important to set this to the full file path, as on L205 the foreach loop
@@ -231,7 +231,7 @@ for(ii in seq_along(samps)){
   tmp_pars <- get_index_pars(chain_comb, samps[ii])
   vl <- viral_load_func(tmp_pars,test_ages,FALSE)
   vl_trajs[ii,]  <- extraDistr::rgumbel(length(vl),vl,tmp_pars["obs_sd"])
-  detect_trajs[ii,] <- c(0,prop_detectable(test_ages[test_ages > 0], tmp_pars, vl[test_ages > 0]))
+  detect_trajs[ii,] <- c(0,prop_detectable(test_ages[test_ages > 0], tmp_pars, vl))
   trajs[ii,] <- pmax(smooth.spline(inc_func_use(tmp_pars,times))$y,0.0000001)
   #trajs[ii,] <- pmax(inc_func_use(get_index_pars(chain_comb, samps[ii]),times),0.0000001)
 }
@@ -391,7 +391,7 @@ p_inc <- ggplot(trajs_quants) +
   export_theme +
   ylab("Relative probability of infection") +
   xlab("Date") +
-  #scale_x_date(limits=as.Date(c("2020-03-01","2020-12-01")),breaks="1 month",expand=c(0,0)) +
+  scale_x_date(limits=as.Date(c("2020-03-01","2020-12-01")),breaks="1 month",expand=c(0,0)) +
   #coord_cartesian(ylim=c(-0.0001,0.005)) +
   scale_y_continuous(expand=c(0,0))  +
   labs(tag="C")

diff --git a/5.bwh/5.fit_ma_single_timepoints_exp.R b/5.bwh/5.fit_ma_single_timepoints_exp.R
@@ -27,7 +27,7 @@ devtools::load_all(paste0(HOME_WD,"lazymcmc"))
 ## Arguments for this run
 set.seed(1)
 n_samp <- 1000
-runname <- "ma_exp_test"
+runname <- "ma_exp"
 run_version <- "exp" ##gp, seir or exp##
 
 ## IMPORTANT - change this flag to TRUE if running the MCMC for the first time

diff --git a/5.bwh/5.fit_ma_single_timepoints_seir.R b/5.bwh/5.fit_ma_single_timepoints_seir.R
@@ -27,7 +27,7 @@ devtools::load_all(paste0(HOME_WD,"lazymcmc"))
 ## Arguments for this run
 set.seed(1)
 n_samp <- 1000
-runname <- "ma_seir_test"
+runname <- "ma_seir_exposed_seed"
 run_version <- "seir" ##gp, seir or exp##
 
 ## IMPORTANT - change this flag to TRUE if running the MCMC for the first time

diff --git a/5.bwh/Figure4.R b/5.bwh/Figure4.R
@@ -39,15 +39,15 @@ obs_times <- c(35, 42, 49, 56, 63, 70, 77, 84, 91, 98, 105, 112, 119, 126,
 seir_timeshift <- 36
 
 ## Run names for the various MCMC chains
-runname_seir <- "ma_seir_test" ## Constrained prior on seed time
+runname_seir <- "ma_seir_constrained" ## Constrained prior on seed time
 runname_seir_prior <- "ma_seir_prior_constrained"
 runname_seir_prior_broad <- "ma_seir_rerun_prior"
 runname_seir_alt <- "ma_seir_rerun" ## Less constrained prior on seed time
-runname_exp <- "ma_exp_test"
+runname_exp <- "ma_exp"
 
 ## Where the MCMC chains are stored
 chainwd <- paste0("~/Documents/GitHub/virosolver_paper/mcmc_chains/5.real_ma_single_timepoint/")
-chainwd_gp <- paste0("~/Documents/GitHub/virosolver_paper/mcmc_chains/4.real_ma_ct/ma_gp_test/")
+chainwd_gp <- paste0("~/Documents/GitHub/virosolver_paper/mcmc_chains/4.real_ma_ct/ma_gp/")
 
 ## MCMC parameters for Ct model fits
 mcmcPars_ct_seir <- c("adaptive_period"=30000)

diff --git a/figures/Figure2.pdf b/figures/Figure2.pdf
diff --git a/figures/Figure2.png b/figures/Figure2.png
diff --git a/figures/ma_gp_test.pdf b/figures/ma_gp_test.pdf
diff --git a/figures/supplement/FigureS4_a.pdf b/figures/supplement/FigureS4_a.pdf
diff --git a/figures/supplement/FigureS4_a.png b/figures/supplement/FigureS4_a.png
diff --git a/figures/supplement/FigureS4_b.pdf b/figures/supplement/FigureS4_b.pdf
diff --git a/figures/supplement/FigureS4_b.png b/figures/supplement/FigureS4_b.png
diff --git a/figures/supplement/biobot.png b/figures/supplement/biobot.png
diff --git a/figures/supplement/nh_data.pdf b/figures/supplement/nh_data.pdf
diff --git a/figures/supplement/nh_posteriors.pdf b/figures/supplement/nh_posteriors.pdf
diff --git a/figures/supplement/nh_posteriors.png b/figures/supplement/nh_posteriors.png
diff --git a/figures/supplement/seeirr_chains.pdf b/figures/supplement/seeirr_chains.pdf
diff --git a/figures/supplement/viral_kinetics_bwh.pdf b/figures/supplement/viral_kinetics_bwh.pdf
diff --git a/figures/supplement/viral_kinetics_bwh.png b/figures/supplement/viral_kinetics_bwh.png
diff --git a/figures/supplement/viral_kinetics_nh.pdf b/figures/supplement/viral_kinetics_nh.pdf
diff --git a/figures/supplement/viral_kinetics_nh.png b/figures/supplement/viral_kinetics_nh.png
diff --git a/pars/massachusetts/partab_exp_model.csv b/pars/massachusetts/partab_exp_model.csv
@@ -1,18 +1,18 @@
 values,names,fixed,lower_bound,upper_bound,steps,lower_start,upper_start
 0.1,beta,0,-0.5,0.5,0.1,-0.2,0.2
 1,overall_prob,0,0,1,0.1,0,1
-1,tshift,1,0,3,0.1,0,10
+0,tshift,1,0,3,0.1,0,10
 5,desired_mode,1,0,7,0.1,0,10
-21.56194409,viral_peak,0,0,40,0.1,15,25
-6.195697774,obs_sd,0,0,25,0.1,1,10
-0.636291677,sd_mod,1,0,1,0.1,0.4,0.6
+19.73598747,viral_peak,0,0,40,0.1,15,25
+5,obs_sd,0,0,25,0.1,1,10
+0.788828786,sd_mod,1,0,1,0.1,0.4,0.6
 14,sd_mod_wane,1,0,14,0.1,0,14
 40,true_0,1,40,100,0.1,40,100
 40,intercept,1,35,100,0.1,35,100
 3,LOD,1,0,10,0.1,0,10
 5,incu,1,0,10,0.1,0,10
-13.72235562,t_switch,0,0,30,0.1,10,30
+13.29211948,t_switch,0,0,30,0.1,10,30
 38,level_switch,0,0,40,0.1,33,40
 1000,wane_rate2,1,0,10000,0.1,10,50
-0.123203205,prob_detect,0,0,1,0.1,0.01,0.1
+0.103299323,prob_detect,0,0,1,0.1,0.01,0.1
 1,t_unit,1,0,1,0.1,0,1
diff --git a/pars/massachusetts/partab_exp_pos_model.csv b/pars/massachusetts/partab_exp_pos_model.csv
@@ -1,18 +1,18 @@
 values,names,fixed,lower_bound,upper_bound,steps,lower_start,upper_start
 -0.1,beta,0,-0.5,0.5,0.1,-0.2,0.2
 1,overall_prob,1,0,1,0.1,0,1
-1,tshift,1,0,3,0.1,0,10
+0,tshift,1,0,3,0.1,0,10
 5,desired_mode,1,0,7,0.1,0,10
-21.56194409,viral_peak,0,0,40,0.1,15,25
-6.195697774,obs_sd,0,0,25,0.1,1,10
-0.636291677,sd_mod,1,0,1,0.1,0.4,0.6
+19.73598747,viral_peak,0,0,40,0.1,15,25
+5,obs_sd,0,0,25,0.1,1,10
+0.788828786,sd_mod,1,0,1,0.1,0.4,0.6
 14,sd_mod_wane,1,0,14,0.1,0,14
 40,true_0,1,40,100,0.1,40,100
 40,intercept,1,35,100,0.1,35,100
 3,LOD,1,0,10,0.1,0,10
 5,incu,1,0,10,0.1,0,10
-13.72235562,t_switch,0,0,30,0.1,10,30
+13.29211948,t_switch,0,0,30,0.1,10,30
 38,level_switch,0,0,40,0.1,33,40
 1000,wane_rate2,1,0,10000,0.1,10,50
-0.123203205,prob_detect,0,0,1,0.1,0.01,0.1
+0.103299323,prob_detect,0,0,1,0.1,0.01,0.1
 1,t_unit,1,0,1,0.1,0,1
diff --git a/pars/massachusetts/partab_gp_model.csv b/pars/massachusetts/partab_gp_model.csv
@@ -1,19 +1,19 @@
 values,names,fixed,lower_bound,upper_bound,steps,lower_start,upper_start
 0.5,overall_prob,0,0,1,0.1,0,1
-1,tshift,1,0,3,0.1,0,10
+0,tshift,1,0,3,0.1,0,10
 5,desired_mode,1,0,7,0.1,0,10
-21.56194409,viral_peak,0,0,40,0.1,15,25
-6.195697774,obs_sd,0,0,25,0.1,1,10
-0.636291677,sd_mod,1,0,1,0.1,0.4,0.6
+19.73598747,viral_peak,0,0,40,0.1,15,25
+5,obs_sd,0,0,25,0.1,1,10
+0.788828786,sd_mod,1,0,1,0.1,0.4,0.6
 14,sd_mod_wane,1,0,14,0.1,0,14
 40,true_0,1,40,100,0.1,40,100
 40,intercept,1,35,100,0.1,35,100
 3,LOD,1,0,10,0.1,0,10
 5,incu,1,0,10,0.1,0,10
-13.72235562,t_switch,0,0,30,0.1,10,30
+13.29211948,t_switch,0,0,30,0.1,10,30
 38,level_switch,0,0,40,0.1,33,40
 1000,wane_rate2,1,0,10000,0.1,10,50
-0.123203205,prob_detect,0,0,1,0.1,0.01,0.1
+0.103299323,prob_detect,0,0,1,0.1,0.01,0.1
 1,t_unit,1,0,1,0.1,0,1
 0,t0,1,0,100,0.1,0,50
 1.5,nu,1,0,100,0.1,2,5

diff --git a/pars/massachusetts/partab_seir_model.csv b/pars/massachusetts/partab_seir_model.csv
@@ -1,18 +1,18 @@
 values,names,fixed,lower_bound,upper_bound,steps,lower_start,upper_start
-1,tshift,1,0,3,0.1,0,10
+0,tshift,1,0,3,0.1,0,10
 5,desired_mode,1,0,7,0.1,0,10
-21.56194409,viral_peak,0,0,40,0.1,15,25
-6.195697774,obs_sd,0,0,25,0.1,1,10
-0.636291677,sd_mod,1,0,1,0.1,0.4,0.6
+19.73598747,viral_peak,0,0,40,0.1,15,25
+5,obs_sd,0,0,25,0.1,1,10
+0.788828786,sd_mod,1,0,1,0.1,0.4,0.6
 14,sd_mod_wane,1,0,14,0.1,0,14
 40,true_0,1,40,100,0.1,40,100
 40,intercept,1,35,100,0.1,35,100
 3,LOD,1,0,10,0.1,0,10
 5,incu,1,0,10,0.1,0,10
-13.72235562,t_switch,0,0,30,0.1,10,30
+13.29211948,t_switch,0,0,30,0.1,10,30
 38,level_switch,0,0,40,0.1,33,40
 1000,wane_rate2,1,0,10000,0.1,10,50
-0.123203205,prob_detect,0,0,1,0.1,0.01,0.1
+0.103299323,prob_detect,0,0,1,0.1,0.01,0.1
 1,t_unit,1,0,1,0.1,0,1
 2.5,R0,0,1,10,0.1,1.5,3
 4,infectious,0,0,25,0.1,5,10

diff --git a/pars/nursing_homes/partab_exp_model.csv b/pars/nursing_homes/partab_exp_model.csv
@@ -1,18 +1,18 @@
 values,names,fixed,lower_bound,upper_bound,steps,lower_start,upper_start
 0.1,beta,0,-0.5,0.5,0.1,-0.2,0.2
 1,overall_prob,0,0,1,0.1,0,1
-1,tshift,1,0,3,0.1,0,10
+0,tshift,1,0,3,0.1,0,10
 5,desired_mode,1,0,7,0.1,0,10
-20.76267783,viral_peak,0,0,40,0.1,15,25
-5.216597448,obs_sd,0,0,25,0.1,1,10
-0.390277615,sd_mod,1,0,1,0.1,0.4,0.6
+20.66182195,viral_peak,0,0,40,0.1,15,25
+5.150558612,obs_sd,0,0,25,0.1,1,10
+0.395291595,sd_mod,1,0,1,0.1,0.4,0.6
 14,sd_mod_wane,1,0,14,0.1,0,14
 40,true_0,1,40,100,0.1,40,100
 40,intercept,1,35,100,0.1,35,100
 3,LOD,1,0,10,0.1,0,10
 5,incu,1,0,10,0.1,0,10
-8.39345933,t_switch,0,0,30,0.1,10,30
+9.383013136,t_switch,0,0,30,0.1,10,30
 33,level_switch,0,0,40,0.1,33,40
 1000,wane_rate2,1,0,10000,0.1,10,50
-0.109859503,prob_detect,0,0,1,0.1,0.01,0.1
+0.109940519,prob_detect,0,0,1,0.1,0.01,0.1
 1,t_unit,1,0,1,0.1,0,1
diff --git a/pars/nursing_homes/partab_exp_pos_model.csv b/pars/nursing_homes/partab_exp_pos_model.csv
@@ -1,18 +1,18 @@
 values,names,fixed,lower_bound,upper_bound,steps,lower_start,upper_start
 -0.1,beta,0,-0.5,0.5,0.1,-0.2,0.2
 1,overall_prob,1,0,1,0.1,0,1
-1,tshift,1,0,3,0.1,0,10
+0,tshift,1,0,3,0.1,0,10
 5,desired_mode,1,0,7,0.1,0,10
-20.76267783,viral_peak,0,0,40,0.1,15,25
-5.216597448,obs_sd,0,0,25,0.1,1,10
-0.390277615,sd_mod,1,0,1,0.1,0.4,0.6
+20.66182195,viral_peak,0,0,40,0.1,15,25
+5.150558612,obs_sd,0,0,25,0.1,1,10
+0.395291595,sd_mod,1,0,1,0.1,0.4,0.6
 14,sd_mod_wane,1,0,14,0.1,0,14
 40,true_0,1,40,100,0.1,40,100
 40,intercept,1,35,100,0.1,35,100
 3,LOD,1,0,10,0.1,0,10
 5,incu,1,0,10,0.1,0,10
-8.39345933,t_switch,0,0,30,0.1,10,30
+9.383013136,t_switch,0,0,30,0.1,10,30
 33,level_switch,0,0,40,0.1,33,40
 1000,wane_rate2,1,0,10000,0.1,10,50
-0.109859503,prob_detect,0,0,1,0.1,0.01,0.1
+0.109940519,prob_detect,0,0,1,0.1,0.01,0.1
 1,t_unit,1,0,1,0.1,0,1
diff --git a/pars/nursing_homes/partab_seir_model.csv b/pars/nursing_homes/partab_seir_model.csv
@@ -1,18 +1,18 @@
 values,names,fixed,lower_bound,upper_bound,steps,lower_start,upper_start
-1,tshift,1,0,3,0.1,0,10
+0,tshift,1,0,3,0.1,0,10
 5,desired_mode,1,0,7,0.1,0,10
-20.76267783,viral_peak,0,0,40,0.1,15,25
-5.216597448,obs_sd,0,0,25,0.1,1,10
-0.390277615,sd_mod,1,0,1,0.1,0.4,0.6
+20.66182195,viral_peak,0,0,40,0.1,15,25
+5.150558612,obs_sd,0,0,25,0.1,1,10
+0.395291595,sd_mod,1,0,1,0.1,0.4,0.6
 14,sd_mod_wane,1,0,14,0.1,0,14
 40,true_0,1,40,100,0.1,40,100
 40,intercept,1,35,100,0.1,35,100
 3,LOD,1,0,10,0.1,0,10
 5,incu,1,0,10,0.1,0,10
-8.39345933,t_switch,0,0,30,0.1,10,30
+9.383013136,t_switch,0,0,30,0.1,10,30
 33,level_switch,0,0,40,0.1,33,40
 1000,wane_rate2,1,0,10000,0.1,10,50
-0.109859503,prob_detect,0,0,1,0.1,0.01,0.1
+0.109940519,prob_detect,0,0,1,0.1,0.01,0.1
 1,t_unit,1,0,1,0.1,0,1
 2.5,R0,0,1,10,0.1,1.5,3
 4,infectious,0,0,25,0.1,5,10

diff --git a/pars/partab_fitted_bwh.csv b/pars/partab_fitted_bwh.csv
diff --git a/pars/partab_fitted_nh.csv b/pars/partab_fitted_nh.csv
diff --git a/scripts/a.fit_viral_load_model.R b/scripts/a.fit_viral_load_model.R
@@ -64,7 +64,7 @@ targeted_dat <- tibble(age=ages_observed, prob=desired_probs)
 
 pars_hinge2 <- parTab$values
 names(pars_hinge2) <- parTab$names
-pars_hinge2["tshift"] <- 1
+
 
 ######################################################
 ## COST FUNCTION
@@ -327,11 +327,11 @@ if(SAVE_PLOTS){
     ggsave("figures/supplement/viral_kinetics_bwh.pdf",p_main,height=6,width=8)
     ggsave("figures/supplement/viral_kinetics_bwh.png",p_main,height=6,width=8)
     parTab$values <- pars1
-    write_csv(parTab,"pars/partab_fitted_bwh.csv")
+    #write_csv(parTab,"pars/partab_fitted_bwh.csv")
   } else {
     ggsave("figures/supplement/viral_kinetics_nh.pdf",p_main,height=6,width=8)
     ggsave("figures/supplement/viral_kinetics_nh.png",p_main,height=6,width=8)
     parTab$values <- pars1
-    write_csv(parTab,"pars/partab_fitted_nh.csv")
+    #write_csv(parTab,"pars/partab_fitted_nh.csv")
   }
 }