<jupyter_start><jupyter_text>4X4 Contact matrix<jupyter_code>#### adjusted matrix
read_excel("../data/modified_matrix.xlsx", sheet="modified") -> adj_contact
adj_contact %<>% rename('60s+' = '60+')
(as.matrix(adj_contact[,-1]) -> adj_contact)<jupyter_output><empty_output><jupyter_text>Incidence data<jupyter_code>options(warn=-1) 
read_excel("../data/db210913.xlsx") -> master_df
master_df %<>% rename(reported=신고일, diagnosis=진단일, onset=발병일, published=보도일,
                      variant=변이, age=연령, route=감염경로)

as.Date(master_df$reported) -> master_df$reported
as.Date(master_df$diagnosis) -> master_df$diagnosis
as.Date(master_df$onset) -> master_df$onset
as.Date(master_df$published) -> master_df$published


master_df %<>%
mutate(age_group = case_when(age < 20 ~ "1",
                             age >= 20 & age < 40 ~ "2",
                             age >= 40 & age < 60 ~ "3",
                             age >= 60 ~ "4"))

options(warn=0)

(Final_date <- as.Date(max(master_df$published)+1))<jupyter_output><empty_output><jupyter_text>Estimating the time delay from the illness onset to report<jupyter_code>master_df %>% filter(!is.na(onset)) %>% 
mutate(dist= as.numeric(published-onset)) %>% filter(dist>=0) -> master_delay

master_delay %>% dplyr::select(onset, dist) %>% filter(dist <= 30)-> dt.delay1
master_delay %>% mutate(distt = as.numeric(Final_date-onset)) %>% 
dplyr::select(onset, distt) %>% na.omit() -> dt.delay2

dt.delay1$dist <- as.numeric(dt.delay1$dist)
dt.delay2$distt <- as.numeric(dt.delay2$distt)
dt.delay1$dist[dt.delay1$dist<=0] <- 0.5

delay.llk <- function(param){
       
a0 = param[1]; b0 = param[2]; 
total.lik.gamma_n <- rep(0,length(dt.delay1))
    
    return(-sum(log((pweibull((dt.delay1$dist+0.5), shape = a0, scale = b0)-
                     pweibull(dt.delay1$dist-0.5, shape = a0, scale = b0))/
                   (pweibull(dt.delay2$distt, shape = a0, scale = b0)))))}

param00 =  c(2, 10)
options(warn=-1)
opt_delay <- optim(param00, fn=delay.llk, method='BFGS', control = list(maxit=100000, trace=0))
param <- opt_delay$par
options(warn=0)<jupyter_output><empty_output><jupyter_text>Daily number of confirmed cases in the Republic of Korea<jupyter_code>as.numeric(Final_date-1-min(master_df$published)) -> time.diff
ttime <- as.data.frame(c(0:time.diff))
colnames(ttime) <- c('t')
ttime %<>% mutate(published = min(master_df$published)+t)
master_df %>% filter(age_group==c("1")) %>% group_by(published) %>% count() %>% as.data.frame() -> data_1
merge(ttime, data_1, by='published', all.x=TRUE) -> data_1
data_1[is.na(data_1)] <- 0

master_df %>% filter(age_group==c("2")) %>% group_by(published) %>% count() %>% as.data.frame() -> data_2
merge(ttime, data_2, by='published', all.x=TRUE) -> data_2
data_2[is.na(data_2)] <- 0

master_df %>% filter(age_group==c("3")) %>% group_by(published) %>% count() %>% as.data.frame() -> data_3
merge(ttime, data_3, by='published', all.x=TRUE) -> data_3
data_3[is.na(data_3)] <- 0

master_df %>% filter(age_group==c("4")) %>% group_by(published) %>% count() %>% as.data.frame() -> data_4
merge(ttime, data_4, by='published', all.x=TRUE) -> data_4
data_4[is.na(data_4)] <- 0<jupyter_output><empty_output><jupyter_text>Variant data<jupyter_code>read.csv("../data/variant_05Jul2021.csv") -> master_prop
master_prop$확진일 <- ymd(master_prop$확진일)
master_prop %<>% filter(확진일 >= as.Date("2020-12-01") & 확진일 < as.Date("2021-06-21"))

read.csv("../data/variant_22Jul2021.csv") -> master_prop2
master_prop2$확진일 <- ymd(master_prop2$확진일)
master_prop2 %<>% filter(확진일 >= as.Date("2021-06-21") & 확진일 < as.Date("2021-07-11"))

rbind(master_prop, master_prop2) -> master_prop

master_prop$확진일 <- strftime(master_prop$확진일, "%V")

master_prop %<>% rename(variant=변이, date=확진일) %>%
mutate(group=case_when(연령 < 20 ~ 1,
                       20 <= 연령 & 연령 < 40 ~ 2,
                       40 <= 연령 & 연령 < 60 ~ 3,
                       60 <= 연령 ~ 4)) %>%
mutate(variant=case_when(variant==c("-") ~ c("conv"),
                         variant==c("알파형(α,영국변이)") ~ c("alpha"),
                         variant==c("델타형(δ,인도변이)") ~ c("delta"),
                         TRUE ~ c("etc"))) %>%
dplyr::select(group, variant, date)
options(warn=-1) 
read_excel("../data/db210913.xlsx") -> master_df
master_df %>% rename(reported=신고일, diagnosis=진단일, onset=발병일, published=보도일,
                     variant=변이, age=연령, route=감염경로) -> master_prop3

as.Date(master_prop3$reported) -> master_prop3$reported
as.Date(master_prop3$diagnosis) -> master_prop3$diagnosis
as.Date(master_prop3$onset) -> master_prop3$onset
as.Date(master_prop3$published) -> master_prop3$published

master_prop3 %<>% filter(diagnosis >= as.Date("2021-07-11"))

master_prop3$published <- strftime(master_prop3$published, "%V")

master_prop3 %<>% rename(date = published) %>%
mutate(group = case_when(age < 20 ~ "1",
                         age >= 20 & age < 40 ~ "2",
                         age >= 40 & age < 60 ~ "3",
                         age >= 60 ~ "4")) %>% 
filter(!is.na(variant)) %>% filter(variant!=c("분석불가")) %>%
mutate(variant=case_when(variant==c("-") ~ c("conv"),
                         variant==c("알파형") ~ c("alpha"),
                         variant==c("델타형") ~ c("delta"),
                         TRUE ~ c("etc"))) %>%
dplyr::select(group, variant, date)

rbind(master_prop, master_prop3) -> master_prop
as.data.frame(unique(master_prop$date)) -> prop_cal
colnames(prop_cal) <- c("date")

master_prop %>% group_by(group) %>% count(variant, date) -> master_group

prop_list <- list()

for (g in 1:(length(unique(master_prop$group)))){  

    master_prop %>% group_by(group) %>% count(variant, date) -> master_group
    
    master_group %<>% filter(group == g) %>% arrange(date)
    unique(master_group$date) %>% as.data.frame() -> group_cal
    colnames(group_cal) <- c("date")

    master_group %<>% spread(date, n) %>% t() %>% as.data.frame()
    colnames(master_group) <- master_group[2,]
    master_group[-(1:2),] -> master_group

    cbind(master_group, group_cal) -> master_group
    merge(master_group, prop_cal, by=c("date"), all.y=TRUE) -> master_group 
    master_group[is.na(master_group)] <- 0

    as.integer(master_group$alpha) -> master_group$alpha; as.integer(master_group$conv) -> master_group$conv
    as.integer(master_group$delta) -> master_group$delta; as.integer(master_group$etc) -> master_group$etc

    master_group %>%
    mutate('Wild-type' = conv/(conv+alpha+etc+delta),
           Alpha = alpha/(conv+alpha+etc+delta),
           Delta = delta/(conv+alpha+etc+delta),
           Etc = etc/(conv+alpha+etc+delta)) -> prop_list[[g]]
}
fig_temp <- list()

for (k in 1:(length(unique(master_prop$group)))){  

    prop_list[[k]] %>% dplyr::select(date, 'Wild-type', Alpha, Delta, Etc) %>% 
    gather(date, group) -> temp
    
    cbind(temp, rep(master_group$date, by=length(unique(master_group$group)))) -> temp
    colnames(temp) <- c("group", "prop", "week")
    temp$group <- factor(temp$group, levels=c('Delta','Alpha','Etc','Wild-type'))
    temp$week <- factor(temp$week, levels=unique(master_prop$date))
    temp -> fig_temp[[k]]
}<jupyter_output><empty_output><jupyter_text>Estimation of parameters<jupyter_code>time_delay <- 9

(start_date <- as.Date("2020-11-07") + time_delay)
(end_date <- as.Date("2020-11-18") + time_delay)
convert_lnorm <- function(mu, sd) {
    tmp <- log((sd / mu)^2 + 1)
    mulog <- log(mu) - 0.5 * tmp; sdlog <- sqrt(tmp)
    list(mulog = mulog, sdlog = sdlog)
}
#### incubation period
inc_fit = list(meanlog=convert_lnorm(4.91, 2.66)$mulog, sdlog=convert_lnorm(4.91, 2.66)$sdlog)
incubation <- function(t){plnorm(t, inc_fit$meanlog, inc_fit$sdlog) - plnorm(t-1, inc_fit$meanlog, inc_fit$sdlog)}


#### right-truncated time delay from illness onset to reporting
rep_fit = list(shape=param[1], scale=param[2])
onsettolabconf <- function(t){pweibull(t,  shape=rep_fit$shape, scale=rep_fit$scale) - 
                              pweibull(t-1, shape=rep_fit$shape, scale=rep_fit$scale)}

#### time delay from infection to reporting
infectiontoreport <- function(t){convolve(onsettolabconf(t),rev(incubation(t)),type = c("open"))}
infectiontoreport <- function(t){convolve(onsettolabconf(1:t),rev(incubation(1:t)),type = c("open"))}

#### generation time distribution
gi_fit = list(shape=2.305, scale=5.452)
generation <- function(t){pweibull(t, shape = gi_fit$shape, scale = gi_fit$scale) - 
                          pweibull(t-1, shape = gi_fit$shape, scale = gi_fit$scale)}
data_1 %<>% rename(group1 = n); data_2 %<>% rename(group2 = n);
data_3 %<>% rename(group3 = n); data_4 %<>% rename(group4 = n);

cbind(data_1, data_2, data_3, data_4) %>% dplyr::select(group1, group2, group3, group4) -> temp_group
cbind(data_1[,1], temp_group) -> master_df
colnames(master_df) <- c("Date", "group1", "group2", "group3", "group4")<jupyter_output><empty_output><jupyter_text>Initial conditions for the Model 1<jupyter_code>burn <- 21 

#### population size
popsize <- c(8932023, 13694461, 16783392, 11939384)

#### initial values for recovered cases
master_df %>% filter(Date >= (start_date-(burn-1))) %>% filter(Date <= end_date) %>% dplyr::select(-Date) ->  init_value

(master_df %>% filter(Date <= (start_date-(burn-1))) %>% dplyr::select(-Date) %>% 
colSums() %>% sum() %>% as.numeric() -> cum_cases)
(appx_cum_cases = 26000)

total_N = sum(init_value$group1)+sum(init_value$group2)+sum(init_value$group3)+sum(init_value$group4)
c(sum(init_value$group1)/total_N, sum(init_value$group2)/total_N,
  sum(init_value$group3)/total_N, sum(init_value$group4)/total_N) -> age_prop
round(appx_cum_cases*age_prop) -> init_recover<jupyter_output><empty_output><jupyter_text>Initial conditions for the Model 2<jupyter_code>(model2_date <- as.Date("2021-07-01")+time_delay)
last_week <- 27
#### vaccine effectivness by vaccines 
## AstraZenca with and without Delta effect
ve_AZ1 = 0.5; ve_AZ2 = 0.77
ve_AZ1_delta = 0.3; ve_AZ2_delta = 0.66

## J&J with and without Delta effect
ve_JJ = 0.5; ve_JJ_delta = 0.3

## mRNA with and without Delta effect
ve_mRNA1 = 0.5; ve_mRNA2 = 0.95
ve_mRNA1_delta = 0.35; ve_mRNA2_delta = 0.88

## vaccine effectivenss against hospitalization
ve_hosp = 0.95
read_excel("../data/vaxx_0926.xlsx") -> vac_prop
vac_prop %<>% rename(age = 연령, dose = 차수, date = 날짜, pfizer = 화이자, JJ = 얀센,
                    moderna = 모더나, AZ = 'AZ-AZ', AZPf = 'AZ-Pf(2차만)') %>%
mutate(mRNA = pfizer + moderna) %>% 
mutate(age = case_when(age==c('12-17세') ~ 1,
                       age==c('18-29세') ~ 2, age==c('30-39세') ~ 2,
                       age==c('40-49세') ~ 3, age==c('50-59세') ~ 3,
                       age==c('60-64세') ~ 4, age==c('65-69세') ~ 4, age==c('70-74세') ~ 4, age==c('75-79세') ~ 4,
                       age==c('80세이상') ~ 4)) %>%
mutate(week = strftime(date, "%V")) %>% dplyr::select(week, age, dose, mRNA, AZ, JJ, AZPf) %>%
mutate(dose = case_when(dose==1 ~ c("dose1"),
                        dose==2 ~ c("dose2")))

vac_prop %>% tail()
vac_prop %>% filter(age==1) -> vac_prop_1

vac_prop %>% dplyr::select(week, age, dose, mRNA) %>% filter(age==2) %>% 
group_by(week, dose) %>% summarise(mRNA = sum(mRNA)) %>% mutate(age=2) -> vac_prop_2_mRNA

vac_prop %>% dplyr::select(week, age, dose, AZ) %>% filter(age==2) %>% 
group_by(week, dose) %>% summarise(AZ = sum(AZ)) %>% mutate(age=2) -> vac_prop_2_AZ

vac_prop %>% dplyr::select(week, age, dose, JJ) %>% filter(age==2) %>% 
group_by(week, dose) %>% summarise(JJ = sum(JJ)) %>% mutate(age=2) -> vac_prop_2_JJ

vac_prop %>% dplyr::select(week, age, dose, AZPf) %>% filter(age==2) %>% 
group_by(week, dose) %>% summarise(AZPf = sum(AZPf)) %>% mutate(age=2) -> vac_prop_2_AZPf

Reduce(function(x,y) merge(x,y,by=c("week", "age", "dose"),all=TRUE),
       list(vac_prop_2_mRNA, vac_prop_2_AZ,vac_prop_2_JJ, vac_prop_2_AZPf)) -> vac_prop_2


vac_prop %>% dplyr::select(week, age, dose, mRNA) %>% filter(age==3) %>% 
group_by(week, dose) %>% summarise(mRNA = sum(mRNA)) %>% mutate(age=3) -> vac_prop_3_mRNA

vac_prop %>% dplyr::select(week, age, dose, AZ) %>% filter(age==3) %>% 
group_by(week, dose) %>% summarise(AZ = sum(AZ)) %>% mutate(age=3) -> vac_prop_3_AZ

vac_prop %>% dplyr::select(week, age, dose, JJ) %>% filter(age==3) %>% 
group_by(week, dose) %>% summarise(JJ = sum(JJ)) %>% mutate(age=3) -> vac_prop_3_JJ

vac_prop %>% dplyr::select(week, age, dose, AZPf) %>% filter(age==3) %>% 
group_by(week, dose) %>% summarise(AZPf = sum(AZPf)) %>% mutate(age=3) -> vac_prop_3_AZPf

Reduce(function(x,y) merge(x,y,by=c("week", "age", "dose"),all=TRUE),
       list(vac_prop_3_mRNA, vac_prop_3_AZ,vac_prop_3_JJ, vac_prop_3_AZPf)) -> vac_prop_3


vac_prop %>% dplyr::select(week, age, dose, mRNA) %>% filter(age==4) %>% 
group_by(week, dose) %>% summarise(mRNA = sum(mRNA)) %>% mutate(age=4) -> vac_prop_4_mRNA

vac_prop %>% dplyr::select(week, age, dose, AZ) %>% filter(age==4) %>% 
group_by(week, dose) %>% summarise(AZ = sum(AZ)) %>% mutate(age=4) -> vac_prop_4_AZ

vac_prop %>% dplyr::select(week, age, dose, JJ) %>% filter(age==4) %>% 
group_by(week, dose) %>% summarise(JJ = sum(JJ)) %>% mutate(age=4) -> vac_prop_4_JJ

vac_prop %>% dplyr::select(week, age, dose, AZPf) %>% filter(age==4) %>% 
group_by(week, dose) %>% summarise(AZPf = sum(AZPf)) %>% mutate(age=4) -> vac_prop_4_AZPf

Reduce(function(x,y) merge(x,y,by=c("week", "age", "dose"),all=TRUE),
       list(vac_prop_4_mRNA, vac_prop_4_AZ,vac_prop_4_JJ, vac_prop_4_AZPf)) -> vac_prop_4

rbind(vac_prop_1, vac_prop_2, vac_prop_3, vac_prop_4) -> vac_prop
vac_prop %<>% arrange(week, age, dose)

vac_prop %>% tail()
master_prop %>% count(variant, date) %>% arrange(date) -> master_group_total

unique(master_group_total$date) %>% as.data.frame() -> group_cal
colnames(group_cal) <- c("date")

master_group_total %<>% spread(date, n) %>% t() %>% as.data.frame()
colnames(master_group_total) <- master_group_total[1,]
master_group_total[-1,] -> master_group_total

cbind(master_group_total, group_cal) -> master_group_total
merge(master_group_total, prop_cal, by=c("date"), all.y=TRUE) -> master_group_total 
master_group_total[is.na(master_group_total)] <- 0

as.integer(master_group_total$alpha) -> master_group_total$alpha
as.integer(master_group_total$conv) -> master_group_total$conv
as.integer(master_group_total$delta) -> master_group_total$delta
as.integer(master_group_total$etc) -> master_group_total$etc

master_group_total %<>%
mutate('Wild-type' = conv/(conv+alpha+etc+delta),
       Alpha = alpha/(conv+alpha+etc+delta),
       Delta = delta/(conv+alpha+etc+delta),
       Etc = etc/(conv+alpha+etc+delta))

master_group_total %<>% filter(date %in% c('09','10','11','12','13','14','15','16','17','18','19',
                                           '20','21','22','23','24','25','26','27')) %>%
dplyr::select(date, 'Wild-type',Alpha, Delta, Etc) %>% rename(week = date) 

vac_prop -> vac_prop_temp
merge(vac_prop, master_group_total, by=c('week'), all.y=TRUE) -> vac_prop
merge(vac_prop_temp, master_group_total, by=c('week'), all=TRUE) -> vac_prop_all

vac_prop %>% head()
vac_prop_all %>% tail()
mRNA_list <- list(); AZ_list <- list(); JJ_list <- list()
age_select <- c(1, 2, 3, 4)


for(i in 1:length(age_select)){
    vac_prop %>% dplyr::select(week, age, dose, AZPf, Delta) %>% filter(age==i) %>% 
    spread(key=dose, AZPf) -> temp_AZPf
    vac_prop %>% dplyr::select(week, age, dose, mRNA, Delta) %>% filter(age==i) %>% 
    spread(key=dose, mRNA) -> temp_mRNA

    merge(temp_mRNA, temp_AZPf, by = c("week", "age", "Delta")) %>%
    rename(dose1 = dose1.x, dose2 = dose2.x, dose_AZPf = dose2.y) %>%
    mutate(prop_1 = (cumsum(dose1)-cumsum(dose2))/popsize[i], 
           prop_2 = (cumsum(dose2)+cumsum(dose_AZPf))/popsize[i]) %>% 
    mutate(prop_mRNA = prop_1*(ve_mRNA1*(1-Delta)+ve_mRNA1_delta*Delta) + 
           prop_2*(ve_mRNA2*(1-Delta)+ve_mRNA2_delta*Delta)) %>%
    dplyr::select(week, age, prop_mRNA) -> mRNA_list[[i]]
    
}

for(i in 1:length(age_select)){
    vac_prop %>% dplyr::select(week, age, dose, AZPf, Delta) %>% filter(age==i) %>% 
    spread(key=dose, AZPf) -> temp_AZPf
    vac_prop %>% dplyr::select(week, age, dose, AZ, Delta) %>% filter(age==i) %>% 
    spread(key=dose, AZ) -> temp_AZ

    merge(temp_AZ, temp_AZPf, by = c("week", "age", "Delta")) %>% 
    rename(dose1 = dose1.x, dose2 = dose2.x, dose_AZPf = dose2.y) %>%
    mutate(prop_1 = (cumsum(dose1)-cumsum(dose2)-cumsum(dose_AZPf))/popsize[i], 
           prop_2 = cumsum(dose2)/popsize[i]) %>% 
    mutate(prop_AZ = prop_1*(ve_AZ1*(1-Delta)+ve_AZ1_delta*Delta) + 
           prop_2*(ve_AZ2*(1-Delta)+ve_AZ2_delta*Delta)) %>%
    dplyr::select(week, age, prop_AZ) -> AZ_list[[i]] -> AZ_list[[i]]
    
}

for(i in 1:length(age_select)){
    vac_prop %>% dplyr::select(week, age, dose, JJ, Delta) %>% filter(age==i) %>% 
    spread(key=dose, JJ) %>%
    mutate(prop_1 = cumsum(dose1)/popsize[i]) %>%
    mutate(prop_JJ = prop_1*(ve_JJ*(1-Delta)+ve_JJ_delta*Delta)) %>%
    dplyr::select(week, age, prop_JJ) -> JJ_list[[i]]   
}
Reduce(function(x,y) merge(x,y,by=c("week", "age"),all=TRUE),
       list(mRNA_list[[1]], AZ_list[[1]], JJ_list[[1]])) %>% 
mutate(prop = prop_mRNA + prop_AZ + prop_JJ) %>% rename(prop1 = prop) %>% 
dplyr::select(week, prop1) -> vac_prop_group1

Reduce(function(x,y) merge(x,y,by=c("week", "age"),all=TRUE),
       list(mRNA_list[[2]], AZ_list[[2]], JJ_list[[2]])) %>% 
mutate(prop = prop_mRNA + prop_AZ + prop_JJ) %>% rename(prop2 = prop) %>% 
dplyr::select(week, prop2) -> vac_prop_group2

Reduce(function(x,y) merge(x,y,by=c("week", "age"),all=TRUE),
       list(mRNA_list[[3]], AZ_list[[3]], JJ_list[[3]])) %>% 
mutate(prop = prop_mRNA + prop_AZ + prop_JJ) %>% rename(prop3 = prop) %>% 
dplyr::select(week, prop3) -> vac_prop_group3
Reduce(function(x,y) merge(x,y,by=c("week", "age"),all=TRUE),
       list(mRNA_list[[4]], AZ_list[[4]], JJ_list[[4]])) %>% 
mutate(prop = prop_mRNA + prop_AZ + prop_JJ) %>% rename(prop4 = prop) %>% 
dplyr::select(week, prop4) -> vac_prop_group4

Reduce(function(x,y) merge(x,y,by=c("week"),all=TRUE),
       list(vac_prop_group1, vac_prop_group2, vac_prop_group3, vac_prop_group4)) -> vac_prop_total

vac_prop_total %>% t() -> temp
colnames(temp) <- temp[1,]
temp[-1,((24-8):(last_week-8))] -> vac_prop6

as.data.frame(vac_prop6)
Start_renewal <- (as.Date("2021-06-13")+time_delay) #### Start from week 21: week 24 - 20 days of burn-in period

master_df %>% filter(Date >= Start_renewal-(burn-1)) %>% 
mutate(Week = strftime(Date, "%V"))-> variant_df

groups <- c("group1", "group2", "group3", "group4")
variant_dff <- list()

for (g in 1:(length(unique(master_prop$group)))){
    
    variant_df %>% dplyr::select(Date, Week, groups[g]) -> variant_temp
    prop_list[[g]] %>% dplyr::select(date, 'Wild-type', Alpha, Delta, Etc) %>% 
    rename(Week = date, Wild='Wild-type') -> prop_temp
    merge(variant_temp, prop_temp, by=c("Week")) %>% rename(case=groups[g]) %>%
    mutate(exp_conv = case*Wild, exp_alpha = case*Alpha,
           exp_delta = case*Delta, exp_etc = case*Etc) %>%
    dplyr::select(Week, Date, case, exp_conv, exp_alpha, exp_delta, exp_etc) -> variant_dff[[g]]
    
    as.Date(variant_dff[[g]]$Date) -> variant_dff[[g]]$Date
    variant_dff[[g]] %<>% arrange(Date) %>% mutate(t = 0:(nrow(variant_dff[[g]])-1)) %>% 
    filter(Date <= model2_date)
}
variant_dff[[2]] %>% filter(Date <= model2_date) %>%
as.data.frame() %>% mutate(adjusted_Date = Date-time_delay,
                                                adjusted_Week = strftime(adjusted_Date, "%V")) %>%
mutate(adjusted_Week2 = as.numeric(adjusted_Week)) %>%
mutate(adjusted_Week2 = case_when(t<=20 ~ 24,
                                  TRUE~ adjusted_Week2)) %>%
mutate(vac_prop = as.numeric(vac_prop6[1,(adjusted_Week2-23)])) %>% dplyr::select(vac_prop) -> temp_vac1

variant_dff[[2]] %>% filter(Date <= model2_date) %>%
as.data.frame() %>% mutate(adjusted_Date = Date-time_delay,
                                                adjusted_Week = strftime(adjusted_Date, "%V")) %>%
mutate(adjusted_Week2 = as.numeric(adjusted_Week)) %>%
mutate(adjusted_Week2 = case_when(t<=20 ~ 24,
                                  TRUE~ adjusted_Week2)) %>%
mutate(vac_prop = as.numeric(vac_prop6[2,(adjusted_Week2-23)])) %>% dplyr::select(vac_prop) -> temp_vac2

variant_dff[[2]] %>% filter(Date <= model2_date) %>%
as.data.frame() %>% mutate(adjusted_Date = Date-time_delay,
                                                adjusted_Week = strftime(adjusted_Date, "%V")) %>%
mutate(adjusted_Week2 = as.numeric(adjusted_Week)) %>%
mutate(adjusted_Week2 = case_when(t<=20 ~ 24,
                                  TRUE~ adjusted_Week2)) %>%
mutate(vac_prop = as.numeric(vac_prop6[3,(adjusted_Week2-23)])) %>% dplyr::select(vac_prop) -> temp_vac3

variant_dff[[2]] %>% filter(Date <= model2_date) %>%
as.data.frame() %>% mutate(adjusted_Date = Date-time_delay,
                                                adjusted_Week = strftime(adjusted_Date, "%V")) %>%
mutate(adjusted_Week2 = as.numeric(adjusted_Week)) %>%
mutate(adjusted_Week2 = case_when(t<=20 ~ 24,
                                  TRUE~ adjusted_Week2)) %>%
mutate(vac_prop = as.numeric(vac_prop6[4,(adjusted_Week2-23)])) %>% dplyr::select(vac_prop) -> temp_vac4

cbind(temp_vac1, temp_vac2, temp_vac3, temp_vac4) -> vac_prop
#### age group version without back-projection
temp_conv <- list(); temp_alpha <- list(); temp_delta <- list(); temp_etc <- list()

for(k in 1:4){
    variant_dff[[k]] %>% arrange(Date) %>% filter(Date <= model2_date) %>% 
    dplyr::select(exp_conv) -> temp_conv[[k]]
    variant_dff[[k]] %>% arrange(Date) %>% filter(Date <= model2_date) %>%
    dplyr::select(exp_alpha) -> temp_alpha[[k]]
    variant_dff[[k]] %>% arrange(Date) %>% filter(Date <= model2_date) %>%
    dplyr::select(exp_delta) -> temp_delta[[k]]
    variant_dff[[k]] %>% arrange(Date) %>% filter(Date <= model2_date) %>%
    dplyr::select(exp_etc) -> temp_etc[[k]]}


do.call(cbind, temp_conv) -> master_conv
do.call(cbind, temp_alpha) -> master_alpha
do.call(cbind, temp_delta) -> master_delta
do.call(cbind, temp_etc) -> master_etc

colnames(master_conv) <- c("group1", "group2", "group3", "group4")
colnames(master_alpha) <- c("group1", "group2", "group3", "group4")
colnames(master_delta) <- c("group1", "group2", "group3", "group4")
colnames(master_etc) <- c("group1", "group2", "group3", "group4")
### fixed parameters for the base model
epsilon = 1/3
sigma = 1/5
omega = 1/18
IFR <- c(0.00001, 0.02, 0.21, 8.36)*0.01
SR <- c(0.001, 0.157, 1.071, 9.129)*0.01

(nrow(init_value) -> lgth)
loglike <- function(param){
    
    #### risk of death among severe cases
    lk_binom <- rep(0, 4)
    
    for (k in 1:4){dbinom(IFR[k]*10000000, size = SR[k]*10000000, prob = param[k], log=TRUE) -> lk_binom[k]}
    return(sum(-lk_binom))
}
    
#### Optim
param0 = c(0.01, 0.127388535031847, 0.196078431372549, 0.915762953225983)

options(warn=-1)
opt_FR <- optim(param0, fn=loglike, method=c("BFGS"), control=list(maxit=10000))
options(warn=0)  

IFR <- opt_FR$par; IFR[1] <- c(0.001*0.01)
time_org <- proc.time()    

loglike <- function(param){
    
    
    llk_total <- rep(0, lgth)


    #### SEIR compartments 
    S = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
    E = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
    I = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
    H = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
    D = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
    R = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))

    pSE = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
    pEI = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
    pIR = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
    pIH = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
    pHR = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
    pHD = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))

    Reported_modelled = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
    llk = data.frame(matrix(0, ncol=ncol(init_value), nrow=lgth))


    #### initail conditions for the calibration
    for (g in 1:ncol(init_value)){         
        S[1,g] <- popsize[g] - (init_value[1,g]*3 + init_value[2,g]*5 + init_recover[g])
        E[1,g] <- init_value[1,g]*3
        I[1,g] <- init_value[2,g]*5
        H[1,g] <- 0  
        D[1,g] <- 0
        R[1,g] <- init_recover[g]     
        
        pSE[1,g] <- 0           
        pEI[1,g] <- 0
        pIR[1,g] <- 0
        pIH[1,g] <- 0
        pHR[1,g] <- 0
        pHD[1,g] <- 0
    }


    #### Dynamics of SEIDR model
    for (t in 2:lgth){
        
        for (g in 1:ncol(init_value)){ 
      
            pSE[t,g] <- S[t-1,g]*(1-exp(-(param[g]*(
                adj_contact[g,1]*I[t-1,1]/(S[t-1,1]+E[t-1,1]+I[t-1,1]+R[t-1,1]) + 
                adj_contact[g,2]*I[t-1,2]/(S[t-1,2]+E[t-1,2]+I[t-1,2]+R[t-1,2]) +
                adj_contact[g,3]*I[t-1,3]/(S[t-1,3]+E[t-1,3]+I[t-1,3]+R[t-1,3]) + 
                adj_contact[g,4]*I[t-1,4]/(S[t-1,4]+E[t-1,4]+I[t-1,4]+R[t-1,4])))))

            pEI[t,g] <- E[t-1,g]*(1-exp(-epsilon))
            pIR[t,g] <- I[t-1,g]*(1-exp(-sigma*(1-SR[g])))
            pIH[t,g] <- I[t-1,g]*(1-exp(-sigma*(SR[g])))
            pHD[t,g] <- H[t-1,g]*(1-exp(-omega*IFR[g]))
            pHR[t,g] <- H[t-1,g]*(1-exp(-omega*(1-IFR[g])))


            S[t,g] <- S[t-1,g] - pSE[t,g] 
            E[t,g] <- E[t-1,g] + pSE[t,g] - pEI[t,g]
            I[t,g] <- I[t-1,g] + pEI[t,g] - pIH[t,g] - pIR[t,g] 
            H[t,g] <- H[t-1,g] + pIH[t,g] - pHD[t,g] - pHR[t,g] 
            D[t,g] <- D[t-1,g] + pHD[t,g] 
            R[t,g] <- R[t-1,g] + pHR[t,g] + pIR[t,g] 

            I[,g][I[,g]<0] <- 0
            pSE[,g][pSE[,g]<0] <- 0

            
            Report <- data.frame(matrix(0, ncol=ncol(init_value), nrow=t))  
            for (tau in 1:(t-1)){Report[tau,g] = (pSE[t-tau+1,g])*infectiontoreport(tau)[tau]}   
            sum(Report[,g]) -> Reported_modelled[t,g]

            llk[t,g] <- -(-Reported_modelled[t,g]+init_value[t,g]*log(Reported_modelled[t,g])-lgamma(init_value[t,g]+1))
        }  
        
        llk_total[t] = sum(llk[t,])

    }

    
    #### reconstructing the NGM
    temp <- matrix(NA, 4, 4)
    NGM <- matrix(NA, 4, 4)
    
    (param[1:4]*param[5]) -> sus_param
    for (m in 1:4){adj_contact[m,]*sus_param[m]-> temp[m,]}
    for (n in 1:4){temp[,n]*1/sigma -> NGM[,n]}
 
    llk_conv <- matrix(0, nrow=ncol(master_conv), ncol=max(variant_dff[[1]]$t)) 
    llk_alpha <- matrix(0, nrow=ncol(master_conv), ncol=max(variant_dff[[1]]$t)) 
    llk_delta <- matrix(0, nrow=ncol(master_conv), ncol=max(variant_dff[[1]]$t)) 
    llk_etc <- matrix(0, nrow=ncol(master_conv), ncol=max(variant_dff[[1]]$t)) 
    

    
    #### for the wild type  
    Cs <- matrix(0, nrow=ncol(master_conv), ncol=max(variant_dff[[1]]$t))  
    for (m in 1:4){        
  
        for (t in burn:max(variant_dff[[1]]$t)){
           
            Css <- rep(0, max(variant_dff[[1]]$t))     

            for (tau in 1:(t-1)){
            Css[tau] = ((variant_dff[[1]]$exp_conv[t-tau+1])*generation(tau)*NGM[m,1])+
                       ((variant_dff[[2]]$exp_conv[t-tau+1])*generation(tau)*NGM[m,2])+
                       ((variant_dff[[3]]$exp_conv[t-tau+1])*generation(tau)*NGM[m,3])+  
                       ((variant_dff[[4]]$exp_conv[t-tau+1])*generation(tau)*NGM[m,4])
            }      

            Cs[m,t] = sum(Css)*(1-vac_prop[t,m])  
            Cs[m,t][Cs[m,t]<=0] <- 1e-5
            llk_conv[m,t] = -Cs[m,t]+master_conv[t+1,m]*log(Cs[m,t])-lgamma(master_conv[t+1,m]+1)
        }
    }          
    llk_conv_total <- -sum(llk_conv)
    
    
    #### for the Alpha
    Cs <- matrix(0, nrow=ncol(master_conv), ncol=max(variant_dff[[1]]$t)) 
     for (m in 1:4){        
  
        for (t in burn:max(variant_dff[[1]]$t)){
           
            Css <- rep(0, max(variant_dff[[1]]$t))     

            for (tau in 1:(t-1)){
            Css[tau] = ((variant_dff[[1]]$exp_alpha[t-tau+1])*generation(tau)*NGM[m,1]*param[6])+
                       ((variant_dff[[2]]$exp_alpha[t-tau+1])*generation(tau)*NGM[m,2]*param[6])+
                       ((variant_dff[[3]]$exp_alpha[t-tau+1])*generation(tau)*NGM[m,3]*param[6])+  
                       ((variant_dff[[4]]$exp_alpha[t-tau+1])*generation(tau)*NGM[m,4]*param[6])
            }      

            Cs[m,t] = sum(Css)*(1-vac_prop[t,m])  
            Cs[m,t][Cs[m,t]<=0] <- 1e-5
            llk_alpha[m,t] = -Cs[m,t]+master_alpha[t+1,m]*log(Cs[m,t])-lgamma(master_alpha[t+1,m]+1)
        }
    }          
    llk_alpha_total <- -sum(llk_alpha)  
    

    #### for the Delta
    Cs <- matrix(0, nrow=ncol(master_conv), ncol=max(variant_dff[[1]]$t))  
     for (m in 1:4){        
  
        for (t in burn:max(variant_dff[[1]]$t)){
           
            Css <- rep(0, max(variant_dff[[1]]$t))  

            for (tau in 1:(t-1)){
            Css[tau] = ((variant_dff[[1]]$exp_delta[t-tau+1])*generation(tau)*NGM[m,1]*param[7])+
                       ((variant_dff[[2]]$exp_delta[t-tau+1])*generation(tau)*NGM[m,2]*param[7])+
                       ((variant_dff[[3]]$exp_delta[t-tau+1])*generation(tau)*NGM[m,3]*param[7])+  
                       ((variant_dff[[4]]$exp_delta[t-tau+1])*generation(tau)*NGM[m,4]*param[7])
            }      

            Cs[m,t] = sum(Css)*(1-vac_prop[t,m]) 
            Cs[Cs[m,t]<=0] <- 1e-5
            llk_delta[m,t] = -Cs[m,t]+master_delta[t+1,m]*log(Cs[m,t])-lgamma(master_delta[t+1,m]+1)
        }
    }          
    llk_delta_total <- -sum(llk_delta)   
    
    
    #### for Etc
    Cs <- matrix(0, nrow=ncol(master_conv), ncol=max(variant_dff[[1]]$t)) 
     for (m in 1:4){        
  
        for (t in burn:max(variant_dff[[1]]$t)){
           
            Css <- rep(0, max(variant_dff[[1]]$t))    

            for (tau in 1:(t-1)){
            Css[tau] = ((variant_dff[[1]]$exp_etc[t-tau+1])*generation(tau)*NGM[m,1]*param[8])+
                       ((variant_dff[[2]]$exp_etc[t-tau+1])*generation(tau)*NGM[m,2]*param[8])+
                       ((variant_dff[[3]]$exp_etc[t-tau+1])*generation(tau)*NGM[m,3]*param[8])+  
                       ((variant_dff[[4]]$exp_etc[t-tau+1])*generation(tau)*NGM[m,4]*param[8])
            }      

            Cs[m,t] = sum(Css)*(1-vac_prop[t,m]) 
            Cs[m,t][Cs[m,t]<=0] <- 1e-5
            llk_etc[m,t] = -Cs[m,t]+master_etc[t+1,m]*log(Cs[m,t])-lgamma(master_etc[t+1,m]+1)
        }
    }          
    llk_etc_total <- -sum(llk_etc)
    
    
    LK_variants <- (llk_conv_total + llk_alpha_total + llk_delta_total + llk_etc_total)
  
    
    return(sum(llk_total[burn:lgth])+LK_variants)
}  

#### Optim
param0 = c(0.00839345001544167, 0.0213549791943258, 0.0225932358342151, 0.0611640144713435,
           0.889137763882033, 0.878319042953008, 1.5475731853412, 0.202334874567077)

options(warn=-1)

opt_est <- optim(param0, fn=loglike, method=c("BFGS"), control=list(maxit=10000), hessian=TRUE)
options(warn=0)  

(proc.time()-time_org)

opt_est

opt_est$par -> est_param
opt_est$hessian -> est_param_hes
saveRDS(est_param,"estimated_ui_param.rds")
saveRDS(est_param_hes,"estimated_ui_CIs.rds")
niter = 5000
#### 95% confidence intervals
library(MASS)

hes_p <- est_param_hes
fisher_info<-solve(hes_p)

est_param_sample <- mvrnorm(n=niter, mu=est_param, Sigma=fisher_info, tol=1e-25, empirical=FALSE, EISPACK=FALSE)<jupyter_output><empty_output><jupyter_text>Reconstructed next generation matrix<jupyter_code>temp <- matrix(NA, 4, 4)
NGM <- matrix(NA, 4, 4)
temp_R <- rep(0, niter)
NGM_list <- list()

for (k in 1:niter){
    for (m in 1:4){
        adj_contact[m,]*est_param_sample[k,m] -> temp[m,]}

    for (n in 1:4){
        temp[,n]*1/sigma -> NGM[,n]}
        
        NGM -> NGM_list[[k]]

    (eigen(NGM)$values[1] -> temp_R[k])
}

#### estimated baseline R and its 95% CIs
c(quantile(temp_R,0.025), quantile(temp_R,0.5), quantile(temp_R,0.975))
#### 95% CIs for each element in the NGM
temp <- rep(0, niter)
NGM -> NGM_lower; NGM-> NGM_upper

for (m in 1:nrow(NGM)){
    for (n in 1:ncol(NGM)){
        for (g in 1:niter){
            NGM_list[[g]][m,n] -> temp[g]
            quantile(temp,0.025) -> NGM_lower[m,n]
            quantile(temp,0.975) -> NGM_upper[m,n]
        }
    }
}<jupyter_output><empty_output><jupyter_text>Numerical simulations  Scenarios & Parameters<jupyter_code>#### duration of simulation
Start_date = as.Date("2021-06-13")
Finish_date = as.Date("2021-12-31")
Inter_date = as.Date("2021-12-26")

#### intervention
Counter_date = as.Date("2021-07-12")
Lift_date1 = as.Date("2021-10-25")
Lift_date2 = as.Date("2021-11-01")
Lift_date3 = as.Date("2021-12-25")
    
#### transformation
counter_date = as.numeric(Counter_date-Start_date)+1
lift_date1 = as.numeric(Lift_date1-Start_date)+1
lift_date2 = as.numeric(Lift_date2-Start_date)+1
lift_date3 = as.numeric(Lift_date3-Start_date)+1
lift_dates <- c(lift_date1, lift_date2, lift_date3)

lgth <- as.numeric(Finish_date-Start_date)
week_inter <- as.numeric(Inter_date-Start_date)
#### Rt
R_base <- as.numeric(quantile(temp_R, 0.5))
Rt_counter = c(0.6, 0.65, 0.7)
#### initial values for R compartment
master_df %>% filter(Date <= as.Date("2021-06-13")) %>% 
dplyr::select(-Date) %>% colSums %>% as.vector() -> init_R

#### initial values for E and I compartments 
master_df %>% filter(Date == as.Date("2021-06-21")) %>% 
dplyr::select(-Date) %>% colSums -> temp 
temp*3 -> init_E
temp*5 -> init_I

sum(init_E)
sum(init_I)

#### initial values for Sv compartment (vaccinated & fully protected) 
round(vac_prop[1,]*popsize) -> init_Sv

as.data.frame(t(vac_prop6)) -> init_vac
init_vac[1,] -> init_vac
init_vac %<>% mutate(prop1 = round(as.numeric(prop1)*popsize[1]), 
                     prop2 = round(as.numeric(prop2)*popsize[2]), 
                     prop3 = round(as.numeric(prop3)*popsize[3]), 
                     prop4 = round(as.numeric(prop4)*popsize[4])) -> init_Sv
prop_list[[1]] %>% filter(date == c("25")) %>%
dplyr::select('Wild-type', Alpha, Delta, Etc) %>% as.vector() -> temp_age1

prop_list[[2]] %>% filter(date == c("25")) %>%
dplyr::select('Wild-type', Alpha, Delta, Etc) %>% as.vector() -> temp_age2

prop_list[[3]] %>% filter(date == c("25")) %>%
dplyr::select('Wild-type', Alpha, Delta, Etc) %>% as.vector() -> temp_age3

prop_list[[4]] %>% filter(date == c("25")) %>%
dplyr::select('Wild-type', Alpha, Delta, Etc) %>% as.vector() -> temp_age4

init_Ew <- as.data.frame(c(init_E[1]*temp_age1[1], init_E[2]*temp_age2[1], 
                           init_E[3]*temp_age3[1], init_E[4]*temp_age4[1]))
init_Ea <- as.data.frame(c(init_E[1]*temp_age1[2], init_E[2]*temp_age2[2], 
                           init_E[3]*temp_age3[2], init_E[4]*temp_age4[2]))
init_Ed <- as.data.frame(c(init_E[1]*temp_age1[3], init_E[2]*temp_age2[3], 
                           init_E[3]*temp_age3[3], init_E[4]*temp_age4[3]))
init_Ee <- as.data.frame(c(init_E[1]*temp_age1[4], init_E[2]*temp_age2[4], 
                           init_E[3]*temp_age3[4], init_E[4]*temp_age4[4]))

init_Iw <- as.data.frame(c(init_I[1]*temp_age1[1], init_I[2]*temp_age2[1], 
                           init_I[3]*temp_age3[1], init_I[4]*temp_age4[1]))
init_Ia <- as.data.frame(c(init_I[1]*temp_age1[2], init_I[2]*temp_age2[2], 
                           init_I[3]*temp_age3[2], init_I[4]*temp_age4[2]))
init_Id <- as.data.frame(c(init_I[1]*temp_age1[3], init_I[2]*temp_age2[3], 
                           init_I[3]*temp_age3[3], init_I[4]*temp_age4[3]))
init_Ie <- as.data.frame(c(init_I[1]*temp_age1[4], init_I[2]*temp_age2[4], 
                           init_I[3]*temp_age3[4], init_I[4]*temp_age4[4]))

init_Ew <- round(init_Ew); init_Ea <- round(init_Ea); init_Ed <- round(init_Ed); init_Ee <- round(init_Ee)
init_Iw <- round(init_Iw); init_Ia <- round(init_Ia); init_Id <- round(init_Id); init_Ie <- round(init_Ie)
vac_scenarios <- function(param){
    
    #### scenario with the increase in allocated vaccines from week 40
    sc1 = 1.1
    sc2 = 1.2
    
    sc1_inc = param[1]
    sc2_inc = param[2]

    start_wk = 40
    start_wk_dose2 = 44


    ####
    vac_prop_all %>% mutate(mRNA_sc1 = case_when(week < start_wk ~ mRNA,
                                                   week >= start_wk & dose == c("dose1") ~ mRNA*sc1_inc,
                                                   week >= start_wk & week < start_wk_dose2 & dose == c("dose2") ~ mRNA,
                                                   week >= start_wk_dose2 & dose == c("dose2") ~ mRNA*sc1_inc)) %>%
    dplyr::select(-c("mRNA")) %>% rename(mRNA = mRNA_sc1)-> vac_prop_all_sc1


    vac_prop_all %>% mutate(mRNA_sc2 = case_when(week < start_wk ~ mRNA,
                                                   week >= start_wk & dose == c("dose1") ~ mRNA*sc2_inc,
                                                   week >= start_wk & week < start_wk_dose2 & dose == c("dose2") ~ mRNA,
                                                   week >= start_wk_dose2 & dose == c("dose2") ~ mRNA*sc2_inc)) %>%
    dplyr::select(-c("mRNA")) %>% rename(mRNA = mRNA_sc2)-> vac_prop_all_sc2


    vac_prop_all %>% mutate(mRNA_total = cumsum(mRNA)) %>% dplyr::select(mRNA_total) -> temp
    temp[nrow(temp),1] -> temp1
    vac_prop_all_sc1 %>% mutate(mRNA_total = cumsum(mRNA)) %>% dplyr::select(mRNA_total) -> temp
    temp[nrow(temp),1] -> temp2
    vac_prop_all_sc2 %>% mutate(mRNA_total = cumsum(mRNA)) %>% dplyr::select(mRNA_total) -> temp
    temp[nrow(temp),1] -> temp3
    
    return(sqrt(mean((temp2/temp1- sc1)^2)) + sqrt(mean((temp3/temp1- sc2)^2)))
}


## Optim
param0 = c(1.7, 2.5)

options(warn=-1)
opt_sc <- optim(param0, fn=vac_scenarios, method=c("BFGS"), control=list(maxit=10000))
options(warn=0)  

opt_sc
sc1_inc = opt_sc$par[1]
sc2_inc = opt_sc$par[2]

start_wk = 40
start_wk_dose2 = 44


####
vac_prop_all %>% mutate(mRNA_sc1 = case_when(week < start_wk ~ mRNA,
                                               week >= start_wk & dose == c("dose1") ~ mRNA*sc1_inc,
                                               week >= start_wk & week < start_wk_dose2 & dose == c("dose2") ~ mRNA,
                                               week >= start_wk_dose2 & dose == c("dose2") ~ mRNA*sc1_inc)) %>%
dplyr::select(-c("mRNA")) %>% rename(mRNA = mRNA_sc1)-> vac_prop_all_sc1


vac_prop_all %>% mutate(mRNA_sc2 = case_when(week < start_wk ~ mRNA,
                                               week >= start_wk & dose == c("dose1") ~ mRNA*sc2_inc,
                                               week >= start_wk & week < start_wk_dose2 & dose == c("dose2") ~ mRNA,
                                               week >= start_wk_dose2 & dose == c("dose2") ~ mRNA*sc2_inc)) %>%
dplyr::select(-c("mRNA")) %>% rename(mRNA = mRNA_sc2)-> vac_prop_all_sc2
vac_prop_all %>% mutate(week_num = as.numeric(week)) %>% dplyr::select(-week) %>% 
rename(week = week_num) %>% dplyr::select(-c("Wild-type", "Alpha", "Etc")) -> vac_prop_model
vac_prop_model -> vac_prop_time

vac_prop_all_sc1 %>% mutate(week_num = as.numeric(week)) %>% dplyr::select(-week) %>% 
rename(week = week_num) %>% dplyr::select(-c("Wild-type", "Alpha", "Etc")) -> vac_prop_model_sc1
vac_prop_model_sc1 -> vac_prop_time_sc1

vac_prop_all_sc2 %>% mutate(week_num = as.numeric(week)) %>% dplyr::select(-week) %>% 
rename(week = week_num) %>% dplyr::select(-c("Wild-type", "Alpha", "Etc")) -> vac_prop_model_sc2
vac_prop_model_sc2 -> vac_prop_time_sc2

(vac_prop_model %>% filter(week == 23) %>% dplyr::select(Delta) %>% unique() %>% as.numeric() -> init_Delta)<jupyter_output><empty_output><jupyter_text>Projection of future waves<jupyter_code>myCluster <- makeCluster(ncore-1, type = "PSOCK")
registerDoParallel(myCluster)

vac_case_final <- list(); vac_severe_final <- list(); vac_deaths_final <- list()
simul_results_list <- list()


simul_results_list <- foreach (k = 1:niter, .packages=c("dplyr","tidyr")) %dopar% {
    
    param <- est_param_sample[k,]
    
    for (p in 1:length(Rt_counter)){
        
        ############################################################
        Rt_counters <- Rt_counter[p]     
        counter_start <- counter_date 
        lift_start <- lift_dates[2]
        ############################################################
        
        #### fixed parameters
        epsilon = 1/3
        sigma = 1/5
        omega = 1/18
        IFR <- opt_FR$par; IFR[1] <- c(0.001*0.01)
        SR_org <- c(0.001, 0.157, 1.071, 9.129)*0.01
        ve_H = 0.95

        lgth <- as.numeric(Finish_date-Start_date)
        llk_total <- rep(0, lgth)
        
        Reported_modelled = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        Reported_modelled_w = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        Reported_modelled_a = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        Reported_modelled_d = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        Reported_modelled_e = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        
        #### compartments 
        S = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        
        Ew = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        Iw = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        Ea = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        Ia = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        Ed = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        Id = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        Ee = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        Ie = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        
        H = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        D = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        R = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))

        Sv = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))      
        eff_N = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth)) 
             
        
        #### rates 
        pSEw = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        pSEa = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        pSEd = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        pSEe = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        
        pEIw = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        pEIa = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        pEId = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        pEIe = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        
        pIRw = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        pIRa = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        pIRd = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        pIRe = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        
        pIHw = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        pIHa = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        pIHd = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        pIHe = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        
        pHR = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        pHD = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
    
        pvSS = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        pvSS_rev = data.frame(matrix(NA, ncol=ncol(init_value), nrow=lgth))
        prop_delta = data.frame(matrix(NA, ncol=1, nrow=lgth))
        
        
        #### initial conditions
        for (g in 1:ncol(init_value)){         
            S[1,g] <- popsize[g] - (init_Ew[g]+init_Iw[g]+init_Ea[g]+init_Ia[g]+
                                    init_Ed[g]+init_Id[g]+init_Ee[g]+init_Ie[g]+
                                    init_R[g]+init_Sv[g])
            
            Ew[1,g] <- init_Ew[g]; Iw[1,g] <- init_Iw[g]
            Ea[1,g] <- init_Ea[g]; Ia[1,g] <- init_Ia[g]         
            Ed[1,g] <- init_Ed[g]; Id[1,g] <- init_Id[g]  
            Ee[1,g] <- init_Ee[g]; Ie[1,g] <- init_Ie[g]         
            
            H[1,g] <- 0
            D[1,g] <- 0 
            R[1,g] <- init_R[g]  
            Sv[1,g] <- init_Sv[g]  
            
            pSEw[1,g] <- 0; pSEa[1,g] <- 0; pSEd[1,g] <- 0; pSEe[1,g] <- 0; 
            
            pEIw[1,g] <- 0; pEIa[1,g] <- 0; pEId[1,g] <- 0; pEIe[1,g] <- 0; 
            pIRw[1,g] <- 0; pIRa[1,g] <- 0; pIRd[1,g] <- 0; pIRe[1,g] <- 0;
            pIHw[1,g] <- 0; pIHa[1,g] <- 0; pIHd[1,g] <- 0; pIHe[1,g] <- 0;          
            pHR[1,g] <- 0; pHD[1,g] <- 0
                 
            pvSS[1,g] <- 0; pvSS_rev[1,g] <- 0;
            
            eff_N[1,g] = S[1,g]+Ew[1,g]+Iw[1,g]+Ea[1,g]+Ia[1,g]+Ed[1,g]+Id[1,g]+Ee[1,g]+Ie[1,g]+R[1,g]+Sv[1,g]
            
            Reported_modelled[1,g] <- 0
            Reported_modelled_w[1,g] <- 0; Reported_modelled_a[1,g] <- 0
            Reported_modelled_d[1,g] <- 0; Reported_modelled_e[1,g] <- 0
            
            prop_delta[1,1] <- init_Delta
        }


        #### Dynamics in rates
        for (t in 2:lgth){
    
            if(t < week_inter){
            #### Vaccinated individuals by the proportion of Delta
            t_mod <- as.integer(t/7)

            if(t_mod <= 1){
                vac_prop_all %>% filter(week==(24+t_mod)) %>% dplyr::select(Delta) %>% 
                                        unique() %>% as.numeric() -> prop_delta[t,1]
            } else {
                prop_delta[t,1] <- sum(Id[(t-1),])/(sum(Iw[(t-1),])+sum(Ia[(t-1),])+sum(Id[(t-1),])+sum(Ie[(t-1),]))
            }

            vac_prop_time %>% dplyr::select(week, Delta) %>% rename(Delta_prev = Delta) %>% distinct()-> vac_prop_model_prev
            merge(vac_prop_model, vac_prop_model_prev, by=c("week"), all.x=TRUE) -> vac_prop_model_temp

            vac_prop_model_temp %>% mutate(Delta = case_when(week==(24+t_mod) ~ as.numeric(prop_delta[t,1]),
                                                             week!=(24+t_mod) ~ as.numeric(Delta_prev))) %>% 
            filter(week <= (24+t_mod)) %>% dplyr::select(-Delta_prev) -> vac_prop_time

            mRNA_list <- list(); AZ_list <- list(); JJ_list <- list()
            for(i in 1:length(age_select)){
                vac_prop_time %>% dplyr::select(week, age, dose, AZPf, Delta) %>% filter(age==i) %>% 
                spread(key=dose, AZPf) -> temp_AZPf
                vac_prop_time %>% dplyr::select(week, age, dose, mRNA, Delta) %>% filter(age==i) %>% 
                spread(key=dose, mRNA) -> temp_mRNA

                merge(temp_mRNA, temp_AZPf, by = c("week", "age", "Delta")) %>%
                rename(dose1 = dose1.x, dose2 = dose2.x, dose_AZPf = dose2.y) %>%
                mutate(prop_1 = (cumsum(dose1)-cumsum(dose2))/popsize[i], 
                       prop_2 = (cumsum(dose2)+cumsum(dose_AZPf))/popsize[i]) %>% 
                mutate(prop_mRNA = prop_1*(ve_mRNA1*(1-Delta)+ve_mRNA1_delta*Delta) + 
                       prop_2*(ve_mRNA2*(1-Delta)+ve_mRNA2_delta*Delta)) %>%
                dplyr::select(week, age, prop_mRNA) -> mRNA_list[[i]]

            }

            for(i in 1:length(age_select)){
                vac_prop_time %>% dplyr::select(week, age, dose, AZPf, Delta) %>% filter(age==i) %>% 
                spread(key=dose, AZPf) -> temp_AZPf
                vac_prop_time %>% dplyr::select(week, age, dose, AZ, Delta) %>% filter(age==i) %>% 
                spread(key=dose, AZ) -> temp_AZ

                merge(temp_AZ, temp_AZPf, by = c("week", "age", "Delta")) %>% 
                rename(dose1 = dose1.x, dose2 = dose2.x, dose_AZPf = dose2.y) %>%
                mutate(prop_1 = (cumsum(dose1)-cumsum(dose2)-cumsum(dose_AZPf))/popsize[i], 
                       prop_2 = cumsum(dose2)/popsize[i]) %>% 
                mutate(prop_AZ = prop_1*(ve_AZ1*(1-Delta)+ve_AZ1_delta*Delta) + 
                       prop_2*(ve_AZ2*(1-Delta)+ve_AZ2_delta*Delta)) %>%
                dplyr::select(week, age, prop_AZ) -> AZ_list[[i]] -> AZ_list[[i]]

            }

            for(i in 1:length(age_select)){
                vac_prop_time %>% dplyr::select(week, age, dose, JJ, Delta) %>% filter(age==i) %>% 
                spread(key=dose, JJ) %>%
                mutate(prop_1 = cumsum(dose1)/popsize[i]) %>%
                mutate(prop_JJ = prop_1*(ve_JJ*(1-Delta)+ve_JJ_delta*Delta)) %>%
                dplyr::select(week, age, prop_JJ) -> JJ_list[[i]]   
            }

            Reduce(function(x,y) merge(x,y,by=c("week", "age"),all=TRUE),
                   list(mRNA_list[[1]], AZ_list[[1]], JJ_list[[1]])) %>% 
            mutate(prop = prop_mRNA + prop_AZ + prop_JJ) %>% rename(prop1 = prop) %>% 
            dplyr::select(week, prop1) -> vac_prop_group1

            Reduce(function(x,y) merge(x,y,by=c("week", "age"),all=TRUE),
                   list(mRNA_list[[2]], AZ_list[[2]], JJ_list[[2]])) %>% 
            mutate(prop = prop_mRNA + prop_AZ + prop_JJ) %>% rename(prop2 = prop) %>% 
            dplyr::select(week, prop2) -> vac_prop_group2

            Reduce(function(x,y) merge(x,y,by=c("week", "age"),all=TRUE),
                   list(mRNA_list[[3]], AZ_list[[3]], JJ_list[[3]])) %>% 
            mutate(prop = prop_mRNA + prop_AZ + prop_JJ) %>% rename(prop3 = prop) %>% 
            dplyr::select(week, prop3) -> vac_prop_group3
            Reduce(function(x,y) merge(x,y,by=c("week", "age"),all=TRUE),
                   list(mRNA_list[[4]], AZ_list[[4]], JJ_list[[4]])) %>% 
            mutate(prop = prop_mRNA + prop_AZ + prop_JJ) %>% rename(prop4 = prop) %>% 
            dplyr::select(week, prop4) -> vac_prop_group4

            Reduce(function(x,y) merge(x,y,by=c("week"),all=TRUE),
                   list(vac_prop_group1, vac_prop_group2, vac_prop_group3, vac_prop_group4)) -> vac_prop_total

            vac_prop_total %>% t() -> temp
            colnames(temp) <- temp[1,]
            temp[-1,(23+t_mod-8):(24+t_mod-8)] -> vac_prop_timet
            
            #### vaccination plans
            as.data.frame(t(vac_prop_timet)) -> diff_vac_org
            diff_vac_org$prop1 <- as.numeric(diff_vac_org$prop1); diff_vac_org$prop2 <- as.numeric(diff_vac_org$prop2)
            diff_vac_org$prop3 <- as.numeric(diff_vac_org$prop3); diff_vac_org$prop4 <- as.numeric(diff_vac_org$prop4)
            cbind(as.data.frame(diff(diff_vac_org$prop1)), as.data.frame(diff(diff_vac_org$prop2)),
                  as.data.frame(diff(diff_vac_org$prop3)), as.data.frame(diff(diff_vac_org$prop4))) -> diff_vac_org
            colnames(diff_vac_org) <- c("diff_1","diff_2","diff_3","diff_4")

            rbind(diff_vac_org, diff_vac_org[nrow(diff_vac_org),]) -> diff_vac_org
            diff_vac_org[1,] -> diff_vac_timet

            
            mRNA_list_org <- list(); AZ_list_org <- list(); JJ_list_org <- list()
            for(i in 1:length(age_select)){
                vac_prop_time %>% dplyr::select(week, age, dose, AZPf, Delta) %>% filter(age==i) %>% 
                spread(key=dose, AZPf) -> temp_AZPf
                vac_prop_time %>% dplyr::select(week, age, dose, mRNA, Delta) %>% filter(age==i) %>% 
                spread(key=dose, mRNA) -> temp_mRNA

                merge(temp_mRNA, temp_AZPf, by = c("week", "age", "Delta")) %>%
                rename(dose1 = dose1.x, dose2 = dose2.x, dose_AZPf = dose2.y) %>%
                mutate(prop_1 = (cumsum(dose1)-cumsum(dose2))/popsize[i], 
                       prop_2 = (cumsum(dose2)+cumsum(dose_AZPf))/popsize[i]) %>% 
                mutate(prop_mRNA = prop_1 + prop_2) %>%
                dplyr::select(week, age, prop_mRNA) -> mRNA_list_org[[i]]
            }

            for(i in 1:length(age_select)){
                vac_prop_time %>% dplyr::select(week, age, dose, AZPf, Delta) %>% filter(age==i) %>% 
                spread(key=dose, AZPf) -> temp_AZPf
                vac_prop_time %>% dplyr::select(week, age, dose, AZ, Delta) %>% filter(age==i) %>% 
                spread(key=dose, AZ) -> temp_AZ

                merge(temp_AZ, temp_AZPf, by = c("week", "age", "Delta")) %>% 
                rename(dose1 = dose1.x, dose2 = dose2.x, dose_AZPf = dose2.y) %>%
                mutate(prop_1 = (cumsum(dose1)-cumsum(dose2)-cumsum(dose_AZPf))/popsize[i], 
                       prop_2 = cumsum(dose2)/popsize[i]) %>% 
                mutate(prop_AZ = prop_1 + prop_2) %>%
                dplyr::select(week, age, prop_AZ) -> AZ_list_org[[i]]
            }

            for(i in 1:length(age_select)){
                vac_prop_time %>% dplyr::select(week, age, dose, JJ, Delta) %>% filter(age==i) %>% 
                spread(key=dose, JJ) %>%
                mutate(prop_1 = cumsum(dose1)/popsize[i]) %>%
                mutate(prop_JJ = prop_1) %>%
                dplyr::select(week, age, prop_JJ) -> JJ_list_org[[i]]   
            }

            Reduce(function(x,y) merge(x,y,by=c("week", "age"),all=TRUE),
                   list(mRNA_list_org[[1]], AZ_list_org[[1]], JJ_list_org[[1]])) %>% 
            mutate(prop = prop_mRNA + prop_AZ + prop_JJ) %>% rename(prop1 = prop) %>% 
            dplyr::select(week, prop1) -> vac_prop_group1

            Reduce(function(x,y) merge(x,y,by=c("week", "age"),all=TRUE),
                   list(mRNA_list_org[[2]], AZ_list_org[[2]], JJ_list_org[[2]])) %>% 
            mutate(prop = prop_mRNA + prop_AZ + prop_JJ) %>% rename(prop2 = prop) %>% 
            dplyr::select(week, prop2) -> vac_prop_group2

            Reduce(function(x,y) merge(x,y,by=c("week", "age"),all=TRUE),
                   list(mRNA_list_org[[3]], AZ_list_org[[3]], JJ_list_org[[3]])) %>% 
            mutate(prop = prop_mRNA + prop_AZ + prop_JJ) %>% rename(prop3 = prop) %>% 
            dplyr::select(week, prop3) -> vac_prop_group3

            Reduce(function(x,y) merge(x,y,by=c("week", "age"),all=TRUE),
                   list(mRNA_list_org[[4]], AZ_list_org[[4]], JJ_list_org[[4]])) %>% 
            mutate(prop = prop_mRNA + prop_AZ + prop_JJ) %>% rename(prop4 = prop) %>% 
            dplyr::select(week, prop4) -> vac_prop_group4

            Reduce(function(x,y) merge(x,y,by=c("week"),all=TRUE),
                   list(vac_prop_group1, vac_prop_group2, vac_prop_group3, vac_prop_group4)) -> vac_prop_total_org

            vac_prop_total_org %>% t() -> temp_org
            colnames(temp_org) <- temp_org[1,]
            temp_org[-1,(23+t_mod-8):(24+t_mod-8)] -> vac_prop_timet_org
            vac_prop_timet_org-vac_prop_timet -> vac_prop_np
            as.data.frame(vac_prop_np[,1]) -> prop_np
            
            
            #### adjusting the SR
            SR <- c(-(SR_org[1]-SR_org[1]*(1-ve_H))/1*prop_np[1,1]+SR_org[1],
                    -(SR_org[2]-SR_org[2]*(1-ve_H))/1*prop_np[2,1]+SR_org[2],
                    -(SR_org[3]-SR_org[3]*(1-ve_H))/1*prop_np[3,1]+SR_org[3],
                    -(SR_org[4]-SR_org[4]*(1-ve_H))/1*prop_np[4,1]+SR_org[4])
            }
            else{
                diff_vac_timet <- as.data.frame(c(0,0,0,0))
            }
            
        
            for (g in 1:ncol(init_value)){                                
                
            #### before 09 July 2021 
                if(t < counter_start){
                    pSEw[t,g] <- S[t-1,g]*(1-exp(-(param[g]*param[5]*(
                    adj_contact[g,1]*(Iw[t-1,1])/eff_N[t-1,1] + 
                    adj_contact[g,2]*(Iw[t-1,2])/eff_N[t-1,2] +
                    adj_contact[g,3]*(Iw[t-1,3])/eff_N[t-1,3] +
                    adj_contact[g,4]*(Iw[t-1,4])/eff_N[t-1,4])))) 

                    pSEa[t,g] <- S[t-1,g]*(1-exp(-(param[g]*param[5]*param[6]*(
                    adj_contact[g,1]*(Ia[t-1,1])/eff_N[t-1,1] + 
                    adj_contact[g,2]*(Ia[t-1,2])/eff_N[t-1,2] +
                    adj_contact[g,3]*(Ia[t-1,3])/eff_N[t-1,3] +
                    adj_contact[g,4]*(Ia[t-1,4])/eff_N[t-1,4]))))                   
                    

                    pSEd[t,g] <- S[t-1,g]*(1-exp(-(param[g]*param[5]*param[7]*(
                    adj_contact[g,1]*(Id[t-1,1])/eff_N[t-1,1] + 
                    adj_contact[g,2]*(Id[t-1,2])/eff_N[t-1,2] +
                    adj_contact[g,3]*(Id[t-1,3])/eff_N[t-1,3] +
                    adj_contact[g,4]*(Id[t-1,4])/eff_N[t-1,4])))) 

                
                    pSEe[t,g] <- S[t-1,g]*(1-exp(-(param[g]*param[5]*param[8]*(
                    adj_contact[g,1]*(Ie[t-1,1])/eff_N[t-1,1] + 
                    adj_contact[g,2]*(Ie[t-1,2])/eff_N[t-1,2] +
                    adj_contact[g,3]*(Ie[t-1,3])/eff_N[t-1,3] +
                    adj_contact[g,4]*(Ie[t-1,4])/eff_N[t-1,4])))) 
                }
                  
                
            #### period with countermeasures (7/12~)
                else if(t >= counter_start & t <= lift_start){
                    
                    t_mod <- as.integer((t-counter_start)/7)
                    Rt_counters_mod <- Rt_counters+0.005*t_mod
                    
                    pSEw[t,g] <- S[t-1,g]*(1-exp(-(param[g]*param[5]*Rt_counters_mod*(
                    adj_contact[g,1]*(Iw[t-1,1])/eff_N[t-1,1] + 
                    adj_contact[g,2]*(Iw[t-1,2])/eff_N[t-1,2] +
                    adj_contact[g,3]*(Iw[t-1,3])/eff_N[t-1,3] +
                    adj_contact[g,4]*(Iw[t-1,4])/eff_N[t-1,4])))) 

                    pSEa[t,g] <- S[t-1,g]*(1-exp(-(param[g]*param[5]*param[6]*Rt_counters_mod*(
                    adj_contact[g,1]*(Ia[t-1,1])/eff_N[t-1,1] + 
                    adj_contact[g,2]*(Ia[t-1,2])/eff_N[t-1,2] +
                    adj_contact[g,3]*(Ia[t-1,3])/eff_N[t-1,3] +
                    adj_contact[g,4]*(Ia[t-1,4])/eff_N[t-1,4]))))                   
                    

                    pSEd[t,g] <- S[t-1,g]*(1-exp(-(param[g]*param[5]*param[7]*Rt_counters_mod*(
                    adj_contact[g,1]*(Id[t-1,1])/eff_N[t-1,1] + 
                    adj_contact[g,2]*(Id[t-1,2])/eff_N[t-1,2] +
                    adj_contact[g,3]*(Id[t-1,3])/eff_N[t-1,3] +
                    adj_contact[g,4]*(Id[t-1,4])/eff_N[t-1,4])))) 

                
                    pSEe[t,g] <- S[t-1,g]*(1-exp(-(param[g]*param[5]*param[8]*Rt_counters_mod*(
                    adj_contact[g,1]*(Ie[t-1,1])/eff_N[t-1,1] + 
                    adj_contact[g,2]*(Ie[t-1,2])/eff_N[t-1,2] +
                    adj_contact[g,3]*(Ie[t-1,3])/eff_N[t-1,3] +
                    adj_contact[g,4]*(Ie[t-1,4])/eff_N[t-1,4])))) 
                } 
                    
                else {
                    pSEw[t,g] <- S[t-1,g]*(1-exp(-(param[g]*param[5]*(
                    adj_contact[g,1]*(Iw[t-1,1])/eff_N[t-1,1] + 
                    adj_contact[g,2]*(Iw[t-1,2])/eff_N[t-1,2] +
                    adj_contact[g,3]*(Iw[t-1,3])/eff_N[t-1,3] +
                    adj_contact[g,4]*(Iw[t-1,4])/eff_N[t-1,4])))) 

                    pSEa[t,g] <- S[t-1,g]*(1-exp(-(param[g]*param[5]*param[6]*(
                    adj_contact[g,1]*(Ia[t-1,1])/eff_N[t-1,1] + 
                    adj_contact[g,2]*(Ia[t-1,2])/eff_N[t-1,2] +
                    adj_contact[g,3]*(Ia[t-1,3])/eff_N[t-1,3] +
                    adj_contact[g,4]*(Ia[t-1,4])/eff_N[t-1,4]))))                   
                    

                    pSEd[t,g] <- S[t-1,g]*(1-exp(-(param[g]*param[5]*param[7]*(
                    adj_contact[g,1]*(Id[t-1,1])/eff_N[t-1,1] + 
                    adj_contact[g,2]*(Id[t-1,2])/eff_N[t-1,2] +
                    adj_contact[g,3]*(Id[t-1,3])/eff_N[t-1,3] +
                    adj_contact[g,4]*(Id[t-1,4])/eff_N[t-1,4])))) 

                
                    pSEe[t,g] <- S[t-1,g]*(1-exp(-(param[g]*param[5]*param[8]*(
                    adj_contact[g,1]*(Ie[t-1,1])/eff_N[t-1,1] + 
                    adj_contact[g,2]*(Ie[t-1,2])/eff_N[t-1,2] +
                    adj_contact[g,3]*(Ie[t-1,3])/eff_N[t-1,3] +
                    adj_contact[g,4]*(Ie[t-1,4])/eff_N[t-1,4]))))  
                }   
                    
            pEIw[t,g] <- Ew[t-1,g]*(1-exp(-epsilon))
            pEIa[t,g] <- Ea[t-1,g]*(1-exp(-epsilon))
            pEId[t,g] <- Ed[t-1,g]*(1-exp(-epsilon))
            pEIe[t,g] <- Ee[t-1,g]*(1-exp(-epsilon))
                    
            pIRw[t,g] <- Iw[t-1,g]*(1-exp(-sigma*(1-SR[g])))
            pIRa[t,g] <- Ia[t-1,g]*(1-exp(-sigma*(1-SR[g])))
            pIRd[t,g] <- Id[t-1,g]*(1-exp(-sigma*(1-SR[g])))
            pIRe[t,g] <- Ie[t-1,g]*(1-exp(-sigma*(1-SR[g])))
                    
            pIHw[t,g] <- Iw[t-1,g]*(1-exp(-sigma*(SR[g])))
            pIHa[t,g] <- Ia[t-1,g]*(1-exp(-sigma*(SR[g])))
            pIHd[t,g] <- Id[t-1,g]*(1-exp(-sigma*(SR[g])))
            pIHe[t,g] <- Ie[t-1,g]*(1-exp(-sigma*(SR[g])))
                
            pHD[t,g] <- H[t-1,g]*(1-exp(-omega*IFR[g]))
            pHR[t,g] <- H[t-1,g]*(1-exp(-omega*(1-IFR[g])))

                
            #### Vaccination
                
            if (t < week_inter){
            nVAC <- (diff_vac_timet[1,g]*popsize[g])/7

            if (nVAC >= 0){pvSS[t,g] <- nVAC
                           pvSS_rev[t,g] <- 0} else
            {nVAC_rev <- -nVAC
             pvSS[t,g] <- 0
             pvSS_rev[t,g] <- nVAC_rev}
            }
                
            else{
             pvSS[t,g] <- 0
             pvSS_rev[t,g] <- 0}
                
                
            #### Dynamics in the model   
                    
            S[t,g] <- S[t-1,g] - pSEw[t,g] - pSEa[t,g] - pSEd[t,g] - pSEe[t,g] - pvSS[t,g] + pvSS_rev[t,g]
            Ew[t,g] <- Ew[t-1,g] + pSEw[t,g] - pEIw[t,g]
            Ea[t,g] <- Ea[t-1,g] + pSEa[t,g] - pEIa[t,g]
            Ed[t,g] <- Ed[t-1,g] + pSEd[t,g] - pEId[t,g]
            Ee[t,g] <- Ee[t-1,g] + pSEe[t,g] - pEIe[t,g]
            Iw[t,g] <- Iw[t-1,g] + pEIw[t,g] - pIHw[t,g] - pIRw[t,g] 
            Ia[t,g] <- Ia[t-1,g] + pEIa[t,g] - pIHa[t,g] - pIRa[t,g] 
            Id[t,g] <- Id[t-1,g] + pEId[t,g] - pIHd[t,g] - pIRd[t,g] 
            Ie[t,g] <- Ie[t-1,g] + pEIe[t,g] - pIHe[t,g] - pIRe[t,g] 
            H[t,g] <- H[t-1,g] + pIHw[t,g] + pIHa[t,g] + pIHd[t,g] + pIHe[t,g]- pHD[t,g] - pHR[t,g] 
            D[t,g] <- D[t-1,g] + pHD[t,g] 
            R[t,g] <- R[t-1,g] + pHR[t,g] + pIRw[t,g] + pIRa[t,g] + pIRd[t,g] + pIRe[t,g]

            Sv[t,g] <- Sv[t-1,g] + pvSS[t,g] - pvSS_rev[t,g]
                
            Iw[,g][Iw[,g]<0] <- 0; Ia[,g][Ia[,g]<0] <- 0; Id[,g][Id[,g]<0] <- 0; Ie[,g][Ie[,g]<0] <- 0 
                
            eff_N[t,g] <- S[t,g]+Ew[t,g]+Iw[t,g]+Ea[t,g]+Ia[t,g]+Ed[t,g]+Id[t,g]+Ee[t,g]+Ie[t,g]+R[t,g]+Sv[t,g]
               
                
            #### the nubmer of reporting cases with the convolution
            Report <- data.frame(matrix(0, ncol=ncol(init_value), nrow=t))  
            for (tau in 1:(t-1)){Report[tau,g] = (pSEw[t-tau+1,g])*infectiontoreport(tau)[tau]}   
            sum(Report[,g]) -> Reported_modelled_w[t,g]
                
            Report <- data.frame(matrix(0, ncol=ncol(init_value), nrow=t))  
            for (tau in 1:(t-1)){Report[tau,g] = (pSEa[t-tau+1,g])*infectiontoreport(tau)[tau]}   
            sum(Report[,g]) -> Reported_modelled_a[t,g]
                
            Report <- data.frame(matrix(0, ncol=ncol(init_value), nrow=t))  
            for (tau in 1:(t-1)){Report[tau,g] = (pSEd[t-tau+1,g])*infectiontoreport(tau)[tau]}   
            sum(Report[,g]) -> Reported_modelled_d[t,g]
                
            Report <- data.frame(matrix(0, ncol=ncol(init_value), nrow=t))  
            for (tau in 1:(t-1)){Report[tau,g] = (pSEe[t-tau+1,g])*infectiontoreport(tau)[tau]}   
            sum(Report[,g]) -> Reported_modelled_e[t,g]
            
            Reported_modelled_w[t,g]+Reported_modelled_a[t,g]+
            Reported_modelled_d[t,g]+Reported_modelled_e[t,g] -> Reported_modelled[t,g]             
            }
            
        }
        
        #### simulation results
        Reported_modelled -> vac_case_final[[p]]
        H -> vac_severe_final[[p]]
        D -> vac_deaths_final[[p]]
        
    }
                    
    list(vac_case_final, vac_severe_final, vac_deaths_final) -> simul_results_list[[k]]
}<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2001306823485485, 'fraction_numerical': 0.036637729860916644, 'mean_word_length': 2.7168382097837402, 'pattern_counts': {'":': 0, '<': 407, '<?xml version=': 0, '>': 660, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 48, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '479273', 'license': 'MIT', 'n_tokens_mistral': 27975, 'n_tokens_neox': 25514, 'n_words': 3397, 'path': 'Open_code_Korea_projection.ipynb', 'repo_name': 'SungmokJung/Projection_COVID_Korea', 'size': '64278'}