Modelling the effect of a vaccination campaign • epidemics

library(epidemics)

Prepare population and initial conditions

Prepare population and contact data.

# load contact and population data from socialmixr::polymod
polymod <- socialmixr::polymod
contact_data <- socialmixr::contact_matrix(
  polymod,
  countries = "United Kingdom",
  age_limits = c(0, 20, 65),
  symmetric = TRUE,
  return_demography = TRUE
)
#> Warning: Automatic country population lookup in `contact_matrix()` was deprecated in
#> socialmixr 0.6.0.
#> When `countries` is given (or a `country` column is present) without
#> `survey_pop`, contact_matrix() currently calls the soft-deprecated `wpp_age()`
#> to look up population data. This automatic lookup will be removed in a future
#> release: callers will then have to supply `survey_pop` whenever `symmetric`,
#> `split`, `per_capita`, `weigh_age`, or `return_demography` is TRUE.
#> ℹ Pass `survey_pop` explicitly to silence this warning, e.g. `survey_pop =
#>   survey_country_population(survey, countries)` or a data frame from the
#>   wpp2024 package.
#> This warning is displayed once per session.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
#> Warning in normalise_weighted_matrix(survey_pop = survey_pop, weighted_matrix = weighted.matrix, : Large differences in the size of the sub-populations with the current age
#> breaks are likely to result in artefacts after making the matrix symmetric.
#> ! Please reconsider the age breaks to obtain more equally sized
#>   sub-populations.
#> ℹ Normalization factors: [0.5 and 2.2]

# prepare contact matrix
contact_matrix <- t(contact_data$matrix)

# prepare the demography vector
demography_vector <- contact_data$demography$population
names(demography_vector) <- rownames(contact_matrix)

Prepare initial conditions for each age group.

# initial conditions
initial_i <- 1e-6
initial_conditions <- c(
  S = 1 - initial_i, E = 0, I = initial_i, R = 0, V = 0
)

# build for all age groups
initial_conditions <- rbind(
  initial_conditions,
  initial_conditions,
  initial_conditions
)

# assign rownames for clarity
rownames(initial_conditions) <- rownames(contact_matrix)

Prepare a population as a population class object.

uk_population <- population(
  name = "UK",
  contact_matrix = contact_matrix,
  demography_vector = demography_vector,
  initial_conditions = initial_conditions
)

Prepare a vaccination campaign

Prepare a vaccination campaign targeting individuals aged over 65.

# prepare a vaccination object
vaccinate_elders <- vaccination(
  name = "vaccinate elders",
  time_begin = matrix(100, nrow(contact_matrix)),
  time_end = matrix(250, nrow(contact_matrix)),
  nu = matrix(c(0.0001, 0, 0))
)

# view vaccination object
vaccinate_elders
#> <vaccination> object
#> 
#>  Vaccination name: 
#> "vaccinate elders"
#> 
#>  Begins at: 
#>      dose_1
#> [1,]    100
#> [2,]    100
#> [3,]    100
#> 
#>  Ends at: 
#>      dose_1
#> [1,]    250
#> [2,]    250
#> [3,]    250
#> 
#>  Vaccination rate: 
#>      dose_1
#> [1,]  1e-04
#> [2,]  0e+00
#> [3,]  0e+00

Note that when the vaccination rate is high, the number of individuals who could potentially transition from ‘susceptible’ to ‘vaccinated’ may be greater than the number of individuals in the ‘susceptible’ compartment: $dS_{i_{t}} = \nu_{i_t}S_{i_{0}} > S_{i_{t}}$ for each demographic group $i$ at time $t$ .

This could realistically occur when there are more doses available than individuals eligible to receive them, for instance towards the end of an epidemic.

epidemics automatically handles such situations by setting $dS_{i_{t}}$ to be the minimum of doses or eligible individuals: : $dS_{i_t} = \text{Min}(\nu_{i_t}S_{i_0}, S_{i_t})$ , such that $S_{i_{t+1}}$ does not take a negative value.

The <vaccination> object can be passed to the vaccination argument of a model function call.

# run an epidemic model using `epidemic`
output <- model_default(
  population = uk_population,
  vaccination = vaccinate_elders,
  time_end = 600, increment = 1.0
)

Note that vaccination modelling is currently only supported for the ‘default’ (model_default()) and ‘Vacamole’ (model_vacamole()) models.