Allow time-varying, age-specific vaccination rates

[pratikunterwegs]

This issue comes from discussions with @sbfnk re: supporting the use of Vacamole for submissions to the ECDC Scenario Hub (by collaborating users).

The aim is to provide more functionality around specifying age-stratified time-varying vaccination rates which would allow <vaccination> objects to more accurately represent real vaccine calendars.

The main open question around the actual implementation is how to represent vaccine calendars in the Vacamole model which has two vaccination compartments (applies to any number of compartments).

[pratikunterwegs]

Continuing this discussion with a small reprex (most of it's setup!) on what a time-varying vaccination functionality could look like.

The <vaccination> class has been refactored to have two elements:

1. A $T \times 2$ matrix of start and end times for each of $T$ intervals such as weeks or phases, and
2. A $T \times N$ column matrix of vaccination rates for each demographic group in each vaccination interval (e.g. doses administered to 65+ in week 1).

Some points for consideration, happy to hear any suggestions:

1. This new <vaccination> implementation needs work before it can be used with the Vacamole model as I need to find a good way of representing the calendar for second (or any N) doses in a model. 3D arrays don't appear to be a good way of doing this as they are more challenging to pass to C++ from R.
2. Here, the class member nu breaks from the usual {epidemics} convention and has vaccination rates for each group in columns rather than rows - should this be standardised so group-specific rates are in rows?
3. Should the Vacamole model be refactored to be a single dose model? This would be easier than refactoring <vaccination> to satisfy point (1) above. The argument in favour is that the $V_1$ compartment does not protect against infection or severe outcomes at the moment. If this compartment is retained, should it also confer some protection?

Reprex for time-varying vaccination

library(ggplot2)
library(epidemics)

# Prepare contact matrix and demography vector
polymod <- socialmixr::polymod
contact_data <- socialmixr::contact_matrix(
  polymod,
  countries = "United Kingdom",
  age.limits = c(0, 60),
  symmetric = TRUE
)

contact_matrix <- t(contact_data$matrix)
demography_vector <- contact_data$demography$population

# make initial conditions - order is important
initial_conditions <- c(
  S = 1 - 1e-6, E = 0,
  I = 1e-6, R = 0, V = 0
)
initial_conditions <- rbind(
  initial_conditions,
  initial_conditions
)

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

# get intervals for vaccination stages/phases/weeks
start_times <- seq(0, 91, 7)
end_times <- start_times + 6

intervals <- cbind(start_times, end_times)

# get dosage for a single age group - extend to multiple
dosage <- runif(nrow(intervals), 1 / 1e3, 4 / 1e3)
dosage <- cbind(dosage, dosage)

vax <- vaccination(
  nu = dosage,
  intervals = intervals
)

vax
#> <vaccination> object
#> Intervals:
#>       start_times end_times
#>  [1,]           0         6
#>  [2,]           7        13
#>  [3,]          14        20
#>  [4,]          21        27
#>  [5,]          28        34
#>  [6,]          35        41
#>  [7,]          42        48
#>  [8,]          49        55
#>  [9,]          56        62
#> [10,]          63        69
#> [11,]          70        76
#> [12,]          77        83
#> [13,]          84        90
#> [14,]          91        97
#> Rates:
#>            dosage      dosage
#>  [1,] 0.003198568 0.003198568
#>  [2,] 0.003311817 0.003311817
#>  [3,] 0.003494860 0.003494860
#>  [4,] 0.002264464 0.002264464
#>  [5,] 0.002723691 0.002723691
#>  [6,] 0.003272980 0.003272980
#>  [7,] 0.001669557 0.001669557
#>  [8,] 0.001807109 0.001807109
#>  [9,] 0.002605438 0.002605438
#> [10,] 0.003692647 0.003692647
#> [11,] 0.002444341 0.002444341
#> [12,] 0.003510983 0.003510983
#> [13,] 0.001281554 0.001281554
#> [14,] 0.002998309 0.002998309

# numbers vaccinated each day in each week
vax$nu %*% diag(demography_vector)
#>            [,1]     [,2]
#>  [1,] 151891.18 40940.11
#>  [2,] 157269.04 42389.63
#>  [3,] 165961.25 44732.49
#>  [4,] 107533.16 28984.03
#>  [5,] 129340.59 34861.92
#>  [6,] 155424.77 41892.53
#>  [7,]  79282.67 21369.52
#>  [8,]  85814.61 23130.11
#>  [9,] 123725.08 33348.33
#> [10,] 175353.64 47264.08
#> [11,] 116074.99 31286.36
#> [12,] 166726.92 44938.87
#> [13,]  60857.49 16403.27
#> [14,] 142381.42 38376.88

# model run time
time_end <- 100L

data <- model_default(
  population = uk_population,
  transmission_rate = 1.3 / 7,
  vaccination = vax,
  time_end = time_end
)

ggplot(data[compartment == "vaccinated", ]) +
  geom_line(
    aes(time, value, col = demography_group)
  ) +
  labs(y = "# Vaccinated")

^{Created on 2024-05-08 with reprex v2.0.2}