Camp population changes for diphtheria model

R/RcppExports.R

@@ -62,15 +62,19 @@
 #' @param rate_interventions A named list of `<rate_intervention>` objects.
 #' @param time_dependence A named list of functions for parameter time
 #' dependence.
+#' @param pop_change_times A numeric vector of times and which the population
+#' of susceptibles changes.
+#' @param pop_change_values An Rcpp List of numeric vectors giving the value of
+#' changes to each demographic group at each change in population.
 #' @param time_end The end time of the simulation.
 #' @param increment The time increment of the simulation.
 #' @return A two element list, where the first element is a list of matrices
 #' whose elements correspond to the numbers of individuals in each compartment
 #' as specified in the initial conditions matrix.
 #' The second list element is a vector of timesteps.
 #' @keywords internal
-.model_diphtheria_cpp <- function(initial_state, transmissibility, infectiousness_rate, recovery_rate, reporting_rate, prop_hosp, hosp_entry_rate, hosp_exit_rate, rate_interventions, time_dependence, time_end = 100.0, increment = 1.0) {
-    .Call(`_epidemics_model_diphtheria_cpp_internal`, initial_state, transmissibility, infectiousness_rate, recovery_rate, reporting_rate, prop_hosp, hosp_entry_rate, hosp_exit_rate, rate_interventions, time_dependence, time_end, increment)
+.model_diphtheria_cpp <- function(initial_state, transmissibility, infectiousness_rate, recovery_rate, reporting_rate, prop_hosp, hosp_entry_rate, hosp_exit_rate, rate_interventions, time_dependence, pop_change_times, pop_change_values, time_end = 100.0, increment = 1.0) {
+    .Call(`_epidemics_model_diphtheria_cpp_internal`, initial_state, transmissibility, infectiousness_rate, recovery_rate, reporting_rate, prop_hosp, hosp_entry_rate, hosp_exit_rate, rate_interventions, time_dependence, pop_change_times, pop_change_values, time_end, increment)
 }
 
 #' @title Run the RIVM Vacamole model

R/check_args_diphtheria.R

@@ -26,6 +26,8 @@
 #' 'exposed' and 'exposed_vaccinated' to the 'infectious' and
 #' 'infectious_vaccinated' compartments;
 #'
+#' - `recovery_rate`: a single number for the recovery rate from the infection;
+#'
 #'  - `reporting_rate`: a single number for the proportion of infectious cases
 #' reported;
 #'
@@ -35,7 +37,14 @@
 #' - `hosp_entry_rate`, `hosp_exit_rate`: two numbers representing the rate of
 #' entry and exit from the 'hospitalised' compartment;
 #'
-#' - `recovery_rate`: a single number for the recovery rate from the infection;
+#' - `rate_interventions`: an Rcpp List giving the interventions on model
+#' parameters;
+#'
+#' - `time_dependence`: an Rcpp List giving the time-dependent effects on model
+#' parameters in the form of R functions;
+#'
+#' - `pop_change_times` and `pop_change_values`: the times and values of changes
+#' in the population of susceptibles;
 #'
 #'  - `time_end`, `increment`: two numbers for the time at which to end the
 #' simulation, and the value by which the simulation time
@@ -101,11 +110,53 @@
     )
   }
 
+  # handle population change mechanic, check for correct demography groups
+  # check for only times being NULL and expect values also NULL
+  if (is.null(mod_args[["pop_change_times"]])) {
+    mod_args[["pop_change_times"]] <- 0
+    mod_args[["pop_change_values"]] <- list(
+      rep(
+        0, length(get_parameter(
+          mod_args[["population"]], "demography_vector"
+        ))
+      )
+    )
+  } else {
+    checkmate::assert_numeric(
+      mod_args[["pop_change_times"]],
+      lower = 0, finite = TRUE,
+      min.len = 1
+    )
+    checkmate::assert_list(
+      mod_args[["pop_change_values"]],
+      any.missing = FALSE,
+      len = length(mod_args[["pop_change_times"]])
+    )
+    invisible(
+      lapply(
+        mod_args[["pop_change_values"]],
+        FUN = function(x) {
+          stopifnot(
+            "`population_change` `values` must be same length as demography" =
+              checkmate::test_numeric(
+                x,
+                len = length(
+                  get_parameter(
+                    mod_args[["population"]], "demography_vector"
+                  )
+                )
+              )
+          )
+        }
+      )
+    )
+  }
+
   # return arguments invisibly
   invisible(mod_args)
 }
 
-#' @title Prepare arguments for the Vacamole epidemic function
+#' @title Prepare arguments for the diphtheria model
 #' @name check_prepare_diphtheria_args
 #' @rdname check_prepare_diphtheria_args
 #' @keywords internal
@@ -117,6 +168,7 @@
 
   # modify initial state by proportion vaccinated
   # NOTE: this assumes that the first column is 'susceptible'
+  # NOTE: there is no explicit vaccinated compartment
   initial_state[, 1] <- initial_state[, 1] * (1 - mod_args[["prop_vaccinated"]])
 
   # return mod args without population and prop_vaccinated

R/model_diphtheria.R

@@ -9,8 +9,10 @@
 #' The model is based on Finger et al. (2019) and is intended to be used in the
 #' context of internally displaced people (IDP) or refugee camps.
 #' This model accommodates age or demographic structure and allows for a
-#' proportion of each demographic group to be vaccinated and to not contribute
-#' to the outbreak.
+#' proportion of each demographic group to be vaccinated at the start of the
+#' outbreak, and thus to not contribute to the outbreak.
+#' The model also allows for changes to the initial population size, to model
+#' influxes or evacuations from camps.
 #'
 #' @param population An object of the `population` class, which holds a
 #' population contact matrix, a demography vector, and the initial conditions
@@ -52,10 +54,16 @@
 #' @param time_end The maximum number of timesteps over which to run the model.
 #' Taken as days, with a default value of 100 days.
 #' @param increment The size of the time increment. Taken as days, with a
+#' @param population_change A two-element list, with elements named `"time"` and
+#' `"values"`, giving the times of population changes, and the corresponding
+#' changes in the population of each demographic group at those times.
+#' `"time"` must be a numeric vector, while `"values"` must be a list of the
+#' length of `"time"`, with each element a numeric vector of the same length as
+#' the number of demographic groups in `population`.
 #' default value of 1 day.
 #' @details
 #'
-#' ## R and Rcpp implementations
+#' ## Rcpp implementations
 #'
 #' `model_diphtheria_cpp()` is a wrapper function for [.model_diphtheria_cpp()],
 #' an internal C++ function that uses Boost _odeint_ solvers for an SEIHR model.
@@ -87,6 +95,15 @@
 #' - Recovery rate (\eqn{\gamma}, `recovery_rate`): 0.333, assuming an
 #' infectious period following symptoms, of 3 days.
 #'
+#' ## Modelling population changes
+#'
+#' This model allows changes to the number of susceptibles in each demographic
+#' group, to represent influxes or evacuations from the camp as would be
+#' expected in humanitarian relief situations.
+#' Users can specify the times and changes (to each demographic group) of
+#' changes using the `population_changes` argument, to examine the effect on
+#' outbreak dynamics.
+#'
 #' @references
 #' Finger, F., Funk, S., White, K., Siddiqui, M. R., Edmunds, W. J., &
 #' Kucharski, A. J. (2019). Real-time analysis of the diphtheria outbreak in
@@ -108,6 +125,7 @@ model_diphtheria_cpp <- function(population,
                                  ),
                                  intervention = NULL,
                                  time_dependence = NULL,
+                                 population_change = NULL,
                                  time_end = 100,
                                  increment = 1) {
   # check class on required inputs
@@ -152,6 +170,18 @@ model_diphtheria_cpp <- function(population,
     )
   )
 
+  # check population change and create
+  checkmate::assert_list(
+    population_change,
+    null.ok = TRUE, len = 2L, types = c("numeric", "list")
+  )
+  if (!is.null(population_change)) {
+    checkmate::assert_names(
+      names(population_change),
+      identical.to = c("time", "values")
+    )
+  }
+
   # check the time end and increment
   # restrict increment to lower limit of 1e-6
   checkmate::assert_number(time_end, lower = 0, finite = TRUE)
@@ -169,6 +199,8 @@ model_diphtheria_cpp <- function(population,
     prop_vaccinated = prop_vaccinated,
     intervention = intervention,
     time_dependence = time_dependence,
+    pop_change_times = population_change[["time"]],
+    pop_change_values = population_change[["values"]],
     time_end = time_end, increment = increment
   )
 

_pkgdown.yml

@@ -28,6 +28,7 @@ articles:
   contents:
   - vacamole
   - ebola_model
+  - diphtheria
   - benchmarking
 
 reference:

inst/WORDLIST

@@ -4,6 +4,7 @@ Altman
 Analogs
 BH
 BMC
+Bazar
 Bjørnstad
 CMD
 COVID
@@ -30,6 +31,7 @@ Katriona
 Klinkenberg
 Korthals
 Krzywinski
+Kutupalong
 LMIC
 Leeuwen
 Leung
@@ -80,6 +82,7 @@ epiparameter
 epirecipes
 etc
 finalsize
+geq
 ggplot
 gh
 github
@@ -88,6 +91,7 @@ kylieainslie
 md
 odeint
 odin
+olds
 org
 packagename
 pkg

inst/include/epidemic_diphtheria.h

@@ -30,6 +30,9 @@ struct epidemic_diphtheria {
   const std::unordered_map<std::string, intervention::rate_intervention>
       interventions;
   const Rcpp::List time_dependence;
+  /// NOTE: not marked const as this interferes with member function
+  /// TODO: diagnose issue and achieve functioning marking as const
+  population::population_change pop_change;
 
   /// @brief Constructor for the diphtheria epidemic struct
   /// @param model_params An unordered map of string-double pairs, with the
@@ -53,11 +56,13 @@ struct epidemic_diphtheria {
       const std::unordered_map<std::string, double>& model_params,
       const std::unordered_map<std::string, intervention::rate_intervention>&
           interventions,
-      const Rcpp::List& time_dependence)
+      const Rcpp::List& time_dependence,
+      const population::population_change& pop_change)
       : model_params(model_params),
         model_params_temp(model_params),
         interventions(interventions),
-        time_dependence(time_dependence) {}
+        time_dependence(time_dependence),
+        pop_change(pop_change) {}
 
   /// @brief Operator for the diphtheria model
   /// @param x The initial state of the population - rows represent age groups
@@ -85,6 +90,11 @@ struct epidemic_diphtheria {
     // 0|1|2|3|4
     // S|E|I|H|R
 
+    // calculate any population changes
+    const Eigen::ArrayXd current_pop_change =
+        pop_change.get_population_change(t);
+
+    // calculate total infections
     const double total_infections = (x.col(2).array()).sum();
 
     // compartmental transitions without accounting for stratified contacts
@@ -99,19 +109,22 @@ struct epidemic_diphtheria {
     Eigen::ArrayXd iToH =
         model_params_temp["prop_hosp"] * model_params_temp["reporting_rate"] *
         model_params_temp["hosp_entry_rate"] * x.col(2).array();
-    Eigen::ArrayXd iToR = model_params_temp["recovery_rate"] * x.col(2).array();
+
+    // recoveries are from the infectious/ed who are NOT hospitalised
+    Eigen::ArrayXd iToR =
+        model_params_temp["recovery_rate"] * (x.col(2).array() - iToH);
     Eigen::ArrayXd hToR =
         model_params_temp["hosp_exit_rate"] * x.col(3).array();
 
     // compartmental changes; note that there are no contacts
     // β: transmissibility; σ: infectiousness rate; γ: recovery rate
     // ν: reporting rate; τ1: 1 / time to hospitalisation;
     // τ2: 1 / time to discharge from hospital; η: prop. hospitalised
-    dxdt.col(0) = -sToE;        // -β*S*I
-    dxdt.col(1) = sToE - eToI;  // β*S*I - σ*E
-    dxdt.col(2) = eToI - iToR;  // σ*E - γ*I
-    dxdt.col(3) = iToH - hToR;  // τ1*η*ν*I - τ2*H
-    dxdt.col(4) = iToR + hToR;  // γ*I + τ2*H
+    dxdt.col(0) = -sToE + current_pop_change;  // -β*S*I + pop movements
+    dxdt.col(1) = sToE - eToI;                 // β*S*I - σ*E
+    dxdt.col(2) = eToI - iToR - iToH;          // σ*E - γ*I - τ1*η*ν*I
+    dxdt.col(3) = iToH - hToR;                 // τ1*η*ν*I - τ2*H
+    dxdt.col(4) = iToR + hToR;                 // γ*I + τ2*H
   }
 };
 

inst/include/population.h

@@ -3,13 +3,11 @@
 #ifndef INST_INCLUDE_POPULATION_H_
 #define INST_INCLUDE_POPULATION_H_
 
-// [[Rcpp::plugins(cpp14)]]
-// [[Rcpp::depends(BH)]]
-// [[Rcpp::depends(RcppEigen)]]
-
 // clang-format off
 #include <Rcpp.h>
 #include <RcppEigen.h>
+
+#include <vector>
 // clang-format on
 
 // add to namespace population
@@ -41,6 +39,53 @@ inline Rcpp::NumericMatrix get_initial_conditions(
   return initial_conditions;
 }
 
+/// @brief Hold parameters for changes to the model population.
+/// This is special functionality for the diphtheria model applied to
+/// humanitarian camps.
+struct population_change {
+  const Rcpp::NumericVector times;
+  const Rcpp::List value;  // note these are actually absolute values
+  const int n_demo_groups;
+
+  /// @brief Constructor for the population change struct
+  /// @param times The times at which the population changes. Changes apply to
+  /// the SUSCEPTIBLES compartment only; the assumption is that the camp is the
+  /// locus of the outbreak.
+  /// @param value The ABSOLUTE change in the population size; an Rcpp List
+  /// of Rcpp NumericVectors, each of the same length as the number of
+  /// demographic groups.
+  population_change(const Rcpp::NumericVector &times, const Rcpp::List &value,
+                    const int &n_demo_groups)
+      : times(times), value(value), n_demo_groups(n_demo_groups) {}
+
+  // member function for population change at time t
+  /// @brief Calculate the population change at time t
+  /// @param t The current timestep
+  /// @return An Eigen Array of population size changes.
+  const Eigen::ArrayXd get_population_change(const double &t);
+};
+
+/// @brief Calculate the population change at time t
+/// @param t The current timestep
+/// @return An Eigen Array of population size changes, which defaults to zeros
+/// if no population change is scheduled
+inline const Eigen::ArrayXd population_change::get_population_change(
+    const double &t) {
+  // empty value for timepoints of no population change
+  Eigen::ArrayXd pop_change(n_demo_groups);
+  pop_change.fill(0.0);
+
+  // crudely check for a match between current time and scheduled pop change
+  for (size_t i = 0; i < times.size(); i++) {
+    if (t >= times[i] && t < (times[i] + 1.0)) {
+      pop_change = Rcpp::as<Eigen::ArrayXd>(value[i]);
+      return pop_change;
+    }
+  }
+
+  return pop_change;
+}
+
 }  // namespace population
 
 #endif  // INST_INCLUDE_POPULATION_H_