diff --git a/R/RcppExports.R b/R/RcppExports.R index 1d1ad629..825e8cb9 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -62,6 +62,10 @@ #' @param rate_interventions A named list of `` 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 @@ -69,8 +73,8 @@ #' 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 diff --git a/R/check_args_diphtheria.R b/R/check_args_diphtheria.R index 2b535d78..75143460 100644 --- a/R/check_args_diphtheria.R +++ b/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 diff --git a/R/model_diphtheria.R b/R/model_diphtheria.R index 08b1d106..fe234464 100644 --- a/R/model_diphtheria.R +++ b/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 ) diff --git a/_pkgdown.yml b/_pkgdown.yml index ab7c5c8f..29d96ef6 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -28,6 +28,7 @@ articles: contents: - vacamole - ebola_model + - diphtheria - benchmarking reference: diff --git a/inst/WORDLIST b/inst/WORDLIST index ec21b8db..5283ee33 100644 --- a/inst/WORDLIST +++ b/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 diff --git a/inst/include/epidemic_diphtheria.h b/inst/include/epidemic_diphtheria.h index 3c6e3506..99c96a26 100644 --- a/inst/include/epidemic_diphtheria.h +++ b/inst/include/epidemic_diphtheria.h @@ -30,6 +30,9 @@ struct epidemic_diphtheria { const std::unordered_map 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& model_params, const std::unordered_map& 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,7 +109,10 @@ 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(); @@ -107,11 +120,11 @@ struct epidemic_diphtheria { // β: 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 } }; diff --git a/inst/include/population.h b/inst/include/population.h index e7c45280..5a03ccb4 100644 --- a/inst/include/population.h +++ b/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 #include + +#include // 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 ×, 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(value[i]); + return pop_change; + } + } + + return pop_change; +} + } // namespace population #endif // INST_INCLUDE_POPULATION_H_ diff --git a/man/check_prepare_diphtheria_args.Rd b/man/check_prepare_diphtheria_args.Rd index 59d17b70..192a2915 100644 --- a/man/check_prepare_diphtheria_args.Rd +++ b/man/check_prepare_diphtheria_args.Rd @@ -27,13 +27,19 @@ vaccinated individuals, respectively; \item \code{infectiousness_rate}: a single number for the transition rate from the 'exposed' and 'exposed_vaccinated' to the 'infectious' and 'infectious_vaccinated' compartments; +\item \code{recovery_rate}: a single number for the recovery rate from the infection; \item \code{reporting_rate}: a single number for the proportion of infectious cases reported; \item \code{prop_hosp}: a single number for the proportion of reported cases that need hospitalisation; \item \code{hosp_entry_rate}, \code{hosp_exit_rate}: two numbers representing the rate of entry and exit from the 'hospitalised' compartment; -\item \code{recovery_rate}: a single number for the recovery rate from the infection; +\item \code{rate_interventions}: an Rcpp List giving the interventions on model +parameters; +\item \code{time_dependence}: an Rcpp List giving the time-dependent effects on model +parameters in the form of R functions; +\item \code{pop_change_times} and \code{pop_change_values}: the times and values of changes +in the population of susceptibles; \item \code{time_end}, \code{increment}: two numbers for the time at which to end the simulation, and the value by which the simulation time is incremented. diff --git a/man/dot-model_diphtheria_cpp.Rd b/man/dot-model_diphtheria_cpp.Rd index 2c9b603c..a412d493 100644 --- a/man/dot-model_diphtheria_cpp.Rd +++ b/man/dot-model_diphtheria_cpp.Rd @@ -15,6 +15,8 @@ hosp_exit_rate, rate_interventions, time_dependence, + pop_change_times, + pop_change_values, time_end = 100, increment = 1 ) @@ -44,6 +46,12 @@ hospital, represented as 1 / time to discharge \eqn{\tau_2}.} \item{time_dependence}{A named list of functions for parameter time dependence.} +\item{pop_change_times}{A numeric vector of times and which the population +of susceptibles changes.} + +\item{pop_change_values}{An Rcpp List of numeric vectors giving the value of +changes to each demographic group at each change in population.} + \item{time_end}{The end time of the simulation.} \item{increment}{The time increment of the simulation.} diff --git a/man/model_diphtheria.Rd b/man/model_diphtheria.Rd index 58513dad..902bd0c4 100644 --- a/man/model_diphtheria.Rd +++ b/man/model_diphtheria.Rd @@ -17,6 +17,7 @@ model_diphtheria_cpp( prop_vaccinated = 0 * get_parameter(population, "demography_vector"), intervention = NULL, time_dependence = NULL, + population_change = NULL, time_end = 100, increment = 1 ) @@ -70,11 +71,18 @@ that is dependent on \code{time} (\code{x} represents a model parameter). See \strong{Details} for more information, as well as the vignette on time- dependence \code{vignette("time_dependence", package = "epidemics")}.} +\item{population_change}{A two-element list, with elements named \code{"time"} and +\code{"values"}, giving the times of population changes, and the corresponding +changes in the population of each demographic group at those times. +\code{"time"} must be a numeric vector, while \code{"values"} must be a list of the +length of \code{"time"}, with each element a numeric vector of the same length as +the number of demographic groups in \code{population}. +default value of 1 day.} + \item{time_end}{The maximum number of timesteps over which to run the model. Taken as days, with a default value of 100 days.} -\item{increment}{The size of the time increment. Taken as days, with a -default value of 1 day.} +\item{increment}{The size of the time increment. Taken as days, with a} } \value{ A \verb{} with the columns "time", "compartment", "age_group", @@ -87,11 +95,13 @@ compartmental ordinary differential equation model with the compartments 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. } \details{ -\subsection{R and Rcpp implementations}{ +\subsection{Rcpp implementations}{ \code{model_diphtheria_cpp()} is a wrapper function for \code{\link[=.model_diphtheria_cpp]{.model_diphtheria_cpp()}}, an internal C++ function that uses Boost \emph{odeint} solvers for an SEIHR model. @@ -119,6 +129,16 @@ individuals are discharged from hospital after 5 days. infectious period following symptoms, of 3 days. } } + +\subsection{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 \code{population_changes} argument, to examine the effect on +outbreak dynamics. +} } \references{ Finger, F., Funk, S., White, K., Siddiqui, M. R., Edmunds, W. J., & diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index f8795272..774b2372 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -38,8 +38,8 @@ BEGIN_RCPP END_RCPP } // model_diphtheria_cpp_internal -Rcpp::List model_diphtheria_cpp_internal(const Eigen::MatrixXd& initial_state, const double& transmissibility, const double& infectiousness_rate, const double& recovery_rate, const double& reporting_rate, const double& prop_hosp, const double& hosp_entry_rate, const double& hosp_exit_rate, const Rcpp::List& rate_interventions, const Rcpp::List& time_dependence, const double& time_end, const double& increment); -RcppExport SEXP _epidemics_model_diphtheria_cpp_internal(SEXP initial_stateSEXP, SEXP transmissibilitySEXP, SEXP infectiousness_rateSEXP, SEXP recovery_rateSEXP, SEXP reporting_rateSEXP, SEXP prop_hospSEXP, SEXP hosp_entry_rateSEXP, SEXP hosp_exit_rateSEXP, SEXP rate_interventionsSEXP, SEXP time_dependenceSEXP, SEXP time_endSEXP, SEXP incrementSEXP) { +Rcpp::List model_diphtheria_cpp_internal(const Eigen::MatrixXd& initial_state, const double& transmissibility, const double& infectiousness_rate, const double& recovery_rate, const double& reporting_rate, const double& prop_hosp, const double& hosp_entry_rate, const double& hosp_exit_rate, const Rcpp::List& rate_interventions, const Rcpp::List& time_dependence, const Rcpp::NumericVector& pop_change_times, const Rcpp::List& pop_change_values, const double& time_end, const double& increment); +RcppExport SEXP _epidemics_model_diphtheria_cpp_internal(SEXP initial_stateSEXP, SEXP transmissibilitySEXP, SEXP infectiousness_rateSEXP, SEXP recovery_rateSEXP, SEXP reporting_rateSEXP, SEXP prop_hospSEXP, SEXP hosp_entry_rateSEXP, SEXP hosp_exit_rateSEXP, SEXP rate_interventionsSEXP, SEXP time_dependenceSEXP, SEXP pop_change_timesSEXP, SEXP pop_change_valuesSEXP, SEXP time_endSEXP, SEXP incrementSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -53,9 +53,11 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const double& >::type hosp_exit_rate(hosp_exit_rateSEXP); Rcpp::traits::input_parameter< const Rcpp::List& >::type rate_interventions(rate_interventionsSEXP); Rcpp::traits::input_parameter< const Rcpp::List& >::type time_dependence(time_dependenceSEXP); + Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type pop_change_times(pop_change_timesSEXP); + Rcpp::traits::input_parameter< const Rcpp::List& >::type pop_change_values(pop_change_valuesSEXP); Rcpp::traits::input_parameter< const double& >::type time_end(time_endSEXP); Rcpp::traits::input_parameter< const double& >::type increment(incrementSEXP); - rcpp_result_gen = Rcpp::wrap(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)); + rcpp_result_gen = Rcpp::wrap(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)); return rcpp_result_gen; END_RCPP } @@ -92,7 +94,7 @@ END_RCPP static const R_CallMethodDef CallEntries[] = { {"_epidemics_model_default_cpp_internal", (DL_FUNC) &_epidemics_model_default_cpp_internal, 15}, - {"_epidemics_model_diphtheria_cpp_internal", (DL_FUNC) &_epidemics_model_diphtheria_cpp_internal, 12}, + {"_epidemics_model_diphtheria_cpp_internal", (DL_FUNC) &_epidemics_model_diphtheria_cpp_internal, 14}, {"_epidemics_model_vacamole_cpp_internal", (DL_FUNC) &_epidemics_model_vacamole_cpp_internal, 20}, {NULL, NULL, 0} }; diff --git a/src/epidemic_diphtheria.cpp b/src/epidemic_diphtheria.cpp index 788318e1..fd993d7d 100644 --- a/src/epidemic_diphtheria.cpp +++ b/src/epidemic_diphtheria.cpp @@ -32,6 +32,10 @@ //' @param rate_interventions A named list of `` 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 @@ -46,6 +50,8 @@ Rcpp::List model_diphtheria_cpp_internal( const double &reporting_rate, const double &prop_hosp, const double &hosp_entry_rate, const double &hosp_exit_rate, const Rcpp::List &rate_interventions, const Rcpp::List &time_dependence, + const Rcpp::NumericVector &pop_change_times, + const Rcpp::List &pop_change_values, const double &time_end = 100.0, // double required by boost solver const double &increment = 1.0) { // initial conditions from input @@ -66,9 +72,13 @@ Rcpp::List model_diphtheria_cpp_internal( rate_interventions_cpp = intervention::rate_intervention_cpp(rate_interventions); + // prepare population changes if any + const population::population_change pop_change( + pop_change_times, pop_change_values, initial_state.rows()); + // create a diphtheria epidemic with parameters epidemics::epidemic_diphtheria this_model( - model_params, rate_interventions_cpp, time_dependence); + model_params, rate_interventions_cpp, time_dependence, pop_change); // prepare storage containers for the observer std::vector x_vec; // is a vector of MatrixXd diff --git a/tests/testthat/test-model_diphtheria.R b/tests/testthat/test-model_diphtheria.R index 6f781c82..176db9ff 100644 --- a/tests/testthat/test-model_diphtheria.R +++ b/tests/testthat/test-model_diphtheria.R @@ -52,11 +52,11 @@ test_that("Diptheria model, basic expectations", { # check for identical numbers of individuals at start and end # NOTE: high tolerance because hospitalised compartment is not directly # linked to infectious compartment per Finger et al. model structure. - # leads to more individuals at final state than initial state + # leads to different individuals at final state than initial state expect_identical( sum(output[output$time == min(output$time), ]$value), sum(output[output$time == max(output$time), ]$value), - tolerance = 100 + tolerance = 1e-6 ) # check that all age groups in the simulation are the same # size as the demography vector @@ -67,3 +67,168 @@ test_that("Diptheria model, basic expectations", { # NOTE: no checks for final state equal to demography vector as # vaccinated individuals are removed from model }) + +# Expectations for the diphtheria model with changed population sizes +test_that("Diphtheria model with population size changes", { + # model runs with population change functionality + # create population change data + p <- list( + time = 70, + values = list( + c(1e4, 2e5, 1e5) + ) + ) + + expect_no_condition( + model_diphtheria_cpp( + population = camp_pop, + population_change = p + ) + ) + + # NOTE: expected final population size is larger than the initial + # but identical to the original + added population + data <- model_diphtheria_cpp( + population = camp_pop, + prop_hosp = 0.08, + population_change = p, + time_end = 200 + ) + + last_value <- aggregate( + value ~ demography_group, + data = data[data$time == max(data$time), ], FUN = "sum" + ) + last_value + + expect_identical( + last_value$value, + camp_pop$demography_vector + p$values[[1]], + tolerance = 1e-6 + ) + + ## Multiple population changes including decreases + p <- list( + time = c(70, 80), + values = list( + c(1e4, 2e5, 1e5), # influx to camp + c(-9e3, -1.5e5, -0.5e5) # evacuation from camp + ) + ) + + expect_no_condition( + model_diphtheria_cpp( + population = camp_pop, + population_change = p + ) + ) + + # NOTE: expected final population size is larger than the initial + # but identical to the original + net added population + data <- model_diphtheria_cpp( + population = camp_pop, + prop_hosp = 0.08, + population_change = p, + time_end = 200 + ) + + last_value <- aggregate( + value ~ demography_group, + data = data[data$time == max(data$time), ], FUN = "sum" + ) + last_value + + expect_identical( + last_value$value, + camp_pop$demography_vector + Reduce(x = p$values, f = `+`), + tolerance = 1e-6 + ) +}) + +# Test for poor specification of the population change mechanic +test_that("Diphtheria model handles population_change errors", { + # badly specified pop change + p <- list( + time = "some time", + values = list( + c(1, 2, 3) + ) + ) + expect_error( + model_diphtheria_cpp( + camp_pop, + population_change = p + ), + regexp = "May only contain the following types: \\{numeric,list\\}" + ) + + # wrong name of first element + p <- list( + timesteps = 10, + values = list( + c(1, 2, 3) + ) + ) + expect_error( + model_diphtheria_cpp( + camp_pop, + population_change = p + ), + regexp = "(Names must be identical)*(time)*(values)" + ) + + # wrong length of values - must be same length as demography groups + p <- list( + time = 10, + values = list( + c(1, 2) + ) + ) + expect_error( + model_diphtheria_cpp( + camp_pop, + population_change = p + ), + regexp = "`population_change` `values` must be same length as demography" + ) +}) + +# Tests for expected failures/input checks on interventions and time-dependence +test_that("Diphtheria model input checks", { + # expect no contacts interventions allowed + expect_error( + model_diphtheria_cpp( + camp_pop, + intervention = list( + contacts = no_contacts_intervention(camp_pop) # a dummy intervention + ) + ), + regexp = "May only contain the following types: \\{rate_intervention\\}" + ) + + # expect failure on rate interventions targeting disallowed parameters + expect_error( + model_diphtheria_cpp( + camp_pop, + intervention = list( + beta = intervention( + type = "rate", + time_begin = 10, time_end = 20, + reduction = 0.1 + ) + ) + ), + regexp = "(must be a subset of)*(but has additional elements)" + ) + + # expect failure on time dependence targeting disallowed parameters + expect_error( + model_diphtheria_cpp( + camp_pop, + time_dependence = list( + beta = function(time, x) x + ) + ), + regexp = "(must be a subset of)*(but has additional elements)" + ) +}) diff --git a/vignettes/diphtheria.Rmd b/vignettes/diphtheria.Rmd new file mode 100644 index 00000000..5f538d66 --- /dev/null +++ b/vignettes/diphtheria.Rmd @@ -0,0 +1,221 @@ +--- +title: "Modelling a diphtheria outbreak in a humanitarian camp setting" +output: + bookdown::html_vignette2: + fig_caption: yes + code_folding: show +pkgdown: + as_is: true +bibliography: references.json +link-citations: true +vignette: > + %\VignetteIndexEntry{Modelling a diphtheria outbreak in a humanitarian camp setting} + %\VignetteEncoding{UTF-8} + %\VignetteEngine{knitr::rmarkdown} +editor_options: + chunk_output_type: console +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + message = FALSE, + warning = FALSE, + fig.width = 5, + fig.height = 4, + dpi = 150 +) +``` + +::: {.alert .alert-warning} +**New to _epidemics_?** It may help to read the ["Get started"](epidemics.html) vignette first! +::: + +This vignette shows how to model an outbreak of diphtheria (or a similar acute directly-transmitted infectious disease) within the setting of a humanitarian aid camp, where the camp population may fluctuate due to external factors, including influxes from crisis affected areas or evacuations to other camps or areas. +In such situations, implementing large-scale public health measures such contact tracing and quarantine, or introducing reactive mass vaccination may be challenging. + +::: {.alert .alert-info} +The [vignette on modelling a vaccination campaign](modelling_vaccination.html) shows how to model the introduction of a mass vaccination campaign for a fixed, stable population. +::: + +_epidemics_ provides a simple SEIHR compartmental model based on @finger2019, in which it is possible to vary the population of each demographic group throughout the model's simulation time and explore the resulting epidemic dynamics. +This baseline model only tracks infections by demographic groups, and does not include variation in contacts between demographic groups (e.g. by age or occupation), as contacts are likely to be less clearly stratified in a camp setting. +The baseline model also does not allow interventions that target social contacts (such as social distancing or quarantine), and does not include a 'vaccinated' compartment. It is therefore suited to analysis of a rapidly spreading infection prior to the introduction of any reactive vaccine campaign. + +However, the model does allow for seasonality in model parameters and interventions on model parameters. +Similarly, the model allows for a proportion of the initial camp population to be considered vaccinated and thus immune from infection. + +```{r setup} +library(epidemics) +library(dplyr) +library(ggplot2) +``` + +## Modelling an outbreak with pre-existing immunity + +We create a population object corresponding to the Kutupalong camp in Cox's Bazar, Bangladesh in 2017-18, rounded to the nearest 100, as described in Additional file 1 provided with @finger2019. +This population has three age groups, < 5 years, 5 -- 14 years, $\geq$ 15 years. +We assume that only one individual is infectious in each age group. + +```{r camp-population} +# three age groups with five compartments SEIHR +n_age_groups <- 3 +demography_vector <- c(83000, 108200, 224600) +initial_conditions <- matrix(0, nrow = n_age_groups, ncol = 5) + +# 1 individual in each group is infectious +initial_conditions[, 1] <- demography_vector - 1 +initial_conditions[, 3] <- rep(1, n_age_groups) + +# camp social contact rates are assumed to be uniform within and between +# age groups +camp_pop <- population( + contact_matrix = matrix(1, nrow = n_age_groups, ncol = n_age_groups), + demography_vector = demography_vector, + initial_conditions = initial_conditions / demography_vector +) + +camp_pop +``` + +We assume, following @finger2019, that 20% of the 5 -- 14 year-olds were vaccinated against (and immune to) diphtheria prior to the start of the outbreak, but that coverage is much lower among other age groups. + +```{r} +# 20% of 5-14 year-olds are vaccinated +prop_vaccinated <- c(0.05, 0.2, 0.05) +``` + +We run the model with its default parameters, assuming that: + +- diphtheria has an $R_0$ of 4.0 and a mean infectious period of 4.5 days, giving a transmissibility ($\beta$) of about 0.889; +- diphtheria has a pre-infectious or incubation period of 3 days, giving an infectiousness rate ($\sigma$) of about 0.33; and +- the recovery rate of diphtheria is about 0.33. + +We also make several assumptions regarding clinical progression and reporting. Specifically: case reporting (that about 3% of infections are reported as cases), the proportion of reported cases needing hospitalisation (1%), the time taken by cases needing hospitalisation to seek and be admitted to hospital (5 days, giving a daily hospitalisation rate of 0.2 among cases), and time spent in hospital (5 days, giving a daily hospitalisation recovery rate of 0.2). + +Finally, we assume there are no interventions or seasonal effects that affect the dynamics of transmission during the outbreak. We then run the model and plot the outcomes. + +```{r diphtheria-basic} +data <- model_diphtheria_cpp( + population = camp_pop, + prop_vaccinated = prop_vaccinated +) +``` + +```{r class.source = 'fold-hide', fig.cap="Model results from a single run showing the number of individuals infectious with diphtheria over 100 days of the outbreak."} +filter(data, compartment == "infectious") |> + ggplot() + + geom_line( + aes(time, value, colour = demography_group) + ) + + scale_y_continuous( + labels = scales::comma + ) + + scale_colour_brewer( + palette = "Dark2", + name = "Age group", + labels = c("<5", "5-15", ">15") + ) + + expand_limits( + x = c(0, 101) + ) + + theme_bw() + + theme( + legend.position = "top" + ) + + labs( + x = "Simulation time (days)", + y = "Individuals infectious" + ) +``` + +## Modelling an outbreak with changing population sizes + +We now model the same outbreak, but an increase in susceptible individuals towards the end of the outbreak, to illustrate the effect that an influx of non-immune individuals could could have on outbreak dynamics. +In this example, we do not assume any prior immunity among new arrivals into the population. + +We prepare a population change schedule as a named list giving the times of each change, and the corresponding changes to each demographic group. + +Note that the model assumes that these changes apply _only to the susceptible compartment_, as we assume that the wider population entering the camp is not yet affected by diphtheria (i.e. no infected or recovered arrival), and that already infected or hospitalised individuals do not leave the camp. + +```{r pop_change} +# susceptibles increase by about 12%, 92%, and 89% of initial sizes +pop_change <- list( + time = 70, + values = list( + c(1e4, 1e5, 2e5) + ) +) +``` + +Here, the population size of the camp increases by about 75% overall, which is similar to reported values in the Kutupalong camp scenario [@finger2019]. + +```{r} +data <- model_diphtheria_cpp( + population = camp_pop, + population_change = pop_change +) + +# summarise population change in susceptibles +data_pop_size <- filter(data, compartment == "susceptible") |> + group_by(time) |> + summarise( + total_susceptibles = sum(value) + ) +``` + +```{r class.source = 'fold-hide', fig.cap="Model results from a single run showing the number of individuals infectious with diphtheria over 100 days of the outbreak, with an increase in the camp population size. Shaded blue region shows the number of individuals susceptible to infection (right-hand side Y axis)."} +ggplot() + + geom_area( + data = data_pop_size, + aes(time, total_susceptibles / 10), + fill = "steelblue", alpha = 0.5 + ) + + geom_line( + data = filter(data, compartment == "infectious"), + aes(time, value, colour = demography_group) + ) + + geom_vline( + xintercept = pop_change$time, + colour = "red", + linetype = "dashed", + linewidth = 0.2 + ) + + annotate( + geom = "text", + x = pop_change$time, + y = 20e3, + label = "Population increase", + angle = 90, + vjust = "inward", + colour = "red" + ) + + scale_y_continuous( + labels = scales::comma, + sec.axis = dup_axis( + trans = function(x) x * 10, + name = "Individuals susceptible" + ) + ) + + scale_colour_brewer( + palette = "Dark2", + name = "Age group (individuals infectious)", + labels = c("<5", "5-15", ">15") + ) + + expand_limits( + x = c(0, 101) + ) + + theme_bw() + + theme( + legend.position = "top" + ) + + labs( + x = "Simulation time (days)", + y = "Individuals infectious" + ) +``` + +This example shows how an increase in the number of susceptibles in the population can lead to a rise in the transmission potential of the infection and therefore extend the duration of an outbreak. + +## References diff --git a/vignettes/ebola_model.Rmd b/vignettes/ebola_model.Rmd index 5a59520e..eae690ea 100644 --- a/vignettes/ebola_model.Rmd +++ b/vignettes/ebola_model.Rmd @@ -41,7 +41,6 @@ The transitions between compartments are based on an Erlang model developed by @ The model does not include demographic variation in contacts, as ebola spreads primarily among contacts caring for infected individuals, making age or demographic structure less important. This model can currently accommodate interventions on model rates only (as there are no demographic groups for contacts interventions). -The only rate in the model that can be influenced by an intervention is the transmission rate $\beta$. ```{r setup} # some initial setup to load necessary packages