99_psa_runs.Rmd

---
title: "PSA Runs"
author: "Carol Liu"
date: "6/1/2020"
output:
  html_document:
    df_print: paged
---


```{r}
memory.limit(size=50000000)

knitr::opts_chunk$set(echo = TRUE)
# Load dependencies, functions and parameters
source("99_dependencies.R")
source("99_model_func.R")
source("99_parm_init_control.R")
source("99_psa_optimizedistr.R")
source("99_psa_parm.R")   #Note this overwrites initial parameters from parm_init_control
source("99_psa_plot.R")
options(scipen=999)
options(digits=4)

```

## Base case
```{r base_results}
## Scenario with seven day screening, two day testing, contact tracing as variable parameter
param <- param.dcm(latent = latent.int,
                   infectious = infectious.int,
                   isolation = isolation,
                   beta_student_to_student = beta_student_to_student.int,
                   beta_on_to_on = beta_on_to_on.int,
                   beta_saf = beta_saf.int,
                   community = community.int,
                   p_asympt_stu = p_asympt_stu.int,
                   p_asympt_saf = p_asympt_saf.int,
                   p_hosp_stu = p_hosp_stu.int,
                   p_hosp_saf = p_hosp_saf.int,
                   p_death_stu = p_death_stu.int,
                   p_death_saf = p_death_saf.int,
                   contacts = contacts.int,
                   p_contacts_reached = p_contacts_reached.int,
                   ili = ili.int,
                   sensitivity = sensitivity.int,
                   eff_npi = eff_npi.int,
                   testing = 0,
                   screening = 0,
                   screening_on=0)

mod_base <- dcm(param, init, control)
mod_base <- mutate_epi(mod_base, I_stu = Isym_on + Isym_off,
                  Icum_stu = Icum_on + Icum_off,
                  P_stu = P_on + P_off,
                  Pcum_stu = Pcum_on + Pcum_off,
                  Qcum_stu = Qcum_on + Qcum_off,
                  Hcum_stu = Hcum_on + Hcum_off,
                  Dcum_stu = Dcum_on + Dcum_off)


```


```{r}
#' ## Four day test interval
## ----four_test_results---------------------------------------------------
## Scenario with seven day screening, two day testing, contact tracing as variable parameter
param <- param.dcm(latent = latent.int,
                   infectious = infectious.int,
                   isolation = isolation,
                   beta_student_to_student = beta_student_to_student.int,
                   beta_on_to_on = beta_on_to_on.int,
                   beta_saf = beta_saf.int,
                   community = community.int,
                   p_asympt_stu = p_asympt_stu.int,
                   p_asympt_saf = p_asympt_saf.int,
                   p_hosp_stu = p_hosp_stu.int,
                   p_hosp_saf = p_hosp_saf.int,
                   p_death_stu = p_death_stu.int,
                   p_death_saf = p_death_saf.int,
                   contacts = contacts.int,
                   p_contacts_reached = p_contacts_reached.int,
                   ili = ili.int,
                   sensitivity = sensitivity.int,
                   eff_npi = eff_npi.int,
                   testing = 1/4,
                   screening = 0,
                   screening_on =0)

mod_4test <- dcm(param, init, control)
mod_4test <- mutate_epi(mod_4test, I_stu = Isym_on + Isym_off,
                  Icum_stu = Icum_on + Icum_off,
                  P_stu = P_on + P_off,
                  Pcum_stu = Pcum_on + Pcum_off,
                  Qcum_stu = Qcum_on + Qcum_off,
                  Hcum_stu = Hcum_on + Hcum_off,
                  Dcum_stu = Dcum_on + Dcum_off)
```

```{r}
#' 
#' ## Thirty day screen
## ----thirty_screen_results-----------------------------------------------
## Scenario with seven day screening, two day testing, contact tracing as variable parameter
param <- param.dcm(latent = latent.int,
                   infectious = infectious.int,
                   isolation = isolation,
                   beta_student_to_student = beta_student_to_student.int,
                   beta_on_to_on = beta_on_to_on.int,
                   beta_saf = beta_saf.int,
                   community = community.int,
                   p_asympt_stu = p_asympt_stu.int,
                   p_asympt_saf = p_asympt_saf.int,
                   p_hosp_stu = p_hosp_stu.int,
                   p_hosp_saf = p_hosp_saf.int,
                   p_death_stu = p_death_stu.int,
                   p_death_saf = p_death_saf.int,
                   contacts = contacts.int,
                   p_contacts_reached = p_contacts_reached.int,
                   ili = ili.int,
                   sensitivity = sensitivity.int,
                   eff_npi = eff_npi.int,
                   testing = 0,
                   screening = 1/30,
                   screening_on=1/30)

mod_30screen <- dcm(param, init, control)
mod_30screen <- mutate_epi(mod_30screen, I_stu = Isym_on + Isym_off,
                  Icum_stu = Icum_on + Icum_off,
                  P_stu = P_on + P_off,
                  Pcum_stu = Pcum_on + Pcum_off,
                  Qcum_stu = Qcum_on + Qcum_off,
                  Hcum_stu = Hcum_on + Hcum_off,
                  Dcum_stu = Dcum_on + Dcum_off)


#' 
```

```{r}
#' ## Combination intervention
## ----combo_results-------------------------------------------------------
## Scenario with seven day screening, two day testing, contact tracing as variable parameter
param <- param.dcm(latent = latent.int,
                   infectious = infectious.int,
                   isolation = isolation,
                   beta_student_to_student = beta_student_to_student.int,
                   beta_on_to_on = beta_on_to_on.int,
                   beta_saf = beta_saf.int,
                   community = community.int,
                   p_asympt_stu = p_asympt_stu.int,
                   p_asympt_saf = p_asympt_saf.int,
                   p_hosp_stu = p_hosp_stu.int,
                   p_hosp_saf = p_hosp_saf.int,
                   p_death_stu = p_death_stu.int,
                   p_death_saf = p_death_saf.int,
                   contacts = contacts.int,
                   p_contacts_reached = p_contacts_reached.int,
                   ili = ili.int,
                   sensitivity = sensitivity.int,
                   eff_npi = eff_npi.int,
                   testing = 1/4,
                   screening = 1/30,
                   screening_on=1/30)

mod_combo <- dcm(param, init, control)
mod_combo <- mutate_epi(mod_combo, I_stu = Isym_on + Isym_off,
                  Icum_stu = Icum_on + Icum_off,
                  P_stu = P_on + P_off,
                  Pcum_stu = Pcum_on + Pcum_off,
                  Qcum_stu = Qcum_on + Qcum_off,
                  Hcum_stu = Hcum_on + Hcum_off,
                  Dcum_stu = Dcum_on + Dcum_off)


#' 
```

```{r}
#' 
#' ## Plots of all scenarios and range
## ----keyresult, echo=F---------------------------------------------------
## Store key results
#Cumulative staff and student infections for each run

# Combine scenarios of interest into list of data frames
list_res<-lapply(list("base"=mod_base,"test4"=mod_4test,"screen30"= mod_30screen,"combo"=mod_combo),as.data.frame)

list_res<-lapply(names(c(list_res)),function(x) {
  list_res[[x]] %>%
    mutate(scenario = x)
})

# Process the output into plots for active and cumulative cases for stadd and student
# Default is median with IQR as upper and lower bounds but can be changed
p1<-psa_plot(list_res[[1]],title1="Active cases-base", title2="Total cases-base")
p2<- psa_plot(list_res[[2]],title1="Active cases-4 day test delay only", title2 = "Total-4 day test delay only")
p3 <-psa_plot(list_res[[3]], title1="Active cases-30day screen interval", title2="Total-30day screen interval")
p4 <- psa_plot(list_res[[4]], title1= "Active cases-Comb. test & screen", title2= "Total-Comb. test&screen")

# Arrange plots for output

png("Plots/5_psa_plots_redo1.png", units="in", width=6, height=5, res=500)
grid.arrange(p1[[1]]+rremove("legend"),p1[[2]]+rremove("legend"),p2[[1]]+rremove("legend"),p2[[2]]+rremove("legend"))
dev.off()

png("Plots/5_psa_plots_redo2.png", units="in", width=6, height=5, res=500)
grid.arrange(p3[[1]]+rremove("legend"),p3[[2]]+rremove("legend"),p4[[1]]+rremove("legend"),p4[[2]]+rremove("legend"))
dev.off()

```

```{r}
#' 
#' ## Totals -- at end of the semester
## ------------------------------------------------------------------------
df_cum<-do.call("rbind",list_res)%>%
      filter(time == max(time)) %>%
      group_by(scenario, run) %>%
      summarize(
        student_n = S_on + E_on + I_on + R_on + P_on + Q_on - Dcum_on +
                    S_off + E_off + I_off + R_off + P_off + Q_off - Dcum_off,
        student_cases = Icum_on + Icum_off,
        student_hosps = Hcum_on + Hcum_off,
        student_isos = Pcum_on + Pcum_off,
        student_quas = Qcum_on + Qcum_off,
        student_deaths = Dcum_on + Dcum_off,
        saf_n = S_saf + E_saf + I_saf + R_saf + P_saf + Q_saf - Dcum_saf,
        saf_cases = Icum_saf,
        saf_hosps = Hcum_saf,
        saf_deaths = Dcum_saf,
        tests = Test
      ) %>%
      ungroup() %>%
      mutate(tests_pc = tests / (student_n + saf_n))

    df_peak <- do.call("rbind",list_res) %>%
      group_by(scenario, run) %>%
      summarize(
        student_cases_peak = max(Isym_on + Isym_off, na.rm = TRUE),
        student_isos_peak = max(P_on + P_off, na.rm = TRUE),
        student_isos_days = sum(P_on + P_off, na.rm = TRUE),
        student_quas_peak = max(Q_on + Q_off, na.rm = TRUE),
        student_quas_days = sum(Q_on + Q_off, na.rm = TRUE),
        saf_cases_peak = max(Isym_saf, na.rm = TRUE),
        ) %>%
      ungroup()

    df_out <- full_join(df_cum, df_peak, by = c("scenario", "run")) %>%
      pivot_longer(
        -c(scenario, run),
        names_to = "measure",
        values_to = "value"
      ) %>%
      group_by(measure, scenario) %>%
      summarize(
        low = quantile(value, 0.025, na.rm = TRUE),
        med = quantile(value, 0.5, na.rm = TRUE),
        high = quantile(value, 0.975, na.rm = TRUE)
      ) %>%
      mutate(value = paste0(
        round(med,digits=0),
        " (", round(low,digits=0), "-",
        round(high, digits=0), ")",sep="")
        ) %>%
      pivot_wider(id_cols = measure, names_from = scenario, values_from = value)

 
var_order <- data.frame(measure=c("student_cases","student_cases_peak","student_hosps","student_deaths","student_isos","student_isos_peak","student_isos_days",
               "student_quas","student_quas_peak","student_quas_days","saf_cases","saf_cases_peak","saf_hosps","saf_deaths","tests","tests_pc")) %>%
            mutate(measure= as.character(measure))

df_out<-var_order %>% left_join(df_out, by="measure")


kable(df_out[,c("measure","base","test4","screen30","combo")], digits = 0, align = "c") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed")) 


#saveRDS(df_out[,c("measure","base","test4","screen30","combo")],"tables/res_scen.RDS")
write.csv(df_out[,c("measure","base","test4","screen30","combo")],"tables/res_scen_redo.csv")


```

```{r}
## Distribution of proportion of cumulative cases due to campus transmission 

propcomm<-do.call("rbind",list_res)%>%
      filter(time == max(time)) %>%
      #select(Icum_on_camp, Icum_off_camp,Icum_stu, Icum_saf_camp, Icum_saf) %>%
      mutate(Icum_stu_camp = Icum_on_camp + Icum_off_camp,
             prop_stu_camp = Icum_stu_camp/Icum_stu,
             prop_saf_camp = Icum_saf_camp/Icum_saf) %>%
      select(scenario, Icum_stu_camp, Icum_stu, prop_stu_camp,prop_saf_camp)%>%
      group_by(scenario) %>%
        summarize(
        low_stu =1-quantile(prop_stu_camp, 0.975, na.rm = TRUE),
        med_stu = 1-quantile(prop_stu_camp, 0.5, na.rm = TRUE),
        high_stu = 1-quantile(prop_stu_camp, 0.025, na.rm = TRUE),
        low_saf = 1-quantile(prop_saf_camp, 0.975, na.rm=TRUE),
        med_saf = 1-quantile(prop_saf_camp, 0.5, na.rm = TRUE),
        high_saf = 1-quantile(prop_saf_camp, 0.025, na.rm=TRUE)
      ) 

propcomm1 <- do.call("rbind",list_res)%>%
      filter(time == max(time)) %>%
      #select(Icum_on_camp, Icum_off_camp,Icum_stu, Icum_saf_camp, Icum_saf) %>%
      mutate(Icum_stu_camp = Icum_on_camp + Icum_off_camp,
             prop_stu_comm = 1-Icum_stu_camp/Icum_stu,
             prop_saf_comm = 1-Icum_saf_camp/Icum_saf) %>%
      select(scenario, prop_stu_comm, prop_saf_comm) 

png("Plots/11222020/supp_boxplot.png", units="in", width=6, height=5, res=800)
propcomm1%>% melt() %>% 
  mutate(scenario = factor(scenario, levels=c("base","screen30","test4","combo"))) %>%
  ggplot() +geom_boxplot(aes(x=scenario, y=value, fill = variable))+
  theme_bw() +ylab("Prop. due to community transmission")+
  scale_fill_manual(values=alpha(c("darkblue","darkred"),0.2),
                    name="",
                    labels=c("Stud","Staff"))+
  ggtitle("Prop. of cumulative cases due to community transmission \nby intervention")
dev.off()

saveRDS(propcomm,"tables/propcomm.RDS")

```

```{r}
#' ## PRCC
## ----prcc----------------------------------------------------------------
results<-data.frame(matrix(ncol=2,nrow=total.set.size))
results[,1]<-as.data.frame(mod_combo) %>% filter(time==116) %>% select(Icum_stu)
results[,2]<-as.data.frame(mod_combo) %>% filter(time==116) %>% select(Icum_saf)
```

```{r}
## Scatter plots for PRCC assumption check


df <- data.frame(matrix(nrow=1000,ncol=18))
for (i in 1:18){
  df[,i] <- mod_combo$param[i]
  colnames(df)[i] <- names(mod_combo$param[i])
}

df <- cbind(df,results)
df<- df %>% select(-isolation)
colnames(df)[18:19]<- c("cum_stud","cum_staff")
df$sim <- seq(1:1000)
df <- df%>% rename("beta_stu_to_stu"="beta_student_to_student")

png("Plots/5_input_output_scatterplots_stud.png", units="in", width=8, height=7, res=700)
df %>% select(-cum_staff) %>% melt(id.vars = c("sim","cum_stud")) %>%
  ggplot(aes(x=value, y= cum_stud)) + geom_point() + facet_wrap(~variable, scales = "free")+
 theme_bw()+theme(plot.title = element_text(size=7),
        axis.text.x = element_text(size=7, angle=45))+ylab("Cumulative student cases")
dev.off()

png("Plots/5_input_output_scatterplots_staff.png", units="in", width=8, height=7, res=700)
df %>% select(-cum_stud) %>% melt(id.vars = c("sim","cum_staff")) %>%
  ggplot(aes(x=value, y= cum_staff)) + geom_point() + facet_wrap(~variable, scales = "free") +theme_bw()+
 theme(plot.title = element_text(size=7),
  axis.text.x = element_text(size=7, angle=45))+ylab("Cumulative staff cases")
dev.off()



```


```{r}
# Actual PRCC test
variable <- data.frame(R0_StoS = R0_student_to_student.int,
                R0_OntoOn = R0_on_to_on.int,
                R0_saf = R0_saf.int,
                comm = community.int,
                eff_npi = eff_npi.int,
                asympt_stu = p_asympt_stu.int,
                asympt_saf = p_asympt_saf.int,
                contacts = contacts.int,
                p_contacts_reached = p_contacts_reached.int,
                ili = ili.int,
                sens = sensitivity.int,
                latent = latent.int,
                infectious = infectious.int)

variable_name <- c("R0 student to student","R0 on campus student to on campus student","R0 students to staff",
                   "Community introduction","Efficacy of NPI","Prop. asymptomatic (students)","Prop. asymptomatic (staff)",
                   "No. contacts per case","Proportion contacts reached","Proportion with ILI symptoms per month",
                   "PCR sensitivity","Latent period (days)", "Infectious period (days)")

pcc_res<-list()

for (i in 1:ncol(results)){
  pcc_res[[i]]<-pcc(variable,results[,i], rank= T, nboot =total.set.size)$PRCC
}

df_res<-as.data.frame(do.call(rbind,pcc_res)) %>% select(original)
df_res$output <-rep(rep(c("Student","Staff"),each =13),1)
df_res$output <-as.factor(df_res$output)
df_res$var<-rep(colnames(variable),2)

png("Plots/5_prcc_redo.png", units="in", width=7, height=5, res=600)
ggplot(data=df_res, aes(x= var,y = original,fill=output)) +
  geom_bar(position="dodge",stat="identity") +
  scale_fill_grey()+
  scale_x_discrete(limits=colnames(variable),labels=variable_name)+
  coord_flip() + theme_classic()+
  ggtitle("Partial rank correlation coefficient of variables")+xlab("Variable")+ylab("")+theme(legend.title=element_blank())
dev.off()

saveRDS(pcc_res,"tables/pcc_res.RDS")
```

```{r}
#' ## Gut checks for sampled distributions
#' 
## ------------------------------------------------------------------------
list.int<-list(R0_student_to_student=R0_student_to_student.int,R0_on_to_on = R0_on_to_on.int,R0_saf = R0_saf.int,
               community = community.int,p_asympt_stu = p_asympt_stu.int,p_asympt_saf = p_asympt_saf.int,contacts = contacts.int,
               ili = ili.int,sensitivity = sensitivity.int, p_contact_reached = p_contacts_reached.int,latent = latent.int,infectious = infectious.int,
               p_hosp_stu = p_hosp_stu.int,p_hosp_saf = p_hosp_saf.int,p_death_stu = p_death_stu.int,p_death_saf=p_death_saf.int,
               eff_npi = eff_npi.int,beta_student_to_student = beta_student_to_student.int,beta_on_to_on = beta_on_to_on.int,beta_saf = beta_saf.int)

par(mfrow = c(2,2))
for (i in 1:length(list.int)){
  p1<-plot(density(x=list.int[[i]]),main=names(list.int[i]))
  print(p1)
}  

summary_parm<-as.data.frame(matrix(data=0,nrow=length(list.int),ncol=3))
colnames(summary_parm) <- c("2.5%","50%","97.5%")

for (i in 1:length(list.int)){
  summary_parm[i,]<-quantile(list.int[[i]],p=c(0.025,0.5,0.975))
  row.names(summary_parm)[i] <- names(list.int[i])
}

summary_parm


#' 
#' Think for the most part, these are aligned with the parameter table apart from a few like the upper bound of sensitivity 
#' 
```

```{r}
#' 
#' ## One time screening
## ----one_time_screen---------------------------------------------------
## Scenario with one time screening across entire semester
param <- param.dcm(latent = latent.int,
                   infectious = infectious.int,
                   isolation = isolation,
                   beta_student_to_student = beta_student_to_student.int,
                   beta_on_to_on = beta_on_to_on.int,
                   beta_saf = beta_saf.int,
                   community = community.int,
                   p_asympt_stu = p_asympt_stu.int,
                   p_asympt_saf = p_asympt_saf.int,
                   p_hosp_stu = p_hosp_stu.int,
                   p_hosp_saf = p_hosp_saf.int,
                   p_death_stu = p_death_stu.int,
                   p_death_saf = p_death_saf.int,
                   contacts = contacts.int,
                   p_contacts_reached = p_contacts_reached.int,
                   ili = ili.int,
                   sensitivity = sensitivity.int,
                   eff_npi = eff_npi.int,
                   testing = 0,
                   screening = 1/116,
                   screening_on=1/116)

mod_116screen <- dcm(param, init, control)
mod_116screen <- mutate_epi(mod_116screen, I_stu = Isym_on + Isym_off,
                        Icum_stu = Icum_on + Icum_off,
                        P_stu = P_on + P_off,
                        Pcum_stu = Pcum_on + Pcum_off,
                        Qcum_stu = Qcum_on + Qcum_off,
                        Hcum_stu = Hcum_on + Hcum_off,
                        Dcum_stu = Dcum_on + Dcum_off)

#' ## Seven day screening interval
## ----seven_day_screen---------------------------------------------------
## Scenario with one time screening across entire semester
param <- param.dcm(latent = latent.int,
                   infectious = infectious.int,
                   isolation = isolation,
                   beta_student_to_student = beta_student_to_student.int,
                   beta_on_to_on = beta_on_to_on.int,
                   beta_saf = beta_saf.int,
                   community = community.int,
                   p_asympt_stu = p_asympt_stu.int,
                   p_asympt_saf = p_asympt_saf.int,
                   p_hosp_stu = p_hosp_stu.int,
                   p_hosp_saf = p_hosp_saf.int,
                   p_death_stu = p_death_stu.int,
                   p_death_saf = p_death_saf.int,
                   contacts = contacts.int,
                   p_contacts_reached = p_contacts_reached.int,
                   ili = ili.int,
                   sensitivity = sensitivity.int,
                   eff_npi = eff_npi.int,
                   testing = 0,
                   screening = 1/7,
                   screening_on=1/7)

mod_7screen <- dcm(param, init, control)
mod_7screen <- mutate_epi(mod_7screen, I_stu = Isym_on + Isym_off,
                            Icum_stu = Icum_on + Icum_off,
                            P_stu = P_on + P_off,
                            Pcum_stu = Pcum_on + Pcum_off,
                            Qcum_stu = Qcum_on + Qcum_off,
                            Hcum_stu = Hcum_on + Hcum_off,
                            Dcum_stu = Dcum_on + Dcum_off)
#' ## Two day testing
## ----two_day_test---------------------------------------------------
## Scenario with one time screening across entire semester
param <- param.dcm(latent = latent.int,
                   infectious = infectious.int,
                   isolation = isolation,
                   beta_student_to_student = beta_student_to_student.int,
                   beta_on_to_on = beta_on_to_on.int,
                   beta_saf = beta_saf.int,
                   community = community.int,
                   p_asympt_stu = p_asympt_stu.int,
                   p_asympt_saf = p_asympt_saf.int,
                   p_hosp_stu = p_hosp_stu.int,
                   p_hosp_saf = p_hosp_saf.int,
                   p_death_stu = p_death_stu.int,
                   p_death_saf = p_death_saf.int,
                   contacts = contacts.int,
                   p_contacts_reached = p_contacts_reached.int,
                   ili = ili.int,
                   sensitivity = sensitivity_2.int,
                   eff_npi = eff_npi.int,
                   testing = 1/2,
                   screening = 0,
                   screening_on=0)

mod_2test <- dcm(param, init, control)
mod_2test <- mutate_epi(mod_2test, I_stu = Isym_on + Isym_off,
                            Icum_stu = Icum_on + Icum_off,
                            P_stu = P_on + P_off,
                            Pcum_stu = Pcum_on + Pcum_off,
                            Qcum_stu = Qcum_on + Qcum_off,
                            Hcum_stu = Hcum_on + Hcum_off,
                            Dcum_stu = Dcum_on + Dcum_off)

#' ## Seven day test interval
## ----seven_day_test---------------------------------------------------
## Scenario with one time screening across entire semester
param <- param.dcm(latent = latent.int,
                   infectious = infectious.int,
                   isolation = isolation,
                   beta_student_to_student = beta_student_to_student.int,
                   beta_on_to_on = beta_on_to_on.int,
                   beta_saf = beta_saf.int,
                   community = community.int,
                   p_asympt_stu = p_asympt_stu.int,
                   p_asympt_saf = p_asympt_saf.int,
                   p_hosp_stu = p_hosp_stu.int,
                   p_hosp_saf = p_hosp_saf.int,
                   p_death_stu = p_death_stu.int,
                   p_death_saf = p_death_saf.int,
                   contacts = contacts.int,
                   p_contacts_reached = p_contacts_reached.int,
                   ili = ili.int,
                   sensitivity = sensitivity_7.int,
                   eff_npi = eff_npi.int,
                   testing = 1/7,
                   screening = 0,
                   screening_on=0)

mod_7test <- dcm(param, init, control)
mod_7test <- mutate_epi(mod_7test, I_stu = Isym_on + Isym_off,
                            Icum_stu = Icum_on + Icum_off,
                            P_stu = P_on + P_off,
                            Pcum_stu = Pcum_on + Pcum_off,
                            Qcum_stu = Qcum_on + Qcum_off,
                            Hcum_stu = Hcum_on + Hcum_off,
                            Dcum_stu = Dcum_on + Dcum_off)

# Base scenario with no control measures at all
param <- param.dcm(latent = latent.int,
                   infectious = infectious.int,
                   isolation = isolation,
                   beta_student_to_student = beta_student_to_student.int,
                   beta_on_to_on = beta_on_to_on.int,
                   beta_saf = beta_saf.int,
                   community = community.int,
                   p_asympt_stu = p_asympt_stu.int,
                   p_asympt_saf = p_asympt_saf.int,
                   p_hosp_stu = p_hosp_stu.int,
                   p_hosp_saf = p_hosp_saf.int,
                   p_death_stu = p_death_stu.int,
                   p_death_saf = p_death_saf.int,
                   contacts = contacts.int,
                   p_contacts_reached = p_contacts_reached.int,
                   ili = ili.int,
                   eff_npi = 0,
                   sensitivity = sensitivity.int,
                   testing = 0,
                   screening = 0,
                   screening_on=0)

mod_nocont <- dcm(param, init, control)
mod_nocont<- mutate_epi(mod_nocont, I_stu = Isym_on + Isym_off,
                        Icum_stu = Icum_on + Icum_off,
                        P_stu = P_on + P_off,
                        Pcum_stu = Pcum_on + Pcum_off,
                        Qcum_stu = Qcum_on + Qcum_off,
                        Hcum_stu = Hcum_on + Hcum_off,
                        Dcum_stu = Dcum_on + Dcum_off)
mod_nocont<-as.data.frame(mod_nocont)

df_cum1<-mod_nocont%>%
      filter(time == max(time)) %>%
      group_by(run) %>%
      summarize(
        student_n = S_on + E_on + I_on + R_on + P_on + Q_on - Dcum_on +
                    S_off + E_off + I_off + R_off + P_off + Q_off - Dcum_off,
        student_cases = Icum_on + Icum_off,
        student_hosps = Hcum_on + Hcum_off,
        student_isos = Pcum_on + Pcum_off,
        student_quas = Qcum_on + Qcum_off,
        student_deaths = Dcum_on + Dcum_off,
        saf_n = S_saf + E_saf + I_saf + R_saf + P_saf + Q_saf - Dcum_saf,
        saf_cases = Icum_saf,
        saf_hosps = Hcum_saf,
        saf_deaths = Dcum_saf,
        tests = Test
      ) %>%
      ungroup() %>%
      mutate(tests_pc = tests / (student_n + saf_n))

    df_peak1 <- mod_nocont %>%
      group_by(run) %>%
      summarize(
        student_cases_peak = max(Isym_on + Isym_off, na.rm = TRUE),
        student_isos_peak = max(P_on + P_off, na.rm = TRUE),
        student_isos_days = sum(P_on + P_off, na.rm = TRUE),
        student_quas_peak = max(Q_on + Q_off, na.rm = TRUE),
        student_quas_days = sum(Q_on + Q_off, na.rm = TRUE),
        saf_cases_peak = max(Isym_saf, na.rm = TRUE),
        ) %>%
      ungroup()

    df_nocont <- full_join(df_cum1, df_peak1, by = c("run")) %>%
      pivot_longer(
        -c(run),
        names_to = "measure",
        values_to = "value"
      ) %>%
      group_by(measure) %>%
      summarize(
        low = quantile(value, 0.025, na.rm = TRUE),
        med = quantile(value, 0.5, na.rm = TRUE),
        high = quantile(value, 0.975, na.rm = TRUE)
      ) %>%
      mutate(value = paste0(
        round(med,digits=0),
        " (", round(low,digits=0), "-",
        round(high, digits=0), ")",sep=""))

#saveRDS(df_nocont,"tables/res_nocont.RDS")


```

```{r}
# Combine all scenarios for cumulative incidence, % reduction and ranges for staff and student
list_res1<-lapply(list(
                      "base"=mod_base,
                      "test2"=mod_2test,
                      "test4"=mod_4test,
                      "test7"=mod_7test,
                      "screen116"=mod_116screen,
                      "screen30"= mod_30screen,
                      "screen7"=mod_7screen,
                      "combo"=mod_combo),as.data.frame)

list_res1<-lapply(names(c(list_res1)),function(x) {
  list_res1[[x]] %>%
    mutate(scenario = x)
})

df_cum<-do.call("rbind",list_res1)%>%
  filter(time == max(time)) %>%
  group_by(scenario, run) %>%
  summarize(
    student_cases = Icum_on + Icum_off,
    saf_cases = Icum_saf,
  ) %>%
  pivot_longer(
    -c(scenario, run),
    names_to = "measure",
    values_to = "value"
  ) %>%
  group_by(measure, scenario) %>%
  summarize(
    low = quantile(value, 0.025, na.rm = TRUE),
    med = quantile(value, 0.5, na.rm = TRUE),
    high = quantile(value, 0.975, na.rm = TRUE)
  ) %>%
  mutate(value = paste0(
    round(med,digits=0),
    " (", round(low,digits=0), "-",
    round(high, digits=0), ")",sep="")
  ) 

df_stu <- df_cum %>% filter(measure=="student_cases") %>%
  mutate(
    perc_red = round((1-med/3068)*100,digits=0)
  )

df_saf <- df_cum %>% filter(measure=="saf_cases") %>%
  mutate(
    perc_red = round((1-med/1063)*100,digits=0)
  )

df_summary <- full_join(df_stu,df_saf)  %>%
    pivot_wider(id_cols = scenario, names_from = measure, values_from = c(value,perc_red))

df_summary

write.csv(df_summary,"tables/res_summary_redo.csv")

#saveRDS(df_summary,"tables/res_summary.RDS")
```