###################################################################### ### Measles Outbreak Model - Cases Averted by Mass Vaccination ### ### Code by Paul Gastanaduy, Prabasaj Paul, Nina Masters ### ### Updated September 2, 2022 ### ###################################################################### # 1 Load packages and set up environment ------------------------------------------------------------------------------------------------------------------- library(tidyverse) library(readxl) library(styler) library(adaptivetau) library(Hmisc) library(ggthemes) library(hrbrthemes) library(gridExtra) # * 1.1 Source Ro, infectious, pre-infectious period parameter definitions & SEIR and SEIRv transitions and rate functions source("Input Parameters/SEIR_and_SEIRv_Transitions_and_Rate_Functions.R") # * 1.2 Specify population size in each age group source("Input Parameters/Ft_McCoy_Initial_Conditions_and_Age_Stratified_Population.R") # Ft. McCoy (22 cases) # * 1.3 Set other model specs seed = 111 #set random seed for reproducibility Nruns = 1000 #define number of model runs # * 1.4 create initial condition named vectors varvec = c(S1 = S1, E1 = E1, I1 = I1, R1 = R1, Cases1 = Cases1, S2 = S2, E2 = E2, I2 = I2, R2 = R2, Cases2 = Cases2, S3 = S3, E3 = E3, I3 = I3, R3 = R3, Cases3 = Cases3, S4 = S4, E4 = E4, I4 = I4, R4 = R4, Cases4 = Cases4, S5 = S5, E5 = E5, I5 = I5, R5 = R5, Cases5 = Cases5) parvec = c(b11 = b11, b12 = b12, b13 = b13, b14 = b14, b15 = b15, b21 = b21, b22 = b22, b23 = b23, b24 = b24, b25 = b25, b31 = b31, b32 = b32, b33 = b33, b34 = b34, b35 = b35, b41 = b41, b42 = b42, b43 = b43, b44 = b44, b45 = b45, b51 = b51, b52 = b52, b53 = b53, b54 = b54, b55 = b55, gamma, theta) # 2 Load in importation and vaccination files --------------------------------------------------------------------------------\ day.zero=as.Date("08/24/2021",format="%m/%d/%Y") # Ft McCoy # vax.doses=read_csv("Input Parameters/Mccoy/vax_mccoy.csv") %>% mutate(date_vax = as.Date(date_vax, format="%m/%d/%Y")) %>% mutate(OBdate = date_vax - day.zero) %>% arrange((OBdate)) %>% filter(OBdate <201) # vax.doses.7daydelay=read_csv("Input Parameters/Mccoy/vax_mccoy_7_day_delay.csv") %>% mutate(date_vax = as.Date(date_vax, format="%m/%d/%Y")) %>% mutate(OBdate = date_vax - day.zero)%>% arrange((OBdate)) %>% filter(OBdate <201) # vax.doses.14daydelay=read_csv("Input Parameters/Mccoy/vax_mccoy_14_day_delay.csv") %>% mutate(date_vax = as.Date(date_vax, format="%m/%d/%Y")) %>% mutate(OBdate = date_vax - day.zero)%>% arrange((OBdate)) %>% filter(OBdate <201) # vax.doses.21daydelay=read_csv("Input Parameters/Mccoy/vax_mccoy_21_day_delay.csv") %>% mutate(date_vax = as.Date(date_vax, format="%m/%d/%Y")) %>% mutate(OBdate = date_vax - day.zero)%>% arrange((OBdate)) %>% filter(OBdate <201) # 2.1 Run model with no vaccination ------------------------------------------------------------------------------------------------\ allruns.novax=c() set.seed(seed) for (myrun in 1:Nruns) { run=ssa.adaptivetau(init.values = varvec, transitions = SEIR_transitions, rateFunc = SEIR_Model, params = parvec, tf = 200, tl.params = list(maxtau=1)) #change to 0, will only be a one row data frame, get the data in the right format allruns.novax=rbind(allruns.novax,cbind(Run=myrun,run)) for (i in 2:nrow(vax.doses)) { last.run=tail(run,1) #define last run last.run[4]=last.run[4]+vax.doses$import_1[i] ##add any importations into I1 last.run[9]=last.run[9]+vax.doses$import_2[i] ##add any importations into I2 last.run[14]=last.run[14]+vax.doses$import_3[i] ##add any importations into I3 last.run[19]=last.run[19]+vax.doses$import_4[i] ##add any importations into I4 last.run[6]=last.run[6]+vax.doses$import_1[i] ##add any importations into Cases1 last.run[11]=last.run[11]+vax.doses$import_2[i] ##add any importations into Cases2 last.run[16]=last.run[16]+vax.doses$import_3[i] ##add any importations into Cases3 last.run[21]=last.run[21]+vax.doses$import_4[i] ##add any importations into Cases4 sub.run=ssa.adaptivetau(init.values = c(S1 = last.run[2], E1 = last.run[3], I1 = last.run[4], R1 = last.run[5], Cases1 = last.run[6], S2 = last.run[7], E2 = last.run[8], I2 = last.run[9], R2 = last.run[10], Cases2 = last.run[11], S3 = last.run[12], E3 = last.run[13], I3 = last.run[14], R3 = last.run[15], Cases3 = last.run[16], S4 = last.run[17], E4 = last.run[18], I4 = last.run[19], R4 = last.run[20], Cases4 = last.run[21], S5 = last.run[22], E5 = last.run[23], I5 = last.run[24], R5 = last.run[25], Cases5 = last.run[26]), transitions = SEIR_transitions,rateFunc = SEIR_Model, params = parvec, tf = vax.doses$OBdate[i] - vax.doses$OBdate[i-1] ,tl.params = list(maxtau=1)) sub.run[,1]=sub.run[,1]+vax.doses$OBdate[i-1] run=rbind(run,sub.run)} allruns.novax=rbind(allruns.novax,cbind(Run=myrun,run)) } allruns.novax <- data.frame(allruns.novax) %>% mutate(I = I1 + I2 + I3 + I4 + I5) %>% mutate(Cases = Cases1 + Cases2 + Cases3 + Cases4 + Cases5) %>% mutate(N = S1 + S2 + S3 + S4 + S5 + E1 + E2 + E3 + E4 + E5 + I1 + I2 + I3 + I4 + I5 + R1 + R2 + R3 + R4 + R5) # * 2.2 post-processing of simulation runs -------------------------------------------------------------------- temp.novax =with(allruns.novax,tapply(Cases,Run,max)) temp1.novax=with(allruns.novax,tapply(Cases1,Run,max)) temp2.novax=with(allruns.novax,tapply(Cases2,Run,max)) temp3.novax=with(allruns.novax,tapply(Cases3,Run,max)) temp4.novax=with(allruns.novax,tapply(Cases4,Run,max)) temp5.novax=with(allruns.novax,tapply(Cases5,Run,max)) #can determine time of outbreak - once E and I are both 0 (would give outbreak duration in days) results.novax <- data.frame(temp.novax,temp1.novax,temp2.novax,temp3.novax,temp4.novax,temp5.novax, row.names(temp.novax)) %>% rename(TotalCases = temp.novax) %>% rename(Cases1 = temp1.novax) %>% rename(Cases2 = temp2.novax) %>% rename(Cases3 = temp3.novax) %>% rename(Cases4 = temp4.novax) %>% rename(Cases5 = temp5.novax) %>% rename(Run = row.names.temp.novax.) median.novax <- median(results.novax$TotalCases) AR_novax = median.novax/(S1 + S2 + S3 + S4 + S5) # write_csv(results.novax, file = "Output Files/mccoy_no_vaccination.csv") # 3: Vaccination Intervention ------------------------------------------------------------------------------------------------------- VE = 0.925 VE_under12m = 0.84 # create initial condition named vectors -------------------------------------------------------------------------------- varvecv = c(S1 = S1, E1 = E1, I1 = I1, R1 = R1, Cases1 = Cases1, S2 = S2, Sv2 = Sv2, E2 = E2, I2 = I2, R2 = R2, Cases2 = Cases2, S3 = S3, Sv3 = Sv3, E3 = E3, I3 = I3, R3 = R3, Cases3 = Cases3, S4 = S4, Sv4 = Sv4, E4 = E4, I4 = I4, R4 = R4, Cases4 = Cases4, S5 = S5, E5 = E5, I5 = I5, R5 = R5, Cases5 = Cases5) parvecv = c(b11 = b11, b12 = b12, b13 = b13, b14 = b14, b15 = b15, b21 = b21, b22 = b22, b23 = b23, b24 = b24, b25 = b25, b31 = b31, b32 = b32, b33 = b33, b34 = b34, b35 = b35, b41 = b41, b42 = b42, b43 = b43, b44 = b44, b45 = b45, b51 = b51, b52 = b52, b53 = b53, b54 = b54, b55 = b55, gamma, theta) # 4: Vaccination with 7 day delay: acting as base case with 7-day delay for immunity to kick in -------------------------------------------------------------------------------------------------------------------------------- allruns.vax.7daydelay=c() set.seed(seed) for (myrun in 1:Nruns) { run=ssa.adaptivetau(init.values = varvecv, transitions = SEIRv_transitions, rateFunc = SEIRv_Model, params = parvecv, tf = 0, tl.params = list(maxtau=1)) allruns.vax.7daydelay=rbind(allruns.vax.7daydelay,cbind(Run=myrun,run)) for (i in 2:nrow(vax.doses.7daydelay)) { last.run=tail(run,1) #define last run #note: vaccination files define cumulative doses on day i as cumulate doses given from the day before (i.e. day i -1) probability_vaxxed = ifelse(vax.doses.7daydelay$cumulative_doses[i] < N_eligible_vax, min(vax.doses.7daydelay$doses[i]/(N_eligible_vax - vax.doses.7daydelay$cumulative_doses[i]),1),0) last.run[4] = last.run[4] +vax.doses.7daydelay$import_1[i] #add any importations into I1 last.run[10] = last.run[10] +vax.doses.7daydelay$import_2[i] #add any importations into I2 last.run[16] = last.run[16] +vax.doses.7daydelay$import_3[i] #add any importations into I3 last.run[22] = last.run[22] +vax.doses.7daydelay$import_4[i] #add any importations into I4 last.run[6] = last.run[6] + vax.doses.7daydelay$import_1[i] #add any importations into Cases1 last.run[12] = last.run[12] +vax.doses.7daydelay$import_2[i] #add any importations into Cases2 last.run[18] = last.run[18] +vax.doses.7daydelay$import_3[i] #add any importations into Cases3 last.run[24] = last.run[24] +vax.doses.7daydelay$import_4[i] #add any importations into Cases4 last.run[8] = last.run[8] + floor(probability_vaxxed*last.run[7]*(1-VE_under12m)) #Sv2 last.run[14] = last.run[14] + floor(probability_vaxxed*last.run[13]*(1-VE)) #Sv3 last.run[20] = last.run[20] + floor(probability_vaxxed*last.run[19]*(1-VE)) #Sv4 last.run[11] = last.run[11] + floor(probability_vaxxed*last.run[7]*VE_under12m) #R2 last.run[17] = last.run[17] + floor(probability_vaxxed*last.run[13]*VE) #R3 last.run[23] = last.run[23] + floor(probability_vaxxed*last.run[19]*VE) #R4 last.run[2] = last.run[2] #S1 last.run[7] = last.run[7] - floor(probability_vaxxed*last.run[7]*VE_under12m) - floor(probability_vaxxed*last.run[7]*(1-VE_under12m)) #S2 last.run[13] = last.run[13] - floor(probability_vaxxed*last.run[13]*VE) - floor(probability_vaxxed*last.run[13]*(1-VE)) #S3 last.run[19] = last.run[19] - floor(probability_vaxxed*last.run[19]*VE) - floor(probability_vaxxed*last.run[19]*(1-VE)) #S4 sub.run=ssa.adaptivetau(init.values = c(S1 = last.run[2], E1 = last.run[3], I1 = last.run[4], R1 = last.run[5], Cases1 = last.run[6], S2 = last.run[7], Sv2 = last.run[8], E2 = last.run[9], I2 = last.run[10], R2 = last.run[11], Cases2 = last.run[12], S3 = last.run[13], Sv3 = last.run[14], E3 = last.run[15], I3 = last.run[16], R3 = last.run[17], Cases3 = last.run[18], S4 = last.run[19], Sv4 = last.run[20], E4 = last.run[21], I4 = last.run[22], R4 = last.run[23], Cases4 = last.run[24], S5 = last.run[25], E5 = last.run[26], I5 = last.run[27], R5 = last.run[28], Cases5 = last.run[29]), transitions = SEIRv_transitions,rateFunc = SEIRv_Model, params = parvecv,tf = vax.doses.7daydelay$OBdate[i] - vax.doses.7daydelay$OBdate[i-1] ,tl.params = list(maxtau=1)) sub.run[,1]=sub.run[,1]+vax.doses.7daydelay$OBdate[i-1] run=rbind(run,sub.run)} allruns.vax.7daydelay=rbind(allruns.vax.7daydelay,cbind(Run=myrun,run)) } # Post-processing of runs ------------------------------------------- ----------------------------------------------------------------------- allruns.vax.7daydelay <- data.frame(allruns.vax.7daydelay) %>% mutate(I = I1 + I2 + I3 + I4 + I5) %>% mutate(Cases = Cases1 + Cases2 + Cases3 + Cases4 + Cases5) %>% mutate(N = S1 + S2 + S3 + S4 + S5 + E1 + E2 + E3 + E4 + E5 + I1 + I2 + I3 + I4 + I5 + R1 + R2 + R3 + R4 + R5) temp.vax.7daydelay=with(allruns.vax.7daydelay,tapply(Cases,Run,max)) temp1.vax.7daydelay=with(allruns.vax.7daydelay,tapply(Cases1,Run,max)) temp2.vax.7daydelay=with(allruns.vax.7daydelay,tapply(Cases2,Run,max)) temp3.vax.7daydelay=with(allruns.vax.7daydelay,tapply(Cases3,Run,max)) temp4.vax.7daydelay=with(allruns.vax.7daydelay,tapply(Cases4,Run,max)) temp5.vax.7daydelay=with(allruns.vax.7daydelay,tapply(Cases5,Run,max)) results.vax.all.7daydelay <- data.frame(temp.vax.7daydelay,temp1.vax.7daydelay,temp2.vax.7daydelay,temp3.vax.7daydelay,temp4.vax.7daydelay,temp5.vax.7daydelay, row.names(temp.vax.7daydelay)) %>% rename(TotalCases = temp.vax.7daydelay) %>% rename(Cases1 = temp1.vax.7daydelay) %>% rename(Cases2 = temp2.vax.7daydelay) %>% rename(Cases3 = temp3.vax.7daydelay) %>% rename(Cases4 = temp4.vax.7daydelay) %>% rename(Cases5 = temp5.vax.7daydelay) %>% rename(Run = row.names.temp.vax.7daydelay.) median.vax.all.7daydelay <- median(results.vax.all.7daydelay$TotalCases) median.vax.all.7daydelay/(S1+S2+S3+S4+S5) # write_csv(results.vax.all.7daydelay, file = "Output Files/mccoy_vaccination_no_delay.csv") # 5: Vaccination with 14 day delay -------------------------------------------------------------------------------------------------------------------------------- allruns.vax.14daydelay=c() set.seed(seed) for (myrun in 1:Nruns) { run=ssa.adaptivetau(init.values = varvecv, transitions = SEIRv_transitions, rateFunc = SEIRv_Model, params = parvecv, tf = 0, tl.params = list(maxtau=1)) allruns.vax.14daydelay=rbind(allruns.vax.14daydelay,cbind(Run=myrun,run)) for (i in 2:nrow(vax.doses.14daydelay)) { last.run=tail(run,1) #define last run probability_vaxxed = ifelse(vax.doses.14daydelay$cumulative_doses[i] < N_eligible_vax, min(vax.doses.14daydelay$doses[i]/(N_eligible_vax - vax.doses.14daydelay$cumulative_doses[i]),1),0) last.run[4] = last.run[4] +vax.doses.14daydelay$import_1[i] #add any importations into I1 last.run[10] = last.run[10] +vax.doses.14daydelay$import_2[i] #add any importations into I2 last.run[16] = last.run[16] +vax.doses.14daydelay$import_3[i] #add any importations into I3 last.run[22] = last.run[22] +vax.doses.14daydelay$import_4[i] #add any importations into I4 last.run[6] = last.run[6] +vax.doses.14daydelay$import_1[i] #add any importations into Cases1 last.run[12] = last.run[12] +vax.doses.14daydelay$import_2[i] #add any importations into Cases2 last.run[18] = last.run[18] +vax.doses.14daydelay$import_3[i] #add any importations into Cases3 last.run[24] = last.run[24] +vax.doses.14daydelay$import_4[i] #add any importations into Cases4 last.run[8] = last.run[8] + floor(probability_vaxxed*last.run[7]*(1-VE_under12m)) #Sv2 last.run[14] = last.run[14] + floor(probability_vaxxed*last.run[13]*(1-VE)) #Sv3 last.run[20] = last.run[20] + floor(probability_vaxxed*last.run[19]*(1-VE)) #Sv4 last.run[11] = last.run[11] + floor(probability_vaxxed*last.run[7]*VE_under12m) #R2 last.run[17] = last.run[17] + floor(probability_vaxxed*last.run[13]*VE) #R3 last.run[23] = last.run[23] + floor(probability_vaxxed*last.run[19]*VE) #R4 last.run[2] = last.run[2] #S1, remove any importations last.run[7] = last.run[7] - floor(probability_vaxxed*last.run[7]*VE_under12m) - floor(probability_vaxxed*last.run[7]*(1-VE_under12m)) #S2 last.run[13] = last.run[13] - floor(probability_vaxxed*last.run[13]*VE) - floor(probability_vaxxed*last.run[13]*(1-VE)) #S3 last.run[19] = last.run[19] - floor(probability_vaxxed*last.run[19]*VE) - floor(probability_vaxxed*last.run[19]*(1-VE)) #S4 sub.run=ssa.adaptivetau(init.values = c(S1 = last.run[2], E1 = last.run[3], I1 = last.run[4], R1 = last.run[5], Cases1 = last.run[6], S2 = last.run[7], Sv2 = last.run[8], E2 = last.run[9], I2 = last.run[10], R2 = last.run[11], Cases2 = last.run[12], S3 = last.run[13], Sv3 = last.run[14], E3 = last.run[15], I3 = last.run[16], R3 = last.run[17], Cases3 = last.run[18], S4 = last.run[19], Sv4 = last.run[20], E4 = last.run[21], I4 = last.run[22], R4 = last.run[23], Cases4 = last.run[24], S5 = last.run[25], E5 = last.run[26], I5 = last.run[27], R5 = last.run[28], Cases5 = last.run[29]), transitions = SEIRv_transitions,rateFunc = SEIRv_Model, params = parvecv,tf = vax.doses.14daydelay$OBdate[i] - vax.doses.14daydelay$OBdate[i-1] ,tl.params = list(maxtau=1)) sub.run[,1]=sub.run[,1]+vax.doses.14daydelay$OBdate[i-1] run=rbind(run,sub.run)} allruns.vax.14daydelay=rbind(allruns.vax.14daydelay,cbind(Run=myrun,run)) } # Post-processing of runs ------------------------------------------- ----------------------------------------------------------------------- allruns.vax.14daydelay <- data.frame(allruns.vax.14daydelay) %>% mutate(I = I1 + I2 + I3 + I4 + I5) %>% mutate(Cases = Cases1 + Cases2 + Cases3 + Cases4 + Cases5) %>% mutate(N = S1 + S2 + S3 + S4 + S5 + E1 + E2 + E3 + E4 + E5 + I1 + I2 + I3 + I4 + I5 + R1 + R2 + R3 + R4 + R5) temp.vax.14daydelay=with(allruns.vax.14daydelay,tapply(Cases,Run,max)) temp1.vax.14daydelay=with(allruns.vax.14daydelay,tapply(Cases1,Run,max)) temp2.vax.14daydelay=with(allruns.vax.14daydelay,tapply(Cases2,Run,max)) temp3.vax.14daydelay=with(allruns.vax.14daydelay,tapply(Cases3,Run,max)) temp4.vax.14daydelay=with(allruns.vax.14daydelay,tapply(Cases4,Run,max)) temp5.vax.14daydelay=with(allruns.vax.14daydelay,tapply(Cases5,Run,max)) results.vax.all.14daydelay <- data.frame(temp.vax.14daydelay,temp1.vax.14daydelay,temp2.vax.14daydelay,temp3.vax.14daydelay,temp4.vax.14daydelay,temp5.vax.14daydelay, row.names(temp.vax.14daydelay)) %>% rename(TotalCases = temp.vax.14daydelay) %>% rename(Cases1 = temp1.vax.14daydelay) %>% rename(Cases2 = temp2.vax.14daydelay) %>% rename(Cases3 = temp3.vax.14daydelay) %>% rename(Cases4 = temp4.vax.14daydelay) %>% rename(Cases5 = temp5.vax.14daydelay) %>% rename(Run = row.names.temp.vax.14daydelay.) median.vax.all.14daydelay <- median(results.vax.all.14daydelay$TotalCases) median.vax.all.14daydelay/(S1+S2+S3+S4+S5) # write_csv(results.vax.all.14daydelay, file = "Output Files/mccoy_vaccination_7_day_delay.csv") # 6: Vaccination with 21 day delay -------------------------------------------------------------------------------------------------------------------------------- allruns.vax.21daydelay=c() set.seed(seed) for (myrun in 1:Nruns) { run=ssa.adaptivetau(init.values = varvecv, transitions = SEIRv_transitions, rateFunc = SEIRv_Model, params = parvecv, tf = 0, tl.params = list(maxtau=1)) allruns.vax.21daydelay=rbind(allruns.vax.21daydelay,cbind(Run=myrun,run)) for (i in 2:nrow(vax.doses.21daydelay)) { last.run=tail(run,1) #define last run probability_vaxxed = ifelse(vax.doses.21daydelay$cumulative_doses[i] < N_eligible_vax, min(vax.doses.21daydelay$doses[i]/(N_eligible_vax - vax.doses.21daydelay$cumulative_doses[i]),1),0) last.run[4] = last.run[4] +vax.doses.21daydelay$import_1[i] #add any importations into I1 last.run[10] = last.run[10] +vax.doses.21daydelay$import_2[i] #add any importations into I2 last.run[16] = last.run[16] +vax.doses.21daydelay$import_3[i] #add any importations into I3 last.run[22] = last.run[22] +vax.doses.21daydelay$import_4[i] #add any importations into I4 last.run[6] = last.run[6] +vax.doses.21daydelay$import_1[i] #add any importations into Cases1 last.run[12] = last.run[12] +vax.doses.21daydelay$import_2[i] #add any importations into Cases2 last.run[18] = last.run[18] +vax.doses.21daydelay$import_3[i] #add any importations into Cases3 last.run[24] = last.run[24] +vax.doses.21daydelay$import_4[i] #add any importations into Cases4 last.run[8] = last.run[8] + floor(probability_vaxxed*last.run[7]*(1-VE_under12m)) #Sv2 last.run[14] = last.run[14] + floor(probability_vaxxed*last.run[13]*(1-VE)) #Sv3 last.run[20] = last.run[20] + floor(probability_vaxxed*last.run[19]*(1-VE)) #Sv4 last.run[11] = last.run[11] + floor(probability_vaxxed*last.run[7]*VE_under12m) #R2 last.run[17] = last.run[17] + floor(probability_vaxxed*last.run[13]*VE) #R3 last.run[23] = last.run[23] + floor(probability_vaxxed*last.run[19]*VE) #R4 last.run[2] = last.run[2] #S1, remove any importations last.run[7] = last.run[7] - floor(probability_vaxxed*last.run[7]*VE_under12m) - floor(probability_vaxxed*last.run[7]*(1-VE_under12m)) #S2 last.run[13] = last.run[13] - floor(probability_vaxxed*last.run[13]*VE) - floor(probability_vaxxed*last.run[13]*(1-VE)) #S3 last.run[19] = last.run[19] - floor(probability_vaxxed*last.run[19]*VE) - floor(probability_vaxxed*last.run[19]*(1-VE)) #S4 sub.run=ssa.adaptivetau(init.values = c(S1 = last.run[2], E1 = last.run[3], I1 = last.run[4], R1 = last.run[5], Cases1 = last.run[6], S2 = last.run[7], Sv2 = last.run[8], E2 = last.run[9], I2 = last.run[10], R2 = last.run[11], Cases2 = last.run[12], S3 = last.run[13], Sv3 = last.run[14], E3 = last.run[15], I3 = last.run[16], R3 = last.run[17], Cases3 = last.run[18], S4 = last.run[19], Sv4 = last.run[20], E4 = last.run[21], I4 = last.run[22], R4 = last.run[23], Cases4 = last.run[24], S5 = last.run[25], E5 = last.run[26], I5 = last.run[27], R5 = last.run[28], Cases5 = last.run[29]), transitions = SEIRv_transitions,rateFunc = SEIRv_Model, params = parvecv,tf = vax.doses.21daydelay$OBdate[i] - vax.doses.21daydelay$OBdate[i-1] ,tl.params = list(maxtau=1)) sub.run[,1]=sub.run[,1]+vax.doses.21daydelay$OBdate[i-1] run=rbind(run,sub.run)} allruns.vax.21daydelay=rbind(allruns.vax.21daydelay,cbind(Run=myrun,run)) } # * 6.2 Post-processing of runs ------------------------------------------- ----------------------------------------------------------------------- allruns.vax.21daydelay <- data.frame(allruns.vax.21daydelay) %>% mutate(I = I1 + I2 + I3 + I4 + I5) %>% mutate(Cases = Cases1 + Cases2 + Cases3 + Cases4 + Cases5) %>% mutate(N = S1 + S2 + S3 + S4 + S5 + E1 + E2 + E3 + E4 + E5 + I1 + I2 + I3 + I4 + I5 + R1 + R2 + R3 + R4 + R5) temp.vax.21daydelay=with(allruns.vax.21daydelay,tapply(Cases,Run,max)) temp1.vax.21daydelay=with(allruns.vax.21daydelay,tapply(Cases1,Run,max)) temp2.vax.21daydelay=with(allruns.vax.21daydelay,tapply(Cases2,Run,max)) temp3.vax.21daydelay=with(allruns.vax.21daydelay,tapply(Cases3,Run,max)) temp4.vax.21daydelay=with(allruns.vax.21daydelay,tapply(Cases4,Run,max)) temp5.vax.21daydelay=with(allruns.vax.21daydelay,tapply(Cases5,Run,max)) results.vax.all.21daydelay <- data.frame(temp.vax.21daydelay,temp1.vax.21daydelay,temp2.vax.21daydelay,temp3.vax.21daydelay,temp4.vax.21daydelay,temp5.vax.21daydelay, row.names(temp.vax.21daydelay)) %>% rename(TotalCases = temp.vax.21daydelay) %>% rename(Cases1 = temp1.vax.21daydelay) %>% rename(Cases2 = temp2.vax.21daydelay) %>% rename(Cases3 = temp3.vax.21daydelay) %>% rename(Cases4 = temp4.vax.21daydelay) %>% rename(Cases5 = temp5.vax.21daydelay) %>% rename(Run = row.names.temp.vax.21daydelay.) median.vax.all.21daydelay <- median(results.vax.all.21daydelay$TotalCases) median.vax.all.21daydelay/(S1+S2+S3+S4+S5) # write_csv(results.vax.all.21daydelay, file = "Output Files/mccoy_vaccination_14_day_delay.csv")