Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Remove infection class #140

Merged
merged 29 commits into from
Nov 21, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
29 commits
Select commit Hold shift + click to select a range
cb34045
Remove infection-class related files
pratikunterwegs Nov 8, 2023
60f97d8
Remove infection-class C++ header
pratikunterwegs Nov 8, 2023
ff0cfcd
Add model args in fn call for default model
pratikunterwegs Nov 8, 2023
cc8bd27
Add model args in fn call for Vacamole model
pratikunterwegs Nov 8, 2023
7767e14
Refactor ebola model to not use infection class
pratikunterwegs Nov 8, 2023
21cf2ca
Use renamed rate params in model headers
pratikunterwegs Nov 8, 2023
3cd502a
Default and Vacamole args checker fns use rates not infection
pratikunterwegs Nov 8, 2023
e07dc3d
Default no time dependence uses renamed rate
pratikunterwegs Nov 8, 2023
048a0f0
Model C++ src files use renamed rates
pratikunterwegs Nov 8, 2023
428513c
Test files do not use infection-class
pratikunterwegs Nov 8, 2023
304f0a2
Update function documentation with direct params
pratikunterwegs Nov 9, 2023
6e0d468
Update vignettes using default model
pratikunterwegs Nov 9, 2023
ec2d390
Correct parameter naming in Vacamole header
pratikunterwegs Nov 9, 2023
f408578
Update Vacamole and {finalsize} vignettes
pratikunterwegs Nov 9, 2023
3f1cf2a
Update ebola vignette
pratikunterwegs Nov 9, 2023
221aab9
Update helper fn get_parameters()
pratikunterwegs Nov 9, 2023
8b31c61
Update NAMESPACE
pratikunterwegs Nov 9, 2023
1df42a6
Update ebola model docs with new parameters
pratikunterwegs Nov 9, 2023
b68aa3e
Update RcppExports
pratikunterwegs Nov 9, 2023
68fc631
Remove ebola model C++ implementation
pratikunterwegs Nov 9, 2023
c75d5fe
Clean up ebola-Cpp related code
pratikunterwegs Nov 9, 2023
587f3f0
Remove test for prob_erlang_cpp
pratikunterwegs Nov 9, 2023
1f6b845
Update parameter uncertainty vignette
pratikunterwegs Nov 9, 2023
3cc7680
Remove infection from time_dep vignette
pratikunterwegs Nov 10, 2023
8bc2606
Update and add function examples
pratikunterwegs Nov 10, 2023
94c39c6
Update documentation and wordlist
pratikunterwegs Nov 10, 2023
347be80
Remove unicode characters from fn docs
pratikunterwegs Nov 10, 2023
f09f514
Make prob_discrete_erlang internal
pratikunterwegs Nov 10, 2023
bb92bc5
Remove <infection> reference from README.Rmd
bahadzie Nov 14, 2023
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 0 additions & 5 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -3,19 +3,15 @@
S3method(c,contacts_intervention)
S3method(c,rate_intervention)
S3method(c,vaccination)
S3method(print,infection)
S3method(print,intervention)
S3method(print,population)
S3method(print,vaccination)
export(as.intervention)
export(as.vaccination)
export(epidemic_size)
export(get_parameter)
export(get_transmission_rate)
export(infection)
export(intervention)
export(is_contacts_intervention)
export(is_infection)
export(is_intervention)
export(is_population)
export(is_rate_intervention)
Expand All @@ -28,7 +24,6 @@ export(model_vacamole_r)
export(new_infections)
export(no_contacts_intervention)
export(no_rate_intervention)
export(no_time_dependence)
export(no_vaccination)
export(population)
export(vaccination)
Expand Down
165 changes: 57 additions & 108 deletions R/RcppExports.R
Original file line number Diff line number Diff line change
Expand Up @@ -9,72 +9,34 @@
#' This function is intended to only be called internally from
#' [model_default_cpp()].
#'
#' Allows heterogeneity in social contact patterns, and variable sizes of
#' demographic groups.
#' Also allows for group-specific initial proportions in each model
#' compartment, as well as group-specific vaccination start dates and
#' vaccination rates, and also group-specific effects of implementing a
#' non-pharmaceutical intervention.
#' The model only allows for single, population-wide rates of
#' transition between the 'susceptible' and 'exposed' compartments, between the
#' 'exposed' and 'infectious' compartments, and in the recovery rate.
#'
#' @param initial_state An `Eigen::MatrixXd` holding the initial state of
#' each demographic-compartmental combination. Rows must represent demographic
#' groups, while columns represent compartments in the order S, E, I, R, V.
#' @param contact_matrix An `Eigen::MatrixXd` holding the population
#' contact matrix.
#' @param beta The transmission rate \eqn{\beta}.
#' @param alpha The rate of transition from exposed to infectious \eqn{\alpha}.
#' @param gamma The recovery rate \eqn{\gamma}.
#' @param time_end The maximum time, defaults to 200.0.
#' @param intervention A non-pharmaceutical intervention applied during the
#' course of the epidemic, with a start and end time, and age-specific effect
#' on contacts. See [intervention()].
#' @param vaccination A vaccination regime followed during the
#' course of the epidemic, with a start and end time, and age-specific effect
#' on the transition of individuals from susceptible to vaccinated.
#' See [vaccination()].
#' @param increment The increment time, defaults to 0.1.
#' @param initial_state A matrix for the initial state of the compartments.
#' @param transmissibility The transmission rate \eqn{\beta} at which
#' unvaccinated and partially vaccinated individuals are infected by the
#' disease.
#' @param infectiousness_rate The rate of transition from exposed to infectious
#' \eqn{\alpha}.
#' @param recovery_rate The recovery rate \eqn{\gamma}.
#' @param contact_matrix The population contact matrix.
#' @param npi_time_begin The start time of any non-pharmaceutical interventions
#' .
#' @param npi_time_end The end time of any non-pharmaceutical interventions.
#' @param npi_cr The reduction in contacts from any non-pharmaceutical
#' interventions.
#' @param vax_time_begin The start time of any vaccination campaigns.
#' @param vax_time_end The end time of any vaccination campaigns.
#' @param vax_nu The vaccination rate of any vaccination campaigns.
#' @param rate_interventions A named list of `<rate_intervention>` objects.
#' @param time_dependence A named list of functions for parameter time
#' dependence.
#' @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 (see [population()]).
#' The second list element is a vector of timesteps.
#' @keywords internal
.model_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_dependence, time_end = 100.0, increment = 1.0) {
.Call(`_epidemics_model_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_dependence, time_end, increment)
}

#' @title Run an SEIR model with Erlang passage times
#'
#' @description Wrapper function for an SEIR compartmental model with Erlang
#' passage times based on Getz and Dougherty (2017) J. Biological Dynamics,
#' and developed to model the West African ebola virus disease outbreak of 2014
#' .
#'
#' @param initial_conditions An integer vector of the initial conditions
#' corresponding to an SEIR model.
#' @param population_size A single integer for the population size.
#' @param beta The transmission rate \eqn{\beta}.
#' @param shape_E A single integer for the shape parameter of the Erlang
#' distribution of passage times through the exposed compartment.
#' @param rate_E A single double for the rate parameter of the Erlang
#' distribution of passage times through the exposed compartment.
#' @param shape_I A single integer for the shape parameter of the Erlang
#' distribution of passage times through the infectious compartment.
#' @param rate_I A single double for the rate parameter of the Erlang
#' distribution of passage times through the infectious compartment.
#' @param time_end A single integer for the maximum simulation time; this is
#' assumed to be in days.
#' @return A list with two elements, `x`, and integer matrix with as many rows
#' as the number of timesteps (given by `time_end`), and four columns, one for
#' each compartment, susceptible, exposed, infectious, recovered, in that order
#' ; and `times`, a vector of the simulation times, taken to be days.
#' This output is intended to be passed to [output_to_df()] to be converted
#' into a data.frame for further analysis.
#' @keywords internal
.model_ebola_cpp <- function(initial_state, population_size, beta, shape_E, rate_E, shape_I, rate_I, time_end) {
.Call(`_epidemics_model_ebola_cpp_internal`, initial_state, population_size, beta, shape_E, rate_E, shape_I, rate_I, time_end)
.model_default_cpp <- function(initial_state, transmissibility, infectiousness_rate, recovery_rate, contact_matrix, npi_time_begin, npi_time_end, npi_cr, vax_time_begin, vax_time_end, vax_nu, rate_interventions, time_dependence, time_end = 100.0, increment = 1.0) {
.Call(`_epidemics_model_default_cpp_internal`, initial_state, transmissibility, infectiousness_rate, recovery_rate, contact_matrix, npi_time_begin, npi_time_end, npi_cr, vax_time_begin, vax_time_end, vax_nu, rate_interventions, time_dependence, time_end, increment)
}

#' @title Run the RIVM Vacamole model
Expand All @@ -86,61 +48,48 @@
#' Model code: https://github.com/kylieainslie/vacamole
#' Manuscript describing the model and its application:
#' https://doi.org/10.2807/1560-7917.ES.2022.27.44.2101090
#' This function is intended to only be called internally from
#' [model_vacamole_cpp()].
#'
#' @details The original model has 8 conceptual compartments - four
#' epidemiological compartments (SEIR), three hospitalisation compartments
#' (H, ICU, ICU2H), and death - see the manuscript in Eurosurveillance.
#' Only infected individuals can enter the hospitalisation or death
#' compartments.
#' Vacamole was implemented as a stand-alone R package, and some versions have
#' been used to generate scenarios for the ECDC Covid-19 Scenario Hub.
#'
#' Individuals from the susceptible compartment may be vaccinated partially
#' or fully (assuming a two dose regimen), with only the second dose reducing
#' their probability of being infected, and of being hospitalised or dying.
#'
#' @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 beta The transmission rate \eqn{\beta} at which unvaccinated and
#' partially vaccinated individuals are infected by the disease.
#' @param beta_v The transmission rate \eqn{\beta_V} at which individuals who
#' have received two vaccine doses are infected by the disease.
#' @param alpha The rate of transition from exposed to infectious \eqn{\alpha}.
#' @param initial_state A matrix for the initial state of the compartments.
#' @param transmissibility The transmission rate \eqn{\beta} at which
#' unvaccinated and partially vaccinated individuals are infected by the
#' disease.
#' @param transmissibility_vax The transmission rate \eqn{\beta_V} at which
#' individuals who have received two vaccine doses are infected by the disease.
#' @param infectiousness_rate The rate of transition from exposed to infectious
#' \eqn{\alpha}.
#' This is common to fully susceptible, partially vaccinated, and fully
#' vaccinated individuals (where fully vaccinated represents two doses).
#' @param omega The mortality rate of fully susceptible and partially
#' vaccinated and unprotected individuals.
#' @param omega_v The mortality rate of individuals who are protected by
#' vaccination.
#' @param eta The hospitalisation rate of fully susceptible and partially
#' @param mortality_rate The mortality rate of fully susceptible and partially
#' vaccinated and unprotected individuals.
#' @param eta_v The hospitalisation rate of individuals who are protected by
#' vaccination.
#' @param gamma The recovery rate \eqn{\gamma}.
#' @param time_end The maximum time. Defaults to 100.0.
#' @param intervention A non-pharmaceutical intervention applied during the
#' course of the epidemic, with a start and end time, and age-specific effect
#' on contacts. See [intervention()].
#' @param vaccination A vaccination regime followed during the
#' course of the epidemic, with a group- and dose-specific start and end time,
#' and age-specific rates of delivery of first and second doses.
#' See [vaccination()].
#' @param increment The increment time, defaults to 0.1.
#' @param mortality_rate_vax The mortality rate of individuals who are
#' protected by vaccination.
#' @param hospitalisation_rate The hospitalisation rate of fully susceptible
#' and partially vaccinated and unprotected individuals.
#' @param hospitalisation_rate_vax The hospitalisation rate of individuals who
#' are protected by vaccination.
#' @param recovery_rate The recovery rate \eqn{\gamma}.
#' @param contact_matrix The population contact matrix.
#' @param npi_time_begin The start time of any non-pharmaceutical interventions
#' .
#' @param npi_time_end The end time of any non-pharmaceutical interventions.
#' @param npi_cr The reduction in contacts from any non-pharmaceutical
#' interventions.
#' @param vax_time_begin The start time of any vaccination campaigns.
#' @param vax_time_end The end time of any vaccination campaigns.
#' @param vax_nu The vaccination rate of any vaccination campaigns.
#' @param rate_interventions A named list of `<rate_intervention>` objects.
#' @param time_dependence A named list of functions for parameter time
#' dependence.
#' @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 (see [population()]).
#' The second list element is a vector of timesteps.
#' @keywords internal
.model_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_dependence, time_end = 100.0, increment = 1.0) {
.Call(`_epidemics_model_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_dependence, time_end, increment)
}

#' @title Compute the discrete probability of the truncated Erlang distribution
#' @name prob_discrete_erlang
#' @rdname prob_discrete_erlang
#'
prob_discrete_erlang_cpp <- function(shape, rate) {
.Call(`_epidemics_prob_discrete_erlang_cpp`, shape, rate)
.model_vacamole_cpp <- function(initial_state, transmissibility, transmissibility_vax, infectiousness_rate, mortality_rate, mortality_rate_vax, hospitalisation_rate, hospitalisation_rate_vax, recovery_rate, contact_matrix, npi_time_begin, npi_time_end, npi_cr, vax_time_begin, vax_time_end, vax_nu, rate_interventions, time_dependence, time_end = 100.0, increment = 1.0) {
.Call(`_epidemics_model_vacamole_cpp_internal`, initial_state, transmissibility, transmissibility_vax, infectiousness_rate, mortality_rate, mortality_rate_vax, hospitalisation_rate, hospitalisation_rate_vax, recovery_rate, contact_matrix, npi_time_begin, npi_time_end, npi_cr, vax_time_begin, vax_time_end, vax_nu, rate_interventions, time_dependence, time_end, increment)
}

54 changes: 24 additions & 30 deletions R/check_args_default.R
Original file line number Diff line number Diff line change
Expand Up @@ -8,14 +8,6 @@
#' `.check_args_model_default()` adds an empty `<intervention>` and
#' `<vaccination>` object if these are missing from the model arguments.
#'
#' `.prepare_args_epidemics_default()` prepares arguments for
#' [.model_default_cpp()], which is the C++ function that solves the default
#' ODE system using a Boost _odeint_ solver.
#' `.prepare_args_epidemics_default()` converts the arguments collected in
#' `mod_args` into simpler structures such as lists and numeric or integer
#' vectors that can be interpreted as C++ types such as `Rcpp::List`,
#' `Rcpp::NumericVector`, or `Eigen::MatrixXd`.
#'
#' @return
#'
#' `.check_args_model_default()` invisibly returns the model arguments passed
Expand All @@ -29,9 +21,9 @@
#' - `initial_state`: the initial conditions modified to represent absolute
#' rather than proportional values;
#'
#' - `beta`, `alpha`, `gamma`: three numbers representing the transmission rate
#' of the infection, the rate of transition from exposed to infectious, and the
#' recovery rate, respectively;
#' - `transmissibility`, `infectiousness_rate`, `recovery_rate`: three numbers
#' representing the transmission rate of the infection, the rate of transition
#' from exposed to infectious, and the recovery rate, respectively;
#'
#' - `contact_matrix`, a numeric matrix for the population contact matrix
#' scaled by the largest real eigenvalue and by the size of each groups;
Expand All @@ -51,19 +43,23 @@
#' is incremented.
#'
#' @keywords internal
#' @details
#' `.prepare_args_model_default()` prepares arguments for
#' [.model_default_cpp()], which is the C++ function that solves the default
#' ODE system using a Boost _odeint_ solver, and for [.ode_model_default()],
#' which is passed to [deSolve::lsoda()] in [model_default_r()].
#'
#' `.prepare_args_model_default()` converts the arguments collected in
#' `mod_args` into simpler structures such as lists and numeric or integer
#' vectors that can be interpreted as C++ types such as `Rcpp::List`,
#' `Rcpp::NumericVector`, or `Eigen::MatrixXd`.
.check_args_model_default <- function(mod_args) {
# check that arguments list has expected names
checkmate::assert_names(
names(mod_args),
type = "unique"
)

# input checking on the infection object
assert_infection(
mod_args[["infection"]],
extra_parameters = "preinfectious_period"
)

# load number of compartments to check initial conditions matrix
compartments_default <- c(
"susceptible", "exposed", "infectious", "recovered", "vaccinated"
Expand All @@ -89,7 +85,9 @@
# check for any other intervention list element names
checkmate::assert_names(
names(mod_args[["intervention"]]),
subset.of = c("beta", "gamma", "alpha", "contacts")
subset.of = c(
"transmissibility", "infectiousness_rate", "recovery_rate", "contacts"
)
)

# if a contacts intervention is passed, check it
Expand All @@ -107,18 +105,18 @@
}

# if there is only an intervention on contacts, add a dummy intervention
# on the transmission rate beta
# on the transmissibility
if (identical(names(mod_args[["intervention"]]), "contacts")) {
mod_args[["intervention"]]$beta <- no_rate_intervention()
mod_args[["intervention"]]$transmissibility <- no_rate_intervention()
}
} else {
# add as a list element named "contacts", and one named "beta"
# add as a list element named "contacts", and one named "transmissibility"
mod_args[["intervention"]] <- list(
contacts = no_contacts_intervention(
mod_args[["population"]]
),
# a dummy intervention on the rate parameter beta
beta = no_rate_intervention()
# a dummy intervention on the rate parameter transmissibility
transmissibility = no_rate_intervention()
)
}

Expand Down Expand Up @@ -163,12 +161,6 @@
get_parameter(mod_args[["population"]], "initial_conditions") *
get_parameter(mod_args[["population"]], "demography_vector")

# calculate infection parameters
gamma <- 1.0 / get_parameter(mod_args[["infection"]], "infectious_period")
alpha <- 1.0 / get_parameter(mod_args[["infection"]], "preinfectious_period")
beta <- get_parameter(mod_args[["infection"]], "r0") /
get_parameter(mod_args[["infection"]], "infectious_period")

# get NPI related times and contact reductions
contact_interventions <- mod_args[["intervention"]][["contacts"]]
npi_time_begin <- get_parameter(contact_interventions, "time_begin")
Expand All @@ -188,7 +180,9 @@
# return selected arguments for internal C++ function
list(
initial_state = initial_state,
beta = beta, alpha = alpha, gamma = gamma,
transmissibility = mod_args$transmissibility,
infectiousness_rate = mod_args$infectiousness_rate,
recovery_rate = mod_args$recovery_rate,
contact_matrix = contact_matrix,
npi_time_begin = npi_time_begin, npi_time_end = npi_time_end,
npi_cr = npi_cr,
Expand Down
Loading
Loading