Add test for stats correctness of high-rate vaccination; update vax v…

…ignette

Rm unused libs

tests/testthat/test-model_default.R

@@ -243,6 +243,22 @@ test_that("Default model: vaccination and stats. correctness", {
  expect_true(
  all(epidemic_size(data_baseline) > epidemic_size(data))
  )
+
+ # expect that high vaccination rates (± 1% per day, e.g. Covid-19 vax)
+ # give statistically correct results (no negative susceptibles at any time)
+ high_rate_vax <- vaccination(
+ time_begin = matrix(200, nrow(contact_matrix)),
+ time_end = matrix(200 + 150, nrow(contact_matrix)),
+ nu = matrix(0.01, nrow = nrow(contact_matrix))
+ )
+ data <- model_default(
+ uk_population,
+ vaccination = high_rate_vax, time_end = 600
+ )
+ checkmate::expect_numeric(
+ data[grepl("susceptible", data$compartment, fixed = TRUE), ]$value,
+ lower = 0
+ )
 })
 
 test_that("Default model: time dependence", {

tests/testthat/test-model_vacamole.R

@@ -343,6 +343,31 @@ test_that("Vacamole model: two dose vaccination and stats. correctness", {
  data[grepl("dose", data$compartment, fixed = TRUE) &
  time %in% seq(50, 55), ]
  )
+
+ # expect that high vaccination rates (± 1% per day, e.g. Covid-19 vax)
+ # give statistically correct results (no negative susceptibles at any time)
+ high_rate_vax <- vaccination(
+ nu = matrix(
+ c(1e-2, 1e-2), # 1% per day
+ nrow = 2, ncol = 2, byrow = TRUE
+ ),
+ time_begin = matrix(
+ c(0, 50),
+ nrow = 2, ncol = 2, byrow = TRUE # second dose given from t = 50 onwards
+ ),
+ time_end = matrix(
+ 100,
+ nrow = 2, ncol = 2
+ )
+ )
+ data <- model_vacamole(
+ uk_population,
+ vaccination = high_rate_vax, time_end = 600
+ )
+ checkmate::expect_numeric(
+ data[grepl("susceptible", data$compartment, fixed = TRUE), ]$value,
+ lower = 0
+ )
 })
 
 test_that("Vacamole model: time dependence", {

vignettes/modelling_vaccination.Rmd

@@ -22,8 +22,6 @@ knitr::opts_chunk$set(
 
 ```{r setup}
 library(epidemics)
-library(dplyr)
-library(ggplot2)
 ```
 
 ## Prepare population and initial conditions
@@ -96,7 +94,15 @@ vaccinate_elders <- vaccination(
 vaccinate_elders
 ```
 
-## Run epidemic model
+::: {.alert .alert-warning}
+**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.
 
 ```{r}
 # run an epidemic model using `epidemic`
@@ -107,42 +113,6 @@ output <- model_default(
 )
 ```
 
-## Prepare data and visualise infections
-
-Plot epidemic over time, showing only the number of individuals in the exposed, infected, and vaccinated compartments.
-
-```{r class.source = 'fold-hide'}
-# 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"
- )
-```
+::: {.alert .alert-info}
+**Note that** vaccination modelling is currently only supported for the 'default' (`model_default()`) and 'Vacamole' (`model_vacamole()`) models.
+:::