In kylieainslie/vacamole: Deterministic Age-Structured SEIR model with Vaccination

knitr::opts_chunk$set(echo = TRUE)

Overview

This vignette provides detailed code to reproduce the main analysis the paper. Briefly, we quantify the benefits of extending COVID-19 vaccination beyond adults by simulating the model (described in the Quick Start Guide) forward in time under different vaccination strategies. Specifically, we compare daily cases, hospital admissions, and intensive care (IC) admissions for vaccination in adults only, those 12 years and above, and those 5 years and above. All of the functions used here are documented, so we encourage users to refer to the helps files.

This vignette demonstrates only how to reproduce the forward simulations. Code to reproduce figures from the main text in the accompanying paper can be found in the inst/extdata/scripts/manuscript folder of the package.

Installation and requirements

vacamole may be installed from github using the devtools package. There are a number of additional packages that we need for this analysis.

# Required to run vacamole 
devtools::install_github("kylieainslie/vacamole")
library(vacamole)
## Required for this analysis
library(deSolve)
library(reshape2)
library(ggplot2)
library(dplyr)
library(stringr)
library(tidyr)
library(readxl)
library(rARPACK)
library(readr)
library(lubridate)
library(kableExtra)

Preparing the data and model inputs

Vaccination schedules

A synthetic vaccination s First we need to load the vaccination schedules for the different vaccination strategies: adults only (18+), extension to adolescents (12+), and extension to adolescents and children (5+). These files are located in the inst/extdata/inputs folder of the package. Each file includes a column for vaccine product and dose for each age group (see Quick Start Guide). Rows indicate the proportion of that age group who received that dose and product on that day.

basis_5plus <- read_csv("../inst/extdata/inputs/vac_schedule_5plus.csv") %>%
  select(-starts_with("...1"))

basis_12plus <- read_csv("../inst/extdata/inputs/vac_schedule_12plus.csv") %>%
  select(-starts_with("...1"))

basis_18plus <- read_csv("../inst/extdata/inputs/vac_schedule_18plus.csv") %>%
  select(-starts_with("...1"))

Next, the vaccination schedules needs to be converted before input into the model. The vaccination schedule data frame is converted into a list with weighted averages of the following quantities: vaccination rate, vaccine effectiveness, and delay to protection based on the number of people in the population who receive each dose and type of vaccine at each time point. The model is designed to incorporate a single vaccine product with a 2-dose regimen that 1) reduces susceptibility to infection, 2) reduces risk of hospitalization if a vaccinated individual is infected, and 3) reduces risk of infecting others (transmission) if a vaccinated person is infected. Therefore we need to provide estimates for vaccine effectiveness against infection, hospitalisation, and transmission (if infected). Since there is more than one vaccine product currently licensed for use against COVID-19 in The Netherlands (the vaccines made by Pfizer/BioNTech, Moderna, AstraZeneca, and Janssen), we incorporate different vaccine products by taking the daily weighted average of the number of people with each vaccine product (and dose), the corresponding delay to protection of each vaccine product, and the vaccine effectiveness against each outcome (Table S1). Janssen is incorporated by using zero for the number of second doses at all time points.

# read in vaccine effectiveness parameters -------------------
# these will be used when converting the vaccine schedule into 
# weighted vaccine effectiveness, delay to protection, and 
# vaccination rate
ve_info <- readRDS("../inst/extdata/inputs/ve_params.rds")

# convert vaccination schedules ------------------------------
# the vaccination schedule is assumed to be cumulative over 
# time

# 5+ --------------------------------------------------------
converted_5plus <- convert_vac_schedule(
  vac_schedule = basis_5plus,
  ve = ve_info[[1]],
  delay = ve_info[[2]],
  hosp_multiplier = ve_info[[3]],
  ve_trans = ve_info[[4]], 
  add_child_vac = TRUE
)

# 12+ -------------------------------------------------------
converted_12plus <- convert_vac_schedule(
  vac_schedule = basis_12plus,
  ve = ve_info[[1]],
  delay = ve_info[[2]],
  hosp_multiplier = ve_info[[3]],
  ve_trans = ve_info[[4]]
)

# 18+ -------------------------------------------------------
converted_18plus <- convert_vac_schedule(
  vac_schedule = basis_18plus,
  ve = ve_info[[1]],
  delay = ve_info[[2]],
  hosp_multiplier = ve_info[[3]],
  ve_trans = ve_info[[4]]
)

Parameter values

The model expects a named list of parameter values (Table 1) describing the rates of transitioning between compartments, the converted vaccine schedule, and contact patterns. The model allows the user to specify a threshold of either cases (use_cases = TRUE) or IC admissions to determine which contact patterns to use. There are four different patterns the user can specify: normal (pre-pandemic), very relaxed, relaxed, and lockdown. These are specified using the c_normal, c_very_relaxed, c_relaxed, and c_lockdown arguments. The thresholds for changing the contact patterns are specified using thresh_n, thresh_l, thresh_m, and thresh_u arguments. If the user does not want the contact patterns to change, they can specify keep_cm_fixed = TRUE and specify the contact patterns to use with the c_start argument. Many of the parameter values used for this analysis are located in the inst/extdata/inputs folder of the package.

params_table <- read_xlsx("parameter_description_table.xlsx") %>%
  rename(Parameter = parameter, Type = type, Description = description) 

params_table %>%
  kbl() %>%
  kable_styling()

# read in transition rates -----------------------------------
transition_rates <- readRDS("../inst/extdata/inputs/transition_rates.rds")

# read in transmission matrices ------------------------------
april_2017 <- readRDS("../inst/extdata/inputs/contact_matrices/converted/transmission_matrix_april_2017.rds")
april_2020 <- readRDS("../inst/extdata/inputs/contact_matrices/converted/transmission_matrix_april_2020.rds")
june_2020 <- readRDS("../inst/extdata/inputs/contact_matrices/converted/transmission_matrix_june_2020.rds")
september_2020 <- readRDS("../inst/extdata/inputs/contact_matrices/converted/transmission_matrix_september_2020.rds")
february_2021 <- readRDS("../inst/extdata/inputs/contact_matrices/converted/transmission_matrix_february_2021.rds")
june_2021 <- readRDS("../inst/extdata/inputs/contact_matrices/converted/transmission_matrix_june_2021.rds")

# create list of transmission matrices -----------------------
cm <- list(april_2017 = april_2017,
           april_2020 = april_2020,
           june_2020 = june_2020,
           september_2020 = september_2020,
           february_2021 = february_2021,
           june_2021 = june_2021)

# define population size (by age group) ----------------------
age_dist <- c(0.10319920, 0.11620856, 0.12740219, 0.12198707, 0.13083463,0.14514332, 0.12092904, 0.08807406, 0.04622194) # source: CBS
n <- 17407585 # Dutch population size
n_vec <- n * age_dist

# create named list of parameter -----------------------------
params <- list(beta = 0.0004 ,  # transmission rate
               beta1 = 0.14,    # amplitude of seasonal forcing
               gamma = 0.5,     # 1/gamma = infectious period
               sigma = 0.5,     # 1/sigma = latent period
               epsilon = 0.01,  # case importation rate
               N = n_vec,       # population size vector (by age group)
               h = transition_rates$h,           # rate from infection to hospital admission
               i1 = transition_rates$i1,         # rate from hospital to IC
               i2 = transition_rates$i2,         # rate from IC back to hospital
               d = transition_rates$d,           # rate from hospital to death
               d_ic = transition_rates$d_ic,     # rate from IC to death
               d_hic = transition_rates$d_hic,   # rate from hospital (after returning from IC) to death
               r = transition_rates$r,           # rate of recovery from hospital
               r_ic = transition_rates$r_ic,     # rate of recovery from hospital (after IC)
               p_report = 1/3,  # case ascertainment rate
               c_start = june_2021,
               c_lockdown = february_2021,
               c_relaxed = june_2020,
               c_very_relaxed = june_2021,
               c_normal = april_2017,
               keep_cm_fixed = FALSE, # if true c_start is used for entire simulation, else not
               vac_inputs = basis_18plus, # converted vaccine schedule - rates per day for doses 1,2 and VE for each day (infection, transmission, hospitalization), delay to protection
               use_cases = TRUE, # if true threshold is for cases, if false threshold is for ic 
               thresh_n = 0.5/100000 * sum(n_vec), # number of cases per 100000: threshold for normal contact matrix 
               thresh_l = 5/100000 * sum(n_vec), # government threshold           
               thresh_m = 14.3/100000 * sum(n_vec), # government threshold
               thresh_u = 100000/100000 * sum(n_vec), # no lockdown unless everyone infected
               no_vac = FALSE,  # if TRUE, no vaccination is assumed
               t_calendar_start = yday(as.Date("2020-01-01")), # calendar start date (ex: if model starts on 31 Jan, then t_calendar_start = 31)
               beta_change = NULL # value to change beta to after normal contact patterns are initiated
)

Initial conditions

The model is partitioned into nine age groups, so users need to supply the initial number of individuals in each age group for each compartment. The initial conditions should be a named list. For the forward simulations we use the number of individuals in each compartment on the last day of the model fit as the initial conditions. Here, we read in previously performed model fits. A separate vignette details the model fit procedure.

# specify file path where model fit outputs are stored --------
file_path <- "../inst/extdata/results/model_fits/"

date_of_fit <- "2021-10-01" # this specifies the suffix on the model fit file

output_from_model_fit <- readRDS(paste0(file_path,"output_from_fits_", date_of_fit, ".rds"))

# get the number of individuals in each state on the last day
# of model fit (22 June 2021) ---------------------------------
# these will be used as the initial conditions
init_cond_22june2021 <- unlist(lapply(unname(output_from_model_fit$`end_date_2021-06-22`), tail, 1))

For the beginning of the forward simulations we use the transmission rate from the last time window of the model fit. However, once non-pharmaceutical interventions (NPIs) are relaxed we need a transmission rate that was not estimated during a time period in which NPIs were in use. Thus, for forward simulations, we use the transmission rate that corresponds to the basic reproduction number of the Delta variant ($R_0 = 4.6$). To incorporate uncertainty, we draw random samples of the reproduction number from a normal distribution with mean 4.6 and standard deviation equal to the square root of the variance from the model fits of the last time window before the forward simulations.

# read in MLE of transmission rate from model fit -------------
beta_mles <- data.frame(beta = readRDS(paste0(file_path, "mles_from_fits_", date_of_fit, ".rds"))) %>%
  mutate(end_date = names(output_from_model_fit))
beta_draws <- readRDS(paste0(file_path, "beta_draws_from_fits_", date_of_fit, ".rds"))

index <- which(beta_mles$end_date == "end_date_2021-06-22")

# determine Rt and beta for forward simulations ---------------
my_sd <- sqrt(0.0000945)  # sqrt(var) of Rt from fits to last time window
S_diag <- diag(init_cond_22june2021[c(2:10)]) # diagonal matrix of susceptibles
rho <- as.numeric(eigs(S_diag %*% params$c_normal, 1)$values) # parameter to convert Rt to beta

# R draws -----------------------------------------------------
rt_for_delta_period <- rnorm(n = 200, mean = 4.6, sd = my_sd)  # R_0=2.3 => beta = 0.0003934816 * 2 = 0.0007869632 (follows from different analysis)

# convert to beta ---------------------------------------------
betas <- c(0.0007869632, (rt_for_delta_period / rho) * params$gamma)

Run model

Now that we have prepared all of the model inputs, we can now run the model for each vaccination scenario for different values of the transmission rate (beta) and the contact patterns using forward_sim_wrap_func(). Results are saved to the current working directory with the file name specified by the tag argument.

# -------------------------------------------------------------
# run models --------------------------------------------------
# -------------------------------------------------------------
todays_date <- Sys.Date()

# 5+ ----------------------------------------------------------
forward_sim_func_wrap(params = params,
                      start_date = "2021-06-22",
                      end_date = "2022-03-31",
                      init_cond = init_cond_22june2021,
                      beta_m = beta_mles[index, 1],
                      vac_inputs = converted_5plus,
                      beta_c = betas,
                      beta_d = beta_draws[[index]][, 1],
                      t_normal = yday(as.Date("2021-11-01")) + 365, # simulation is in days since start
                      contact_matrices = cm,
                      tag = paste0("results_5plus_", todays_date))

# 12+ --------------------------------------------------------
forward_sim_func_wrap(params = params,
                      start_date = "2021-06-22",
                      end_date = "2022-03-31",
                      init_cond = init_cond_22june2021,
                      beta_m = beta_mles[index, 1],
                      vac_inputs = converted_12plus,
                      beta_c = betas,
                      beta_d = beta_draws[[index]][, 1],
                      t_normal = yday(as.Date("2021-11-01")) + 365,
                      contact_matrices = cm,
                      tag = paste0("results_12plus_", todays_date))

# 18+ -------------------------------------------------------
forward_sim_func_wrap(params = params,
                      start_date = "2021-06-22",
                      end_date = "2022-03-31",
                      init_cond = init_cond_22june2021,
                      beta_m = beta_mles[index,1],
                      vac_inputs = converted_18plus,
                      beta_c = betas,
                      beta_d = beta_draws[[index]][,1],
                      t_normal = yday(as.Date("2021-11-01")) + 365,
                      contact_matrices = cm,
                      tag = paste0("results_18plus_",todays_date))

Plot results

The output from the model runs needs quite a bit of data wrangling to get it into a form that's easy to plot. For simplicity, the code to perform the data wrangling is located in inst/extdata/scripts/manuscript. Here, we will just load the already wrangled results and show the code to make figure 1 from the paper. The scripts to make the other figures are located in inst/extdata/scripts/manuscript. ```r

load wrangled results -------------------------------------

all_res_for_plot <- readRDS("../inst/extdata/results/all_res_for_plot.rds")

Figure 1 (main) -------------------------------------------

fig1a <- ggplot(data = all_res_for_plot %>% filter(outcome == "Daily Cases", date >= as.Date("2021-11-01"), Immunity == "No Waning") %>% group_by(Scenario, age_group2, date, outcome) %>% summarise_at(.vars = c("mle", "lower", "upper"), .funs = "sum"), aes(x = date, y = mle, fill = age_group2,linetype = Scenario)) + geom_ribbon(aes(ymin = lower, ymax = upper, fill = age_group2), alpha = 0.3) + geom_line(aes(color = age_group2), size = 1) + scale_linetype_manual(values = c("dotted", "dashed", "solid")) + labs(y = "Cases per day", x = "Date of infection") + ylim(0,NA) + scale_x_date(date_breaks = "2 weeks", date_labels = "%d %b %Y") + theme(legend.position = "bottom", panel.background = element_blank(), axis.text.x = element_text(angle = 45, hjust = 1, size = 14), axis.text.y = element_text(size = 14), strip.text.x = element_text(size = 14), legend.text = element_text(size = 14), legend.title = element_text(size = 14), legend.box="vertical", axis.title=element_text(size=14,face="bold")) + guides(fill=guide_legend("Age Group"), colour = guide_legend("Age Group"), linetype = guide_legend("Strategy"))

Figure 1 (inset) ------------------------------------------

bar plot of percent difference ----------------------------

Table 1 ---------------------------------------------------

table1 <- all_res_for_plot %>% filter(outcome != "Daily Deaths", Immunity == "No Waning") %>% group_by(Scenario, age_group2, outcome) %>% summarise_at(.vars = c("mle", "lower", "upper"), .funs = "sum")

calculate percent difference

table1a <- table1 %>% group_by(age_group2, outcome) %>% mutate(abs_diff = mle - mle[Scenario == "Vaccination of 18+"], abs_diff_lower = lower - lower[Scenario == "Vaccination of 18+"], abs_diff_upper = upper - upper[Scenario == "Vaccination of 18+"], perc_diff = (mle * 100)/mle[Scenario == "Vaccination of 18+"] - 100, perc_diff_lower = (lower * 100)/lower[Scenario == "Vaccination of 18+"] - 100, perc_diff_upper = (upper * 100)/upper[Scenario == "Vaccination of 18+"] - 100) %>% mutate_if(is.numeric, round, 1) %>% as.data.frame()

figure 1 inset bar chart ---------------------------------

fig1_inset <- ggplot(data = table1a %>% filter(Scenario != "Vaccination of 18+"), aes(x = outcome, y = abs(perc_diff), fill = age_group2)) + geom_bar(stat = "Identity", position = position_dodge()) + geom_errorbar(aes(ymin = abs(perc_diff_upper), ymax = abs(perc_diff_lower), width = 0.2), position = position_dodge(0.9)) + labs(x = "", y = "Percent Reduction (%)", fill = "Age Group") + scale_x_discrete(labels = c("Daily\nCases", "Hospital\nAdmissions", "IC\n Admissions")) + facet_wrap(~Scenario, nrow = 2) + theme(legend.position = "bottom", panel.background = element_blank(), axis.text.x = element_text(size = 12), #angle = 45, hjust = 1, axis.text.y = element_text(size = 14), strip.text.x = element_text(size = 14), legend.text = element_text(size = 14), legend.title = element_text(size = 14), axis.title=element_text(size=14,face = "bold")) + guides(fill=guide_legend(nrow=2,byrow=TRUE))

fig1 <- fig1a + annotation_custom(ggplotGrob(fig1_inset), xmin = as.Date("2022-02-01"), xmax = as.Date("2022-04-04"), ymin = 15000, ymax = 100000) fig1

kylieainslie/vacamole documentation built on Oct. 15, 2024, 10:17 a.m.