Updates for resubmission

code/fit_posteriors.R

@@ -20,7 +20,8 @@ ct_fit <- sampling(ct_model,
 		wr_midpoint=prior_pars$wr_midpoint,
 		sigma_prior=prior_pars$sigma_prior,
 		lambda=prior_pars$lambda,
-		fpmean=prior_pars$fpmean), 
+		fpmean=prior_pars$fpmean,
+		priorsd=prior_pars$priorsd), 
 	iter=1000, chains=4) # iter=2000
 # , control = list(adapt_delta=0.85)
 # , control = list(adapt_delta=0.99)

code/fit_posteriors.stan

@@ -15,8 +15,8 @@ data{
   int<lower=0> n_id;       // Number of infections
   real<lower=0> lod;        // Limit of detection
   int<lower=0> id[N];       // Vector marking which datum belongs to which id
-  int<lower=0> category[n_id];  // Immune categories
-  int<lower=0> max_category;     // Highest immune category
+  int<lower=0> category[n_id];  // Analysis categories
+  int<lower=0> max_category;     // Highest analysis category
   int<lower=0> n_adj;        // Number of adjustment variables 
   int<lower=0> adjust[n_id,n_adj];  // Adjustment categories
   int<lower=0> max_adjust;     // Highest adjustment category, across all adjustment variables
@@ -29,6 +29,7 @@ data{
   real sigma_prior[2];      // Prior observation noise Cauchy scale
   real<lower=0, upper=1> lambda; // Mixture probability (~1-sensitivity)
   real<lower=0> fpmean;     // False positive mean Ct
+  real<lower=0> priorsd;    // Standard deviation for the lognormal priors
 }
 
 
@@ -128,19 +129,16 @@ transformed parameters{
 
     dp[i] = exp(log_dp_mean + 
       log_dpadj[category[i]] + 
-      // log_dpadj_adjust[adjust[i]] + 
       dpadj_add + 
       log_dp_sd*dp_raw[i])*dp_midpoint;
 
     wp[i] = exp(log_wp_mean + 
       log_wpadj[category[i]] + 
-      // log_wpadj_adjust[adjust[i]] + 
       wpadj_add + 
       log_wp_sd*wp_raw[i])*wp_midpoint;
 
     wr[i] = exp(log_wr_mean + 
       log_wradj[category[i]] + 
-      // log_wradj_adjust[adjust[i]] + 
       wradj_add + 
       log_wr_sd*wr_raw[i])*wr_midpoint;
 
@@ -160,28 +158,27 @@ model{
 
   tp ~ normal(tp_prior[1], tp_prior[2]);
 
-  log_dp_mean ~ normal(0, 0.25); // hierarchical mean
-  log_dp_sd ~ normal(0, 0.25) T[0,]; // sd of the individual-level draws
-  log_dpadj_cat_raw ~ normal(0, 0.25); 
+  log_dp_mean ~ normal(0, priorsd); // hierarchical mean
+  log_dp_sd ~ normal(0, priorsd) T[0,]; // sd of the individual-level draws
+  log_dpadj_cat_raw ~ normal(0, priorsd); 
   for(indexA in 1:(max_adjust-1)){
-    log_dpadj_adjust_raw[indexA] ~ normal(0, 0.25);   
+    log_dpadj_adjust_raw[indexA] ~ normal(0, priorsd);
   }
-  // log_dpadj_adjust_raw ~ normal(0, 0.25); 
   dp_raw ~ std_normal(); // for generating individual-level guesses 
 
-  log_wp_mean ~ normal(0, 0.25); // hierarchical mean
-  log_wp_sd ~ normal(0, 0.25) T[0,]; // sd of the individual-level draws
-  log_wpadj_cat_raw ~ normal(0, 0.25); 
+  log_wp_mean ~ normal(0, priorsd); // hierarchical mean
+  log_wp_sd ~ normal(0, priorsd) T[0,]; // sd of the individual-level draws
+  log_wpadj_cat_raw ~ normal(0, priorsd); 
   for(indexA in 1:(max_adjust-1)){
-    log_wpadj_adjust_raw[indexA] ~ normal(0, 0.25); 
+    log_wpadj_adjust_raw[indexA] ~ normal(0, priorsd); 
   }
   wp_raw ~ std_normal(); // for generating individual-level guesses
 
-  log_wr_mean ~ normal(0, 0.25); // hierarchical mean
-  log_wr_sd ~ normal(0, 0.25) T[0,]; // sd of the individual-level draws
-  log_wradj_cat_raw ~ normal(0, 0.25); 
+  log_wr_mean ~ normal(0, priorsd); // hierarchical mean
+  log_wr_sd ~ normal(0, priorsd) T[0,]; // sd of the individual-level draws
+  log_wradj_cat_raw ~ normal(0, priorsd); 
   for(indexA in 1:(max_adjust-1)){
-    log_wradj_adjust_raw[indexA] ~ normal(0, 0.25); 
+    log_wradj_adjust_raw[indexA] ~ normal(0, priorsd); 
   }
   wr_raw ~ std_normal(); // for generating individual-level guesses
 

code/fit_posteriors_preamble.R

@@ -77,6 +77,7 @@ prior_pars <- list(
 	category=catlist,
 	max_category=max(catlist),
 	adjust=adjmat,
-	max_adjust=max(max(adjmat))
+	max_adjust=max(max(adjmat)),
+	priorsd=run_pars$priorsd
 )	
 

code/make_inf_table.R

@@ -0,0 +1,168 @@
+# FIRST INFECTION =============================================================
+
+# By variant:
+ct_dat_refined %>% 
+	group_by(InfectionEvent) %>% 
+	slice(1) %>% 
+	filter(InfNum==1) %>%
+	mutate(LineageBroad=case_when(
+		LineageBroad=="BA.1"~"BA.1/BA.2", 
+		LineageBroad=="BA.2"~"BA.1/BA.2", 
+		LineageBroad=="BA.4"~"BA.4/BA.5", 
+		LineageBroad=="BA.5"~"BA.4/BA.5",
+		LineageBroad=="Other"~"Other/Unspecified",
+		LineageBroad=="None"~"Other/Unspecified",
+		TRUE~LineageBroad 
+		)) %>% 
+	mutate(LineageBroad=factor(LineageBroad, levels=c("Alpha","Delta","BA.1/BA.2","BA.4/BA.5","Other/Unspecified"))) %>% 
+	group_by(LineageBroad) %>% 
+	summarise(N=n()) %>% 
+	ungroup() %>% 
+	mutate(pct=round(100*N/sum(N),1)) %>% 
+	print(n=Inf)
+
+# By variant + vax status:
+ct_dat_refined %>% 
+	group_by(InfectionEvent) %>% 
+	slice(1) %>% 
+	filter(InfNum==1) %>%
+	mutate(LineageBroad=case_when(
+		LineageBroad=="BA.1"~"BA.1/BA.2", 
+		LineageBroad=="BA.2"~"BA.1/BA.2", 
+		LineageBroad=="BA.4"~"BA.4/BA.5", 
+		LineageBroad=="BA.5"~"BA.4/BA.5",
+		LineageBroad=="Other"~"Other/Unspecified",
+		LineageBroad=="None"~"Other/Unspecified",
+		TRUE~LineageBroad 
+		)) %>% 
+	mutate(LineageBroad=factor(LineageBroad, levels=c("Alpha","Delta","BA.1/BA.2","BA.4/BA.5","Other/Unspecified"))) %>% 
+	mutate(VaccinationStatus=case_when(is.na(VaccinationStatus)~"Not Reported",TRUE~VaccinationStatus)) %>% 
+	mutate(VaccinationStatus=factor(VaccinationStatus, levels=c("Not Vaccinated","Fully Vaccinated","Not Reported"))) %>% 
+	group_by(LineageBroad, VaccinationStatus) %>% 
+	summarise(N=n()) %>% 
+	ungroup() %>% 
+	mutate(pct=round(100*N/sum(N),1)) %>% 
+	print(n=Inf)
+
+
+# By variant + boost status:
+ct_dat_refined %>% 
+	group_by(InfectionEvent) %>% 
+	slice(1) %>% 
+	filter(InfNum==1) %>%
+	mutate(LineageBroad=case_when(
+		LineageBroad=="BA.1"~"BA.1/BA.2", 
+		LineageBroad=="BA.2"~"BA.1/BA.2", 
+		LineageBroad=="BA.4"~"BA.4/BA.5", 
+		LineageBroad=="BA.5"~"BA.4/BA.5",
+		LineageBroad=="Other"~"Other/Unspecified",
+		LineageBroad=="None"~"Other/Unspecified",
+		TRUE~LineageBroad 
+		)) %>% 
+	mutate(LineageBroad=factor(LineageBroad, levels=c("Alpha","Delta","BA.1/BA.2","BA.4/BA.5","Other/Unspecified"))) %>% 
+	mutate(BoosterStatus=case_when(is.na(BoosterStatus)~"Not Reported",TRUE~BoosterStatus)) %>% 
+	mutate(BoosterStatus=factor(BoosterStatus, levels=c("Not Boosted","Boosted","Not Reported"))) %>% 
+	group_by(LineageBroad, BoosterStatus) %>% 
+	summarise(N=n()) %>% 
+	ungroup() %>% 
+	mutate(pct=round(100*N/sum(N),1)) %>% 
+	print(n=Inf)
+
+
+# By age group:
+ct_dat_refined %>% 
+	group_by(InfectionEvent) %>% 
+	slice(1) %>% 
+	filter(InfNum==1) %>%
+	group_by(AgeGrp) %>% 
+	summarise(N=n()) %>% 
+	ungroup() %>% 
+	mutate(pct=round(100*N/sum(N),1)) %>% 
+	print(n=Inf)
+
+
+
+
+
+# SECOND INFECTION ===========================================================
+
+# By variant:
+ct_dat_refined %>% 
+	group_by(InfectionEvent) %>% 
+	slice(1) %>% 
+	filter(InfNum==2) %>%
+	mutate(LineageBroad=case_when(
+		LineageBroad=="BA.1"~"BA.1/BA.2", 
+		LineageBroad=="BA.2"~"BA.1/BA.2", 
+		LineageBroad=="BA.4"~"BA.4/BA.5", 
+		LineageBroad=="BA.5"~"BA.4/BA.5",
+		LineageBroad=="Other"~"Other/Unspecified",
+		LineageBroad=="None"~"Other/Unspecified",
+		TRUE~LineageBroad 
+		)) %>% 
+	mutate(LineageBroad=factor(LineageBroad, levels=c("Alpha","Delta","BA.1/BA.2","BA.4/BA.5","Other/Unspecified"))) %>% 
+	group_by(LineageBroad) %>% 
+	summarise(N=n()) %>% 
+	ungroup() %>% 
+	mutate(pct=round(100*N/sum(N),1)) %>% 
+	print(n=Inf)
+
+# By variant + vax status:
+ct_dat_refined %>% 
+	group_by(InfectionEvent) %>% 
+	slice(1) %>% 
+	filter(InfNum==2) %>%
+	mutate(LineageBroad=case_when(
+		LineageBroad=="BA.1"~"BA.1/BA.2", 
+		LineageBroad=="BA.2"~"BA.1/BA.2", 
+		LineageBroad=="BA.4"~"BA.4/BA.5", 
+		LineageBroad=="BA.5"~"BA.4/BA.5",
+		LineageBroad=="Other"~"Other/Unspecified",
+		LineageBroad=="None"~"Other/Unspecified",
+		TRUE~LineageBroad 
+		)) %>% 
+	mutate(LineageBroad=factor(LineageBroad, levels=c("Alpha","Delta","BA.1/BA.2","BA.4/BA.5","Other/Unspecified"))) %>% 
+	mutate(VaccinationStatus=case_when(is.na(VaccinationStatus)~"Not Reported",TRUE~VaccinationStatus)) %>% 
+	mutate(VaccinationStatus=factor(VaccinationStatus, levels=c("Not Vaccinated","Fully Vaccinated","Not Reported"))) %>% 
+	group_by(LineageBroad, VaccinationStatus) %>% 
+	summarise(N=n()) %>% 
+	ungroup() %>% 
+	mutate(pct=round(100*N/sum(N),1)) %>% 
+	print(n=Inf)
+
+
+# By variant + boost status:
+ct_dat_refined %>% 
+	group_by(InfectionEvent) %>% 
+	slice(1) %>% 
+	filter(InfNum==2) %>%
+	mutate(LineageBroad=case_when(
+		LineageBroad=="BA.1"~"BA.1/BA.2", 
+		LineageBroad=="BA.2"~"BA.1/BA.2", 
+		LineageBroad=="BA.4"~"BA.4/BA.5", 
+		LineageBroad=="BA.5"~"BA.4/BA.5",
+		LineageBroad=="Other"~"Other/Unspecified",
+		LineageBroad=="None"~"Other/Unspecified",
+		TRUE~LineageBroad 
+		)) %>% 
+	mutate(LineageBroad=factor(LineageBroad, levels=c("Alpha","Delta","BA.1/BA.2","BA.4/BA.5","Other/Unspecified"))) %>% 
+	mutate(BoosterStatus=case_when(is.na(BoosterStatus)~"Not Reported",TRUE~BoosterStatus)) %>% 
+	mutate(BoosterStatus=factor(BoosterStatus, levels=c("Not Boosted","Boosted","Not Reported"))) %>% 
+	group_by(LineageBroad, BoosterStatus) %>% 
+	summarise(N=n()) %>% 
+	ungroup() %>% 
+	mutate(pct=round(100*N/sum(N),1)) %>% 
+	print(n=Inf)
+
+
+# By age group:
+ct_dat_refined %>% 
+	group_by(InfectionEvent) %>% 
+	slice(1) %>% 
+	filter(InfNum==2) %>%
+	group_by(AgeGrp) %>% 
+	summarise(N=n()) %>% 
+	ungroup() %>% 
+	mutate(pct=round(100*N/sum(N),1)) %>% 
+	print(n=Inf)
+

code/run_analysis.R

@@ -32,7 +32,7 @@ source("code/set_run_pars.R")
 # Run the analysis
 # =============================================================================
 
-for(run_pars_index in 11:11){
+for(run_pars_index in 16:17){
 
 	run_pars <- run_pars_list[[run_pars_index]]
 
@@ -55,7 +55,7 @@ for(run_pars_index in 11:11){
 # Individual-level comparisons (Spearman correlations)
 # =============================================================================
 
-for(run_pars_index in 10:10){
+for(run_pars_index in 14:14){
 
 	run_pars <- run_pars_list[[run_pars_index]]
 
@@ -72,7 +72,7 @@ for(run_pars_index in 10:10){
 # Post-hoc figure generation
 # =============================================================================
 
-for(run_pars_index in 10:10){
+for(run_pars_index in 14:14){
 
 	run_pars <- run_pars_list[[run_pars_index]]