Modelling the effect of a vaccination campaign

Code

library(epidemics)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(ggplot2)

Prepare population and initial conditions

Prepare population and contact data.

Code

# 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
)
#> Removing participants that have contacts without age information. To change this behaviour, set the 'missing.contact.age' option

# 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.

Code

# 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.

Code

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.

Code

# 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

Run epidemic model

Code

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

Prepare data and visualise infections

Plot epidemic over time, showing only the number of individuals in the exposed, infected, and vaccinated compartments.

Code

# plot figure of epidemic curve
filter(output, compartment %in% c("exposed", "infectious")) %>%
  ggplot(
    aes(
      x = time,
      y = value,
      col = demography_group,
      linetype = compartment
    )
  ) +
  geom_line() +
  scale_y_continuous(
    labels = scales::comma
  ) +
  scale_colour_brewer(
    palette = "Dark2",
    name = "Age group"
  ) +
  expand_limits(
    y = c(0, 500e3)
  ) +
  coord_cartesian(
    expand = FALSE
  ) +
  theme_bw() +
  theme(
    legend.position = "top"
  ) +
  labs(
    x = "Simulation time (days)",
    linetype = "Compartment",
    y = "Individuals"
  )