add {epidemics} summative assessments

[avallecam]

# load packages -----------------------------------------------------------

library(epidemics)
library(socialmixr)
#> 
#> Attaching package: 'socialmixr'
#> The following object is masked from 'package:utils':
#> 
#>     cite
library(tidyverse)

# parameters --------------------------------------------------------------

preinfectious_period <- 3.0
infectious_period <- 7.0
basic_reproduction <- 1.46

infectiousness_rate <- 1.0/preinfectious_period
recovery_rate <- 1.0/infectious_period
transmission_rate <- basic_reproduction/infectious_period

# paste from tutorial -----------------------------------------------------

polymod <- socialmixr::polymod
contact_data <- socialmixr::contact_matrix(
  polymod,
  countries = "Italy",
  age.limits = c(0, 20, 40),
  symmetric = TRUE
)
#> Using POLYMOD social contact data. To cite this in a publication, use the 'get_citation()' function
#> Removing participants without age information. To change this behaviour, set the 'missing.participant.age' option
#> 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)

# initial conditions: one in every 1 million is infected
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 <- matrix(
  rep(initial_conditions, dim(contact_matrix)[1]),
  ncol = 5, byrow = TRUE
)
rownames(initial_conditions) <- rownames(contact_matrix)

# prepare the population to model as affected by the epidemic
uk_population <- population(
  name = "UK",
  contact_matrix = contact_matrix,
  demography_vector = demography_vector,
  initial_conditions = initial_conditions
)

uk_population
#> <population> object
#> 
#>  Population name:
#> "UK"
#> 
#>  Demography 
#> [0,20): 11,204,261 (20%)
#> [20,40): 16,305,622 (30%)
#> 40+: 31,298,598 (50%)
#> 
#>  Contact matrix 
#>                  
#> contact.age.group    [0,20)  [20,40)      40+
#>           [0,20)  14.306502 2.255442 2.897234
#>           [20,40)  3.282356 9.126904 4.182137
#>           40+      8.093293 8.027601 8.900621

# intervention: close school ------------------------------------------------------------

rownames(contact_matrix)
#> [1] "[0,20)"  "[20,40)" "40+"

close_schools <- intervention(
  name = "School closure",
  type = "contacts",
  time_begin = 200,
  time_end = 200 + 100,
  reduction = matrix(c(0.5, 0.01, 0.01)) # group specific effect
)

close_schools
#> <contacts_intervention> object
#> 
#>  Intervention name:
#> "School closure"
#> 
#>  Begins at: 
#> [1] 200
#> 
#>  Ends at: 
#> [1] 300
#> 
#>  Reduction: 
#>              Interv. 1
#> Demo. grp. 1      0.50
#> Demo. grp. 2      0.01
#> Demo. grp. 3      0.01

# intervention: mask mandate ----------------------------------------------

mask_mandate <- intervention(
  name = "Mask mandate",
  type = "rate",
  time_begin = 200,
  time_end = 200 + 200,
  reduction = 0.163 # rate impact all infectious (reduce infectiousness)
)

mask_mandate
#> <rate_intervention> object
#> 
#>  Intervention name:
#> "Mask mandate"
#> 
#>  Begins at: 
#>      [,1]
#> [1,]  200
#> 
#>  Ends at: 
#>      [,1]
#> [1,]  400
#> 
#>  Reduction: 
#> Interv. 1 
#>     0.163


# intervention: vaccination -----------------------------------------------

vaccinate <- vaccination(
  name = "vaccinate all",
  time_begin = matrix(200, nrow(contact_matrix)),
  time_end = matrix(200 + 150, nrow(contact_matrix)),
  nu = matrix(c(0.01, 0.0, 0.01))
)


# run baseline ------------------------------------------------------------

baseline <- model_default(
  # population
  population = uk_population,
  # rates
  transmission_rate = transmission_rate,
  infectiousness_rate = infectiousness_rate,
  recovery_rate = recovery_rate,
  # time
  time_end = 600,increment = 1
)

# run intervention --------------------------------------------------------

intervention_school <- model_default(
  # population
  population = uk_population,
  # rates
  transmission_rate = transmission_rate,
  infectiousness_rate = infectiousness_rate,
  recovery_rate = recovery_rate,
  # intervention
  intervention = list(contacts = close_schools),
  # time
  time_end = 600,increment = 1
)

intervention_mask <- model_default(
  # population
  population = uk_population,
  # rates
  transmission_rate = transmission_rate,
  infectiousness_rate = infectiousness_rate,
  recovery_rate = recovery_rate,
  # intervention
  intervention = list(transmission_rate = mask_mandate),
  # time
  time_end = 600,increment = 1
)

intervention_vaccine <- model_default(
  # population
  population = uk_population,
  # rates
  transmission_rate = transmission_rate,
  infectiousness_rate = infectiousness_rate,
  recovery_rate = recovery_rate,
  # intervention
  vaccination = vaccinate,
  # time
  time_end = 600,increment = 1
)


# intervention: combine ---------------------------------------------------

intervention_mask_vaccine <- model_default(
  # population
  population = uk_population,
  # rates
  transmission_rate = transmission_rate,
  infectiousness_rate = infectiousness_rate,
  recovery_rate = recovery_rate,
  # intervention
  intervention = list(contacts = close_schools),
  # intervention = list(transmission_rate = mask_mandate),
  vaccination = vaccinate,
  # time
  time_end = 600,increment = 1
)

# # get new infections ------------------------------------------------------
# 
# baseline_infections <- new_infections(baseline,by_group = FALSE)
# intervention_infections <- new_infections(intervention,by_group = FALSE)
# 
# 
# # get finalsize -----------------------------------------------------------
# 
# epidemic_size(baseline)
# epidemic_size(intervention)

# visualize ---------------------------------------------------------------

baseline$intervention_type <- "baseline"
intervention_school$intervention_type <- "close_schools"
intervention_mask$intervention_type <- "mask_mandate"
intervention_vaccine$intervention_type <- "vaccinate"
intervention_mask_vaccine$intervention_type <- "mask + vaccine"

output <- bind_rows(
  baseline,
  intervention_school,
  intervention_mask,
  intervention_vaccine,
  intervention_mask_vaccine
)

output %>% 
  filter(compartment == "infectious") %>% 
  ggplot(
    aes(x = time, 
        y = value,
        linetype = intervention_type,
        colour = intervention_type
    )
  ) +
  stat_summary(
    fun = "sum",
    geom = "line",
    linewidth = 1
  ) +
  scale_y_continuous(labels = scales::comma) +
  geom_vline(
    xintercept = c(
      close_schools$time_begin,
      close_schools$time_end
    ),
    linetype = 2
  ) +
  geom_vline(
    xintercept = c(
      mask_mandate$time_begin,
      mask_mandate$time_end
    ),
    linetype = 3
  ) +
  geom_vline(
    xintercept = c(
      vaccinate$time_begin,
      vaccinate$time_end
    ),
    linetype = 4
  ) +
  labs(
    x = "Simulation time (days)",
    colour = "Age group",
    y = "Individuals"
  ) +
  theme_bw()

^{Created on 2024-04-26 with reprex v2.1.0}