Skip to content

Commit

Permalink
Diphtheria model cpp accepts param vecs, WIP #166
Browse files Browse the repository at this point in the history
  • Loading branch information
@pratikunterwegs
pratikunterwegs committed Mar 11, 2024
1 parent 46d380a commit 6875275
Show file tree
Hide file tree
Showing 2 changed files with 204 additions and 81 deletions.
201 changes: 134 additions & 67 deletions R/model_diphtheria.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,21 +11,10 @@
#' This model accommodates age or demographic structure and allows for a
#' 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
#' of each demographic group. See [population()].
#' @param transmissibility A single number for the rate at which individuals
#' move from the susceptible to the exposed compartment upon contact with an
#' infectious individual. Often denoted as \eqn{\beta}, with
#' \eqn{\beta = R_0 / \text{infectious period}}.
#' @param infectiousness_rate A single number for the rate at which individuals
#' move from the exposed to the infectious compartment. Often denoted as
#' \eqn{\sigma}, with \eqn{\sigma = 1.0 / \text{pre-infectious period}}.
#' This value does not depend upon the number of infectious individuals in the
#' population.
#' The model also allows for changes to the number of susceptibles in each age
#' group to model influxes or evacuations from camps.
#'
#' @inheritParams model_default
#' @param reporting_rate A single number for the proportion of infectious cases
#' that is reported; this is a precursor to hospitalisation as only reported
#' cases are hospitalised.
Expand All @@ -37,30 +26,18 @@
#' @param hosp_exit_rate A single number for the rate at which individuals are
#' discharged from hospital to enter the 'recovered' compartment.
#' This is calculated as 1 / time to discharge, denoted \eqn{\tau_2}.
#' @param recovery_rate A single number for the recovery rate, denoted
#' \eqn{\gamma}.
#' @param prop_vaccinated A numeric vector of the same length as the number of
#' demographic groups indicated the proportion of each group that is vaccinated.
#' These individuals are not included in the model dynamics.
#' @param intervention A named list of `<rate_intervention>` objects
#' representing optional pharmaceutical or non-pharmaceutical interventions
#' applied to the model parameters listed above.
#' @param time_dependence A named list where each name
#' is a model parameter, and each element is a function with
#' the first two arguments being the current simulation `time`, and `x`, a value
#' that is dependent on `time` (`x` represents a model parameter).
#' See **Details** for more information, as well as the vignette on time-
#' dependence \code{vignette("time_dependence", package = "epidemics")}.
#' @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
#'
#' ## Rcpp implementations
Expand All @@ -71,8 +48,13 @@
#' ## Model parameters
#'
#' This model only allows for single, population-wide rates transitions between
#' compartments. The default values are taken from Finger et al. (2019) where
#' possible.
#' compartments per model run.
#'
#' However, model parameters may be passed as numeric vectors. These vectors
#' must follow Tidyverse recycling rules: all vectors must have the same length,
#' or, vectors of length 1 will be recycled to the length of any other vector.
#'
#' The default values are taken from Finger et al. (2019) where possible:
#'
#' - Transmissibility (\eqn{\beta}, `transmissibility`): 0.8888889, assuming an
#' \eqn{R_0} of 4.0 and a total infectious period of 4.5 days.
Expand Down Expand Up @@ -109,8 +91,57 @@
#' Kucharski, A. J. (2019). Real-time analysis of the diphtheria outbreak in
#' forcibly displaced Myanmar nationals in Bangladesh. BMC Medicine, 17, 58.
#' \doi{10.1186/s12916-019-1288-7}.
#' @return A `<data.frame>` with the columns "time", "compartment", "age_group",
#' and "value".
#' @return A `data.frame` with the columns "time", "compartment", "age_group",
#' "value", and "run", giving the number of individuals per demographic group
#' in each compartment at each timestep in long (or "tidy") format, with "run"
#' indicating the unique parameter combination.
#' @examples
#' # create a dummy camp population with three age groups
#' # diphtheria model is SEIHR
#' # assume that most are susceptible, some infectious
#' # values taken from supplementary material in Finger et al. for the
#' # Kutupalong camp, rounded to the nearest 100
#' n_age_groups <- 3
#' demography_vector <- c(83000, 108200, 224600)
#' initial_conditions <- matrix(0, nrow = n_age_groups, ncol = 5)
#'
#' # set susceptibles and infectious
#' initial_conditions[, 1] <- demography_vector - 1
#' initial_conditions[, 3] <- rep(1, n_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
#' )
#'
#' # assume younger age groups are vaccinated
#' prop_vaccinated <- c(0.2, 0.10, 0.1)
#'
#' # run model for single, default parameter set
#' data <- model_diphtheria_cpp(
#' camp_pop,
#' prop_vaccinated = prop_vaccinated
#' )
#' head(data)
#' tail(data)
#'
#' # run model with increase in population
#' # create population change data
#' p <- list(
#' time = 70,
#' values = list(
#' c(1e4, 2e5, 1e5)
#' )
#' )
#'
#' data <- model_diphtheria_cpp(
#' camp_pop,
#' prop_vaccinated = prop_vaccinated,
#' population_change = p
#' )
#' head(data)
#' tail(data)
#' @export
model_diphtheria_cpp <- function(population,
transmissibility = 4.0 / 4.5,
Expand All @@ -130,18 +161,37 @@ model_diphtheria_cpp <- function(population,
increment = 1) {
# check class on required inputs
checkmate::assert_class(population, "population")
checkmate::assert_number(transmissibility, lower = 0, finite = TRUE)
checkmate::assert_number(infectiousness_rate, lower = 0, finite = TRUE)
checkmate::assert_number(recovery_rate, lower = 0, finite = TRUE)
checkmate::assert_numeric(transmissibility, lower = 0, finite = TRUE)
checkmate::assert_numeric(infectiousness_rate, lower = 0, finite = TRUE)
checkmate::assert_numeric(recovery_rate, lower = 0, finite = TRUE)
# reporting rate and prop_hosp are expected to be proportions bounded 0 - 1
checkmate::assert_number(reporting_rate, lower = 0, upper = 1)
checkmate::assert_number(prop_hosp, lower = 0, upper = 1)
checkmate::assert_numeric(reporting_rate, lower = 0, upper = 1)
checkmate::assert_numeric(prop_hosp, lower = 0, upper = 1)
# hospital entry and exit rate are very likely to be bounded 0 - 1 but
# are allowed to be > 1 for now
checkmate::assert_number(hosp_entry_rate, lower = 0, finite = TRUE)
checkmate::assert_number(hosp_exit_rate, lower = 0, finite = TRUE)
checkmate::assert_numeric(hosp_entry_rate, lower = 0, finite = TRUE)
checkmate::assert_numeric(hosp_exit_rate, lower = 0, finite = TRUE)

# check the time end and increment
# restrict increment to lower limit of 1e-6
checkmate::assert_integerish(time_end, lower = 0)
checkmate::assert_number(increment, lower = 1e-6, finite = TRUE)

# check all vector lengths are equal or 1L
params <- list(
transmissibility = transmissibility,
infectiousness_rate = infectiousness_rate,
recovery_rate = recovery_rate,
reporting_rate = reporting_rate,
prop_hosp = prop_hosp,
hosp_entry_rate = hosp_entry_rate,
hosp_exit_rate = hosp_exit_rate,
time_end = time_end
)
assert_recyclable(params)

# check that prop_vaccinated is the same length as demography_vector
# TODO: treat this as a composable element for vectorisation
checkmate::assert_numeric(
prop_vaccinated,
lower = 0, upper = 1,
Expand Down Expand Up @@ -182,43 +232,60 @@ model_diphtheria_cpp <- function(population,
)
}

# check the time end and increment
# restrict increment to lower limit of 1e-6
checkmate::assert_number(time_end, lower = 0, finite = TRUE)
checkmate::assert_number(increment, lower = 1e-6, finite = TRUE)

# collect population, infection, and model arguments passed as `...`
model_arguments <- list(
population = population,
transmissibility = transmissibility,
infectiousness_rate = infectiousness_rate,
recovery_rate = recovery_rate,
reporting_rate = reporting_rate,
prop_hosp = prop_hosp,
hosp_entry_rate = hosp_entry_rate, hosp_exit_rate = hosp_exit_rate,
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
# collect population, infection
# combine parameters and composable elements into a list of lists of mod args
params <- .recycle_vectors(params)
params <- .transpose_base(params)
model_arguments <- lapply(
params, function(x) {
c(
list(
population = population,
prop_vaccinated = prop_vaccinated,
intervention = intervention,
time_dependence = time_dependence,
pop_change_times = population_change[["time"]],
pop_change_values = population_change[["values"]],
increment = increment
),
x
)
}
)

# prepare checked arguments for function
# this necessary as check_args adds intervention and vaccination
# if missing
model_arguments <- .prepare_args_model_diphtheria(
.check_args_model_diphtheria(model_arguments)
# cross-check model arguments
# TODO: simplify to only check composable elements
model_arguments <- lapply(
model_arguments, function(l) {
.prepare_args_model_diphtheria(
.check_args_model_diphtheria(l)
)
}
)

# get compartment names
compartments <- c(
"susceptible", "exposed", "infectious", "hospitalised", "recovered"
)

# run model over arguments
output <- do.call(.model_diphtheria_cpp, model_arguments)
# run model over arguments, prepare output, and return list
# assign run number as list index - this is the specific combination of
# parameters
output <- Map(
model_arguments, seq_along(model_arguments),
f = function(l, i) {
output_ <- .output_to_df(
do.call(.model_diphtheria_cpp, l),
population = population,
compartments = compartments
)
output_$run <- i
output_
}
)
output <- data.table::rbindlist(output)

# prepare output and return
output_to_df(output, population, compartments)
# convert to data.frame and return
# TODO: parameters need to be pulled along
data.table::setDF(output)[]
}
Loading

0 comments on commit 6875275

Please sign in to comment.