Implement interventions on model rates

NAMESPACE

@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
-S3method(c,intervention)
+S3method(c,contacts_intervention)
+S3method(c,rate_intervention)
 S3method(c,vaccination)
 S3method(print,infection)
 S3method(print,intervention)
@@ -20,13 +21,16 @@ export(get_model_names)
 export(get_parameter)
 export(infection)
 export(intervention)
+export(is_contacts_intervention)
 export(is_infection)
 export(is_intervention)
 export(is_population)
+export(is_rate_intervention)
 export(is_vaccination)
 export(model_library)
 export(new_infections)
-export(no_intervention)
+export(no_contacts_intervention)
+export(no_rate_intervention)
 export(no_vaccination)
 export(population)
 export(vaccination)

R/RcppExports.R

@@ -41,8 +41,8 @@
 #' as specified in the initial conditions matrix (see [population()]).
 #' The second list element is a vector of timesteps.
 #' @keywords internal
-.epidemic_default_cpp <- function(initial_state, beta, alpha, gamma, contact_matrix, npi_time_begin, npi_time_end, npi_cr, vax_time_begin, vax_time_end, vax_nu, time_end = 100.0, increment = 0.1) {
-    .Call(`_epidemics_epidemic_default_cpp_internal`, initial_state, beta, alpha, gamma, contact_matrix, npi_time_begin, npi_time_end, npi_cr, vax_time_begin, vax_time_end, vax_nu, time_end, increment)
+.epidemic_default_cpp <- function(initial_state, beta, alpha, gamma, contact_matrix, npi_time_begin, npi_time_end, npi_cr, vax_time_begin, vax_time_end, vax_nu, rate_interventions, time_end = 100.0, increment = 0.1) {
+    .Call(`_epidemics_epidemic_default_cpp_internal`, initial_state, beta, alpha, gamma, contact_matrix, npi_time_begin, npi_time_end, npi_cr, vax_time_begin, vax_time_end, vax_nu, rate_interventions, time_end, increment)
 }
 
 #' @title Run an SEIR model with Erlang passage times
@@ -132,8 +132,8 @@
 #' as specified in the initial conditions matrix (see [population()]).
 #' The second list element is a vector of timesteps.
 #' @keywords internal
-.epidemic_vacamole_cpp <- function(initial_state, beta, beta_v, alpha, omega, omega_v, eta, eta_v, gamma, contact_matrix, npi_time_begin, npi_time_end, npi_cr, vax_time_begin, vax_time_end, vax_nu, time_end = 100.0, increment = 0.1) {
-    .Call(`_epidemics_epidemic_vacamole_cpp_internal`, initial_state, beta, beta_v, alpha, omega, omega_v, eta, eta_v, gamma, contact_matrix, npi_time_begin, npi_time_end, npi_cr, vax_time_begin, vax_time_end, vax_nu, time_end, increment)
+.epidemic_vacamole_cpp <- function(initial_state, beta, beta_v, alpha, omega, omega_v, eta, eta_v, gamma, contact_matrix, npi_time_begin, npi_time_end, npi_cr, vax_time_begin, vax_time_end, vax_nu, rate_interventions, time_end = 100.0, increment = 0.1) {
+    .Call(`_epidemics_epidemic_vacamole_cpp_internal`, initial_state, beta, beta_v, alpha, omega, omega_v, eta, eta_v, gamma, contact_matrix, npi_time_begin, npi_time_end, npi_cr, vax_time_begin, vax_time_end, vax_nu, rate_interventions, time_end, increment)
 }
 
 #' @title Compute the discrete probability of the truncated Erlang distribution

R/check_args_default.R

@@ -76,12 +76,52 @@
   # add null intervention and vaccination if these are missing
   # if not missing, check that they conform to expectations
   if (!"intervention" %in% names(mod_args)) {
-    mod_args[["intervention"]] <- no_intervention(
-      mod_args[["population"]]
+    # add as a list element named "contacts", and one named "beta"
+    mod_args[["intervention"]] <- list(
+      contacts = no_contacts_intervention(
+        mod_args[["population"]]
+      ),
+      # a dummy intervention on the rate parameter beta
+      beta = no_rate_intervention()
     )
   } else {
-    assert_intervention(mod_args[["intervention"]], mod_args[["population"]])
+    # if interventions are passed, check for the types and names
+    stopifnot(
+      "`intervention` must be a list of <intervention>s" =
+        checkmate::test_list(
+          mod_args[["intervention"]],
+          types = c("contacts_intervention", "rate_intervention"),
+          names = "unique"
+        )
+    )
+
+    # check for any other intervention list element names
+    checkmate::assert_names(
+      names(mod_args[["intervention"]]),
+      subset.of = c("beta", "gamma", "alpha", "contacts")
+    )
+
+    # if a contacts intervention is passed, check it
+    if ("contacts" %in% names(mod_args[["intervention"]])) {
+      # check the intervention on contacts
+      assert_intervention(
+        mod_args[["intervention"]][["contacts"]], "contacts",
+        mod_args[["population"]]
+      )
+    } else {
+      # if not contacts intervention is passed, add a dummy one
+      mod_args[["intervention"]]$contacts <- no_contacts_intervention(
+        mod_args[["population"]]
+      )
+    }
+
+    # if there is only an intervention on contacts, add a dummy intervention
+    # on the transmission rate beta
+    if (identical(names(mod_args[["intervention"]]), "contacts")) {
+      mod_args[["intervention"]]$beta <- no_rate_intervention()
+    }
   }
+
   if (!"vaccination" %in% names(mod_args)) {
     mod_args[["vaccination"]] <- no_vaccination(
       mod_args[["population"]]
@@ -125,9 +165,15 @@
     get_parameter(mod_args[["infection"]], "infectious_period")
 
   # get NPI related times and contact reductions
-  npi_time_begin <- get_parameter(mod_args[["intervention"]], "time_begin")
-  npi_time_end <- get_parameter(mod_args[["intervention"]], "time_end")
-  npi_cr <- get_parameter(mod_args[["intervention"]], "contact_reduction")
+  contact_interventions <- mod_args[["intervention"]][["contacts"]]
+  npi_time_begin <- get_parameter(contact_interventions, "time_begin")
+  npi_time_end <- get_parameter(contact_interventions, "time_end")
+  npi_cr <- get_parameter(contact_interventions, "reduction")
+
+  # get other interventions if any
+  rate_interventions <- mod_args[["intervention"]][
+    setdiff(names(mod_args[["intervention"]]), "contacts")
+  ]
 
   # get vaccination related times and rates
   vax_time_begin <- get_parameter(mod_args[["vaccination"]], "time_begin")
@@ -143,6 +189,7 @@
     npi_cr = npi_cr,
     vax_time_begin = vax_time_begin, vax_time_end = vax_time_end,
     vax_nu = vax_nu,
+    rate_interventions = rate_interventions,
     time_end = mod_args$time_end,
     increment = mod_args$increment
   )

R/check_args_vacamole.R

@@ -108,11 +108,53 @@
   # add null intervention if this is missing
   # if not missing, check that it conforms to expectations
   if (!"intervention" %in% names(mod_args)) {
-    mod_args[["intervention"]] <- no_intervention(
-      mod_args[["population"]]
+    # add as a list element named "contacts", and one named "beta"
+    mod_args[["intervention"]] <- list(
+      contacts = no_contacts_intervention(
+        mod_args[["population"]]
+      ),
+      # a dummy intervention on the rate parameter beta
+      beta = no_rate_intervention()
     )
   } else {
-    assert_intervention(mod_args[["intervention"]], mod_args[["population"]])
+    # if interventions are passed, check for the types and names
+    stopifnot(
+      "`intervention` must be a list of <intervention>s" =
+        checkmate::test_list(
+          mod_args[["intervention"]],
+          types = c("contacts_intervention", "rate_intervention"),
+          names = "unique"
+        )
+    )
+
+    # check for any other intervention list element names
+    checkmate::assert_names(
+      names(mod_args[["intervention"]]),
+      subset.of = c(
+        "beta", "beta_v", "gamma", "alpha", "eta", "eta_v",
+        "omega", "omega_v", "contacts"
+      )
+    )
+
+    # if a contacts intervention is passed, check it
+    if ("contacts" %in% names(mod_args[["intervention"]])) {
+      # check the intervention on contacts
+      assert_intervention(
+        mod_args[["intervention"]][["contacts"]], "contacts",
+        mod_args[["population"]]
+      )
+    } else {
+      # if not contacts intervention is passed, add a dummy one
+      mod_args[["intervention"]]$contacts <- no_contacts_intervention(
+        mod_args[["population"]]
+      )
+    }
+
+    # if there is only an intervention on contacts, add a dummy intervention
+    # on the transmission rate beta
+    if (identical(names(mod_args[["intervention"]]), "contacts")) {
+      mod_args[["intervention"]]$beta <- no_rate_intervention()
+    }
   }
 
   # return arguments invisibly
@@ -158,9 +200,15 @@
     (1.0 - get_parameter(mod_args[["infection"]], "hosp_reduction_vax"))
 
   # get NPI related times and contact reductions
-  npi_time_begin <- get_parameter(mod_args[["intervention"]], "time_begin")
-  npi_time_end <- get_parameter(mod_args[["intervention"]], "time_end")
-  npi_cr <- get_parameter(mod_args[["intervention"]], "contact_reduction")
+  contact_interventions <- mod_args[["intervention"]][["contacts"]]
+  npi_time_begin <- get_parameter(contact_interventions, "time_begin")
+  npi_time_end <- get_parameter(contact_interventions, "time_end")
+  npi_cr <- get_parameter(contact_interventions, "reduction")
+
+  # get other interventions if any
+  rate_interventions <- mod_args[["intervention"]][
+    setdiff(names(mod_args[["intervention"]]), "contacts")
+  ]
 
   # get vaccination related times and rates
   vax_time_begin <- get_parameter(mod_args[["vaccination"]], "time_begin")
@@ -179,6 +227,7 @@
     npi_cr = npi_cr,
     vax_time_begin = vax_time_begin, vax_time_end = vax_time_end,
     vax_nu = vax_nu,
+    rate_interventions = rate_interventions,
     time_end = mod_args$time_end,
     increment = mod_args$increment
   )

R/epidemic_default.R

@@ -85,7 +85,7 @@ epidemic_default_cpp <- function(population,
 
   # check class add intervention and vaccination if not NULL
   if (!is.null(intervention)) {
-    checkmate::assert_class(intervention, "intervention")
+    checkmate::assert_list(intervention, types = "intervention")
     model_arguments[["intervention"]] <- intervention
   }
   if (!is.null(vaccination)) {
@@ -142,16 +142,25 @@ epidemic_default_cpp <- function(population,
     cr = params[["npi_cr"]]
   )
 
+  # modify parameters
+  infection_params <- params[c("beta", "alpha", "gamma")]
+
+  infection_params <- intervention_on_rates(
+    t = t,
+    interventions = params[["rate_interventions"]],
+    parameters = infection_params
+  )
+
   # modify the vaccination rate depending on the regime
   # the number of doses is already checked before passing
   current_nu <- params[["vax_nu"]] *
     ((params[["vax_time_begin"]] < t) &
       (params[["vax_time_end"]] > t))
 
   # calculate transitions
-  sToE <- (params[["beta"]] * y[, 1] * contact_matrix_ %*% y[, 3])
-  eToI <- params[["alpha"]] * y[, 2]
-  iToR <- params[["gamma"]] * y[, 3]
+  sToE <- (infection_params[["beta"]] * y[, 1] * contact_matrix_ %*% y[, 3])
+  eToI <- infection_params[["alpha"]] * y[, 2]
+  iToR <- infection_params[["gamma"]] * y[, 3]
   sToV <- current_nu * y[, 1]
 
   # define compartmental changes
@@ -194,7 +203,10 @@ epidemic_default_r <- function(population,
 
   # check class add intervention and vaccination if not NULL
   if (!is.null(intervention)) {
-    checkmate::assert_class(intervention, "intervention")
+    checkmate::assert_list(
+      intervention,
+      types = c("intervention", "list")
+    )
     model_arguments[["intervention"]] <- intervention
   }
   if (!is.null(vaccination)) {

R/epidemic_vacamole.R

@@ -120,7 +120,10 @@ epidemic_vacamole_cpp <- function(population,
 
   # check class add intervention and vaccination if not NULL
   if (!is.null(intervention)) {
-    checkmate::assert_class(intervention, "intervention")
+    checkmate::assert_list(
+      intervention,
+      types = c("contacts_intervention", "rate_intervention")
+    )
     model_arguments[["intervention"]] <- intervention
   }
 
@@ -181,6 +184,20 @@ epidemic_vacamole_cpp <- function(population,
     cr = params[["npi_cr"]]
   )
 
+  # modify parameters
+  infection_params <- params[
+    c(
+      "beta", "beta_v", "alpha", "eta", "eta_v",
+      "omega", "omega_v", "gamma"
+    )
+  ]
+
+  infection_params <- intervention_on_rates(
+    t = t,
+    interventions = params[["rate_interventions"]],
+    parameters = infection_params
+  )
+
   # modify the vaccination rate depending on the regime
   # the number of doses is already checked before passing
   # columns are: 1: rate of first dose, 2: rate of second dose
@@ -276,7 +293,10 @@ epidemic_vacamole_r <- function(population,
 
   # check class add intervention and vaccination if not NULL
   if (!is.null(intervention)) {
-    checkmate::assert_class(intervention, "intervention")
+    checkmate::assert_list(
+      intervention,
+      types = c("contacts_intervention", "rate_intervention")
+    )
     model_arguments[["intervention"]] <- intervention
   }
 

R/input_check_helpers.R

@@ -209,31 +209,47 @@ assert_vaccination <- function(x, doses, population) {
 #' function is for internal use in argument checking functions.
 #'
 #' @param x A [intervention] object.
+#' @param type A string for the type of intervention to check for. May be one of
+#' `"contacts"` or `"rate"`.
 #' @param population An optional argument which is a [population] object.
 #' When present, this is used to check whether the intervention object `x` has
-#' corresponding values of `contact_reduction` for each demographic group in
+#' corresponding values of `reduction` for each demographic group in
 #' `population`.
 #'
 #' @keywords internal
 #'
 #' @return Silently returns the `intervention` object `x`.
 #' Primarily called for its side effects of throwing errors when `x` does not
 #' meet certain requirements.
-assert_intervention <- function(x, population) {
+assert_intervention <- function(x, type = c("contacts", "rate"),
+                                population) {
   # check for input class
-  checkmate::assert_class(x, "intervention")
+  type <- match.arg(type, several.ok = FALSE)
+  switch(type,
+    contacts = checkmate::assert_class(x, "contacts_intervention"),
+    rate = checkmate::assert_class(x, "rate_intervention")
+  )
 
-  # check that the length of contact_reduction is the same as the number
-  # of demographic groups
-  n_demo_groups <- NULL
-  if (!missing(population)) {
-    n_demo_groups <- length(get_parameter(population, "demography_vector"))
+  if (type == "contacts") {
+    # check the population
+    # check that the length of reduction is the same as the number
+    # of demographic groups
+    n_demo_groups <- NULL
+    if (!missing(population)) {
+      n_demo_groups <- length(get_parameter(population, "demography_vector"))
+    }
+    checkmate::assert_matrix(
+      get_parameter(x, "reduction"),
+      mode = "numeric",
+      nrows = n_demo_groups
+    )
+  } else if (type == "rate") {
+    checkmate::assert_numeric(
+      get_parameter(x, "reduction"),
+      lower = 0.0, upper = 1.0,
+      len = length(get_parameter(x, "time_begin"))
+    )
   }
-  checkmate::assert_matrix(
-    get_parameter(x, "contact_reduction"),
-    mode = "numeric",
-    nrows = n_demo_groups
-  )
 
   # invisibly return x
   invisible(x)