epiverse-trace / epidemics

A library of published compartmental epidemic models, and classes to represent demographic structure, non-pharmaceutical interventions, and vaccination regimes, to compose epidemic scenarios.
https://epiverse-trace.github.io/epidemics/
Other
8 stars 2 forks source link

Add Vacamole model to {epidemics} #54

Closed pratikunterwegs closed 1 year ago

pratikunterwegs commented 1 year ago

This PR aims to add an initial version of the Vacamole model to {epidemics}. For more information on the Vacamole model and its adaptation for {epidemics}, see the Discussion issue #44.

A good place to start with reviewing this PR is to read the vignette titled vacamole.Rmd; this is written with {pkgdown} presentation in mind.

  1. Fixes #50 by adding the Vacamole model function details to the model library JSON,
  2. Fixes #51 by adding a vignette on using Vacamole for modelling the effectiveness of vaccination in reducing deaths and hospitalisations,
  3. Fixes #52 by adding the Vacamole model in a header file,
  4. Fixes #53 and fixes #55 by updating vaccination() to allow multi-dose vaccination, as well as implementing a c() method for vaccination objects that allows constructing multi-dose vaccinations by concatenating multiple vaccination objects,
  5. Fixes #56 by editing epidemic_size() to count deaths in epidemic size by default if a "dead" compartment is found.
  6. Adds safer testing of the add_to_library() function using {withr}, taking on this dependency in "Suggests" - this needs review,
  7. Adds spellchecking after rebasing on main following merging of PR #60, and updates the WORDLIST.

Please see PR #57 that contains changes that will be added on top of this PR - these importantly include updates to helper functions that are relevent to the changes in this PR.

pratikunterwegs commented 1 year ago

Add benchmarking output here; approx mean of 23 milliseconds per run. If runtime scales linearly, 200 scenarios with 1000 runs each would take 1h15 mins.

library(epidemics)

polymod <- socialmixr::polymod
contact_data <- socialmixr::contact_matrix(
  polymod,
  countries = "United Kingdom",
  age.limits = c(0, 60),
  symmetric = TRUE
)
#> Using POLYMOD social contact data. To cite this in a publication, use the 'cite' function
#> Removing participants that have contacts without age information. To change this behaviour, set the 'missing.contact.age' option
contact_matrix <- t(contact_data$matrix)
demography_vector <- contact_data$demography$population

# // 0| 1| 2|3| 4|5| 6|7| 8|9|10
# // S|V1|V2|E|EV|I|IV|H|HV|D|R

# make initial conditions - order is important
initial_conditions <- c(
  S = 1 - 1e-6,
  V1 = 0, V2 = 0,
  E = 0, EV = 0,
  I = 1e-6, IV = 0,
  H = 0, HV = 0, D = 0, R = 0
)
initial_conditions <- rbind(
  initial_conditions,
  initial_conditions
)

# create a population
uk_population <- population(
  name = "UK population",
  contact_matrix = matrix(1, 2, 2),
  demography_vector = 67e6 * c(0.4, 0.6),
  initial_conditions = initial_conditions
)

# prepare a two dose vaccination regime for three age groups
double_vaccination <- vaccination(
  name = "double_vaccination",
  nu = matrix(
    1e-3,
    nrow = 3, ncol = 2
  ),
  time_begin = matrix(
    00,
    nrow = 3, ncol = 2
  ),
  time_end = matrix(
    200,
    nrow = 3, ncol = 2
  )
)

# make infection class for Vacamole model
# note extra arguments passed as ...
infect <- infection(
  name = "covid", r0 = 5, infectious_period = 10,
  preinfectious_period = 5,
  eta = 1 / 1000, omega = 1 / 1000,
  susc_reduction_vax = 0.5,
  hosp_reduction_vax = 0.7,
  mort_reduction_vax = 0.9
)

microbenchmark::microbenchmark(
  "Vacacmole" = epidemic(
    model_name = "vacamole",
    population = uk_population,
    infection = infect,
    vaccination = double_vaccination,
    time_end = 400, increment = 1
  ),
  times = 1000
)
#> Warning in microbenchmark::microbenchmark(Vacacmole = epidemic(model_name
#> = "vacamole", : less accurate nanosecond times to avoid potential integer
#> overflows
#> Unit: milliseconds
#>       expr      min       lq   mean   median       uq      max neval
#>  Vacacmole 15.17328 22.15613 22.888 23.14522 23.74511 30.98645  1000

Created on 2023-05-24 by the reprex package (v2.0.1)

pratikunterwegs commented 1 year ago

Thanks @Bisaloo and @bahadzie for your feedback - I have created issues #58 and #59 to discuss changes to the infection and vaccination classes. @Bisaloo could you please see if the open conversations can be resolved for now, or if you would prefer another solution. Will also wait for advice on inviting original Vacamole devs to review.