Correct vaccination implementation

R/check_args_default.R

@@ -74,6 +74,10 @@
   vax_time_end <- mod_args[["vaccination"]][["time_end"]]
   vax_nu <- mod_args[["vaccination"]][["nu"]]
 
+  # calculate vax_nu as the NUMBER of vaccinations per timestep
+  # rather than a fraction/rate. See issue #198 for more details.
+  vax_nu <- vax_nu * mod_args[["population"]][["demography_vector"]]
+
   # remove processed model arguments
   mod_args[c("population", "intervention", "vaccination")] <- NULL
   mod_args[c("param_set", "scenario")] <- NULL

R/check_args_vacamole.R

@@ -94,6 +94,10 @@
   vax_time_end <- mod_args[["vaccination"]][["time_end"]]
   vax_nu <- mod_args[["vaccination"]][["nu"]]
 
+  # calculate vax_nu as the NUMBER of vaccinations per timestep
+  # rather than a fraction/rate. See issue #198 for more details.
+  vax_nu <- vax_nu * mod_args[["population"]][["demography_vector"]]
+
   # remove processed model arguments
   mod_args[c("population", "intervention", "vaccination")] <- NULL
   mod_args[c("param_set", "scenario")] <- NULL

inst/include/model_default.h

@@ -122,16 +122,20 @@ struct epidemic_default {
         model_params_temp.at("infectiousness_rate") * x.col(1).array();
     Eigen::ArrayXd iToR =
         model_params_temp.at("recovery_rate") * x.col(2).array();
-    Eigen::ArrayXd sToV = vax_nu_current.col(0).array() * x.col(0).array();
+
+    /// NOTE: vax_nu_current holds COUNTS, not rates. See issue #198 for details
+    Eigen::ArrayXd sToV = vax_nu_current.col(0).array();
 
     // compartmental changes accounting for contacts (for dS and dE)
-    // β: transmission_rate; ν: vaccination rate; σ: infectiousness rate
+    // β: transmission_rate; ν: vaccination COUNTS, see issue #198
+    /// NOTE: see issue #198 for more on vaccination
+    // σ: infectiousness rate
     // γ: recovery rate
-    dxdt.col(0) = -sToE - sToV;  // -β*S*contacts*I - ν*S
+    dxdt.col(0) = -sToE - sToV;  // -β*S*contacts*I - ν
     dxdt.col(1) = sToE - eToI;   // β*S*contacts*I - σ*E
     dxdt.col(2) = eToI - iToR;   // σ*E - γ*I
     dxdt.col(3) = iToR;          // γ*I
-    dxdt.col(4) = sToV;          // ν*S
+    dxdt.col(4) = sToV;          // ν
   }
 };
 

inst/include/model_vacamole.h

@@ -142,10 +142,11 @@ struct epidemic_vacamole {
                           x.col(0).array() *
                           (cm_temp * (x.col(5) + x.col(6))).array();
 
+    /// NOTE: vax_nu_current holds COUNTS, not rates. See issue #198 for details
     // Susceptible to vaccinated with one dose
-    Eigen::ArrayXd sToV1 = vax_nu_current.col(0).array() * x.col(0).array();
+    Eigen::ArrayXd sToV1 = vax_nu_current.col(0).array();
     // Vaccinated one dose to vaccinated with two doses
-    Eigen::ArrayXd v1ToV2 = vax_nu_current.col(1).array() * x.col(1).array();
+    Eigen::ArrayXd v1ToV2 = vax_nu_current.col(1).array();
 
     // Vaccinated one dose to exposed - same as susceptible to exposed
     Eigen::ArrayXd v1ToE = model_params_temp.at("transmission_rate") *
@@ -200,13 +201,14 @@ struct epidemic_vacamole {
 
     // compartmental changes accounting for contacts
     // β: transmission_rate; βv: transmission_rate for doubly vaccinated;
-    // ν1, ν2: vaccination rate first, second dose;
+    // ν1, ν2: vaccination COUNTS for first, second dose;
+    /// NOTE: see issue #198 for more on vaccination
     // σ: infectiousness rate; γ: recovery rate;
     // η: hospitalisation rate; η_v: hosp. rate doubly vaccinated
     // ω: mortality rate; ω_v: mort. rate doubly vaccinated
-    dxdt.col(0) = -sToE - sToV1;           // -β*S*contacts*(I+Iv) - ν1*S
-    dxdt.col(1) = sToV1 - v1ToE - v1ToV2;  // ν1*S -β*V1*contacts*(I+Iv) - ν2*V
-    dxdt.col(2) = v1ToV2 - v2ToEv;         // ν2*V - βv*V2*contacts*(I+Iv)
+    dxdt.col(0) = -sToE - sToV1;           // -β*S*contacts*(I+Iv) - ν1
+    dxdt.col(1) = sToV1 - v1ToE - v1ToV2;  // ν1 -β*V1*contacts*(I+Iv) - ν2
+    dxdt.col(2) = v1ToV2 - v2ToEv;         // ν2 - βv*V2*contacts*(I+Iv)
     dxdt.col(3) = sToE + v1ToE - eToI;     // β*(S+V1)*contacts*(I+Iv) - σE
     dxdt.col(4) = v2ToEv - evToIv;         // βv*V2*contacts*(I+Iv) - αE
     dxdt.col(5) = eToI - iToH - iToR - iToD;        // σE - I(γ + η + ω)

inst/include/vaccination.h

@@ -14,18 +14,6 @@
 // add to namespace vaccination
 namespace vaccination {
 
-/// @brief Get the group specific vaccination rate from a `vaccination` object
-/// @param vaccination A `vaccination` object in the form of an Rcpp List
-/// @return An Eigen Array of group-specific vaccination rates, with dimensions
-/// matching those the `vaccination` class member `nu`, in which each element
-/// `i,j` gives the group- and dose-specific vaccination rate, in row `i` and
-/// column `j` respectively.
-inline Eigen::MatrixXd get_nu(const Rcpp::List &vaccination) {
-  // convert to Eigen array, then resize to match the nu matrix
-  Eigen::MatrixXd nu(Rcpp::as<Eigen::MatrixXd>(vaccination["nu"]));
-  return nu;
-}
-
 /// @brief Get the vaccination rate at a time in the vaccination regime
 /// @param t A double value giving the current time, which is compared against
 /// the vaccination start and end times.

tests/testthat/_snaps/model_vacamole.md

@@ -0,0 +1,33 @@
+# Vacamole model: two dose vaccination and stats. correctness
+
+    Code
+      data[grepl("dose", data$compartment, fixed = TRUE) & time %in% seq(50, 55), ]
+    Output
+           time demography_group         compartment       value
+          <num>           <char>              <char>       <num>
+       1:    50           [0,60) vaccinated_one_dose 2370379.149
+       2:    50              60+ vaccinated_one_dose  638905.568
+       3:    50           [0,60) vaccinated_two_dose    3957.270
+       4:    50              60+ vaccinated_two_dose    1066.626
+       5:    51           [0,60) vaccinated_one_dose 2394121.139
+       6:    51              60+ vaccinated_one_dose  645305.113
+       7:    51           [0,60) vaccinated_two_dose   27700.882
+       8:    51              60+ vaccinated_two_dose    7466.380
+       9:    52           [0,60) vaccinated_one_dose 2417863.063
+      10:    52              60+ vaccinated_one_dose  651704.650
+      11:    52           [0,60) vaccinated_two_dose   51444.480
+      12:    52              60+ vaccinated_two_dose   13866.132
+      13:    53           [0,60) vaccinated_one_dose 2441604.917
+      14:    53              60+ vaccinated_one_dose  658104.178
+      15:    53           [0,60) vaccinated_two_dose   75188.064
+      16:    53              60+ vaccinated_two_dose   20265.882
+      17:    54           [0,60) vaccinated_one_dose 2465346.699
+      18:    54              60+ vaccinated_one_dose  664503.696
+      19:    54           [0,60) vaccinated_two_dose   98931.632
+      20:    54              60+ vaccinated_two_dose   26665.630
+      21:    55           [0,60) vaccinated_one_dose 2489088.407
+      22:    55              60+ vaccinated_one_dose  670903.205
+      23:    55           [0,60) vaccinated_two_dose  122675.184
+      24:    55              60+ vaccinated_two_dose   33065.377
+           time demography_group         compartment       value
+

tests/testthat/_snaps/vaccination.md

@@ -79,3 +79,27 @@
            dose_1 dose_2
       [1,]  0.001  0.001
 
+# Epidemic model with vaccination
+
+    Code
+      data_vaccination[data_vaccination$compartment == "vaccinated" &
+        data_vaccination$time < 5, ]
+    Output
+           time demography_group compartment     value
+          <num>           <char>      <char>     <num>
+       1:     0           [0,40)  vaccinated    0.0000
+       2:     0          [40,65)  vaccinated    0.0000
+       3:     0              65+  vaccinated    0.0000
+       4:     1           [0,40)  vaccinated    0.0000
+       5:     1          [40,65)  vaccinated    0.0000
+       6:     1              65+  vaccinated  967.3058
+       7:     2           [0,40)  vaccinated    0.0000
+       8:     2          [40,65)  vaccinated    0.0000
+       9:     2              65+  vaccinated 1934.6116
+      10:     3           [0,40)  vaccinated    0.0000
+      11:     3          [40,65)  vaccinated    0.0000
+      12:     3              65+  vaccinated 2901.9174
+      13:     4           [0,40)  vaccinated    0.0000
+      14:     4          [40,65)  vaccinated    0.0000
+      15:     4              65+  vaccinated 3869.2232
+

tests/testthat/test-model_vacamole.R

@@ -38,12 +38,12 @@ uk_population <- population(
 double_vaccination <- vaccination(
   name = "double_vaccination",
   nu = matrix(
-    1e-3,
-    nrow = 2, ncol = 2
+    c(1e-3, 5e-4), # fewer second doses
+    nrow = 2, ncol = 2, byrow = TRUE
   ),
   time_begin = matrix(
-    00,
-    nrow = 2, ncol = 2
+    c(0, 50),
+    nrow = 2, ncol = 2, byrow = TRUE # second dose given from t = 50 onwards
   ),
   time_end = matrix(
     100,
@@ -336,6 +336,13 @@ test_that("Vacamole model: two dose vaccination and stats. correctness", {
   expect_true(
     all(epidemic_size(data_baseline) > epidemic_size(data))
   )
+
+  # expect snapshot for range from t = 50:55
+  # second dose begins at t = 50
+  expect_snapshot(
+    data[grepl("dose", data$compartment, fixed = TRUE) &
+      time %in% seq(50, 55), ]
+  )
 })
 
 test_that("Vacamole model: time dependence", {

tests/testthat/test-vaccination.R

@@ -132,6 +132,11 @@ data <- model_default(
 )
 
 test_that("Epidemic model with vaccination", {
+  # Expect snapshot for vaccinated compartment
+  expect_snapshot(
+    data_vaccination[data_vaccination$compartment == "vaccinated" &
+      data_vaccination$time < 5, ]
+  )
   # expect that only the last age group is vaccinated
   total_vaccinated <- data_vaccination[data_vaccination$compartment ==
     "vaccinated" & data_vaccination$time ==