| Title: | Support Code for the Introduction to Practical Disease Modelling Course |
|---|---|
| Description: | A collection of functions and classes that are used to illustrate fundamental concepts of disease modelling as part of a physical course series taught jointly by the University of Copenhagen and University of Sydney. The package is not intended to be used for other purposes; you are welcome to do so, but please note that the package may change without notice in order to suit changes to the associated taught course. All functions, arguments and classes are therefore subject to change without following any form of deprecation process. |
| Authors: | Matt Denwood [aut, cre], Carsten Kirkeby [rev], Michael Ward [rev] |
| Maintainer: | Matt Denwood <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 0.5.2-1 |
| Built: | 2025-02-07 17:26:02 UTC |
| Source: | https://github.com/ku-awdc/IPDMR |
A super-simple continuous-time two-compartment model
ab_continuous(rate = 0.1, time_points = seq(0, 10, by = 0.1))ab_continuous(rate = 0.1, time_points = seq(0, 10, by = 0.1))
rate |
the rate of progression from A to B per unit time (must be positive) |
time_points |
a vector of time points to include in the data frame returned |
A super-simple discrete-time two-compartment model
ab_discrete(rate = 0.1, d_time = 0.1, max_time = 10)ab_discrete(rate = 0.1, d_time = 0.1, max_time = 10)
rate |
the rate of progression from A to B per unit time (must be positive) |
d_time |
the desired time step (delta time) |
max_time |
the desired maximum time point (must be greater than the time step) |
This is utility function taking a vector of compartments and rates, and returning a matrix of calculated (for deterministic) or sampled (for stochastic) values for the change in compartment size from each input compartment (row) to each of the output compartments represented by the input rates (columns). When only a single rate is input then the function just returns the usual compartment * 1-exp(-rated_time) (for deterministic) or rbinom(length(compartment), compartment, 1-exp(-rated_time)) (for stochastic). Where there are multiple rates, then the function corrects for the overlap in the competing rates and (for stochastic) uses a multinomial distribution instead of a binomial distribution to ensure that the sum of each row of output is less than or equal to the corresponding input compartment size.
Multiple rows are only relevant where there are multiple sub-compartments, i.e. movement from the E state. In all other cases, a matrix with a single row is returned. Likewise, multiple columns are only relevant where there are multiple rates - in all other cases, a matrix with a single column is returned. If you know that you have a single compartment and/or a single rate, then you can either explicitly select the first row and/or column or use as.numeric() on the output to remove the dimensions attribute (see the bottom of the examples).
apply_rates( compartments, rates, d_time, update_type = c("deterministic", "stochastic") )apply_rates( compartments, rates, d_time, update_type = c("deterministic", "stochastic") )
compartments |
a vector of compartment sizes on which to apply all rates |
rates |
a vector of rates representing movement from each of the compartments |
d_time |
the size of the time step |
update_type |
deterministic or stochastic |
## A single compartment and rate: I <- 10; gamma <- 0.1; d_time <- 1 apply_rates(I, gamma, d_time, "deterministic") ### Is the same as: I * (1 - exp(-gamma*d_time)) ### Stochastic version uses rbinom: apply_rates(I, gamma, d_time, "stochastic") ## A single compartment and vector of rates: cull <- 0.01 apply_rates(I, c(gamma,cull), d_time, "deterministic") ## Due to competing rates, this is NOT the same as: I * (1 - exp(-gamma*d_time)) I * (1 - exp(-cull*d_time)) ### Stochastic version uses rmultinom: apply_rates(I, c(gamma,cull), d_time, "stochastic") ## A vector of compartments and single rate: E <- c(0,1,2); omega <- 0.1 apply_rates(E, omega, d_time, "deterministic") ### Is the same as: apply_rates(E[1], omega, d_time, "deterministic") apply_rates(E[2], omega, d_time, "deterministic") apply_rates(E[3], omega, d_time, "deterministic") ### Is the same as: E * (1 - exp(-omega*d_time)) ### Stochastic version uses rbinom: apply_rates(E, omega, d_time, "stochastic") ## A vector of compartments and vector of rates: repl <- 0.001 apply_rates(E, c(omega,repl), d_time, "deterministic") apply_rates(E, c(omega,repl), d_time, "stochastic") ## If the rate vector has names, then these are carried over to the ## output column names: apply_rates(I, c(toR=gamma,toS=cull), d_time, "deterministic") ## If you only have a single rate (and/or single compartment) then you ## probably want to remove the dimensions attribute before using, i.e.: apply_rates(E, omega, d_time, "deterministic")[,1] apply_rates(I, c(gamma,cull), d_time, "deterministic")[1,] ## Or: as.numeric(apply_rates(E, omega, d_time, "deterministic")) as.numeric(apply_rates(I, c(gamma,cull), d_time, "deterministic"))## A single compartment and rate: I <- 10; gamma <- 0.1; d_time <- 1 apply_rates(I, gamma, d_time, "deterministic") ### Is the same as: I * (1 - exp(-gamma*d_time)) ### Stochastic version uses rbinom: apply_rates(I, gamma, d_time, "stochastic") ## A single compartment and vector of rates: cull <- 0.01 apply_rates(I, c(gamma,cull), d_time, "deterministic") ## Due to competing rates, this is NOT the same as: I * (1 - exp(-gamma*d_time)) I * (1 - exp(-cull*d_time)) ### Stochastic version uses rmultinom: apply_rates(I, c(gamma,cull), d_time, "stochastic") ## A vector of compartments and single rate: E <- c(0,1,2); omega <- 0.1 apply_rates(E, omega, d_time, "deterministic") ### Is the same as: apply_rates(E[1], omega, d_time, "deterministic") apply_rates(E[2], omega, d_time, "deterministic") apply_rates(E[3], omega, d_time, "deterministic") ### Is the same as: E * (1 - exp(-omega*d_time)) ### Stochastic version uses rbinom: apply_rates(E, omega, d_time, "stochastic") ## A vector of compartments and vector of rates: repl <- 0.001 apply_rates(E, c(omega,repl), d_time, "deterministic") apply_rates(E, c(omega,repl), d_time, "stochastic") ## If the rate vector has names, then these are carried over to the ## output column names: apply_rates(I, c(toR=gamma,toS=cull), d_time, "deterministic") ## If you only have a single rate (and/or single compartment) then you ## probably want to remove the dimensions attribute before using, i.e.: apply_rates(E, omega, d_time, "deterministic")[,1] apply_rates(I, c(gamma,cull), d_time, "deterministic")[1,] ## Or: as.numeric(apply_rates(E, omega, d_time, "deterministic")) as.numeric(apply_rates(I, c(gamma,cull), d_time, "deterministic"))
Autoplot methods associated with IPDMR model functions
## S3 method for class 'ipdmr_ct' autoplot(object, ...)## S3 method for class 'ipdmr_ct' autoplot(object, ...)
object |
an object returned by an IPDMR model function |
... |
additional arguments (currently ignored) |
This is a general-purpose between-group model spread class that takes a list of within-group spread models as input, and provides methods to update/run (and save/reset) these models while taking care of between-group (density-dependent) spread by calculating and inputting trans_external for each of the groups as appropriate based on the beta_matrix, which can be modified as needed (see examples).
The input models can either be created using the make_group function, or they can be your own R6 (or C++) class that implement the necessary methods and fields (i.e. reset, save, update, I, state, and trans_external). The models need not all be the same subtype, i.e. you can mix SIR/SEIR/SI models and even R6 and C++ implementations.
Note that the between-group model stores a reference to (i.e. shallow copy of) the input models, so you can see the state of the input groups being modified as the between-group model runs. You can also change the state of these groups at any time between updates/runs, either by modifying the groups used to create the between-group model class or by using the (read-only) active binding groups to retrieve one or more reference to a group. Of course you can also change the beta_matrix at any time between updates/runs.
beta_matrixa matrix representing between-group (density-dependent) spread
groupsa list of the internally-stored within-group models (this can be used to directly check and/or modify their state and/or parameters)
timethe current time point
transmission_betweenthe between-group transmission type (frequency or density)
statea data frame of the current state of each group
new()
Create a new between-group model
BetweenGroupClass$new(groups)
groupsa list of within-group models as either R6 or C++ classes (see e.g. make_group)
A new between-group model object
update()
Update the state of each group for a single time point, including setting trans_external based on the value of the beta_matrix and vector of I obtained from the groups (density-dependent spread).
BetweenGroupClass$update(d_time)
d_timethe desired time step (delta time)
self, invisibly
save()
Instruct each group to save the current state and parameter values for later retrieval using reset()
BetweenGroupClass$save()
self, invisibly
reset()
Instruct each group to reset the current state and parameter values to their last saved state
BetweenGroupClass$reset()
self, invisibly
run()
Update the state of each group for a given number of time points
BetweenGroupClass$run(add_time, d_time)
add_timethe additional time to add to the current time of the model
d_timethe desired time step (delta time)
a data frame of the model state at each (new) time point
run_until()
Update the state of each group for several time points, until a stopping criterion is met. NOTE: this method is a stub: it has not yet been implemented!
BetweenGroupClass$run_until(criterion_fun, d_time)
criterion_funa function taking an input data frame of states and returning a logical scalar indicating if the simulation should be stopped (TRUE) or continue to be updated (FALSE)
d_timethe desired time step (delta time)
a data frame of the model state at each (new) time point
print()
Print method showing the number of groups and current time
BetweenGroupClass$print()
self, invisibly
clone()
The objects of this class are cloneable with this method.
BetweenGroupClass$clone(deep = FALSE)
deepWhether to make a deep clone.
## Create a list of groups: groups <- list( make_group("stochastic", group_name="Group1"), make_group("stochastic", group_name="Group2") ) ## Modify the starting conditions as needed: groups[[1]]$S <- 100 groups[[1]]$I <- 0 groups[[2]]$S <- 99 groups[[2]]$I <- 1 ## etc ## Create a between-group model: model <- BetweenGroupClass$new(groups) ## Modify beta_matrix as needed: model$beta_matrix <- matrix(c(0,0.1,0.01,0), nrow=2, ncol=2) ## Then run the model: output <- model$run(21, d_time=0.1) ggplot2::autoplot(output) ## Change the beta_matrix to e.g. remove transmission from Group1->Group2: model$beta_matrix <- matrix(c(0,0,0.01,0), nrow=2, ncol=2) ## And run for more time points: output2 <- model$run(21, d_time=0.1) ggplot2::autoplot(output2) ## To combine time points: all_output <- dplyr::bind_rows(output, output2) summary(all_output)## Create a list of groups: groups <- list( make_group("stochastic", group_name="Group1"), make_group("stochastic", group_name="Group2") ) ## Modify the starting conditions as needed: groups[[1]]$S <- 100 groups[[1]]$I <- 0 groups[[2]]$S <- 99 groups[[2]]$I <- 1 ## etc ## Create a between-group model: model <- BetweenGroupClass$new(groups) ## Modify beta_matrix as needed: model$beta_matrix <- matrix(c(0,0.1,0.01,0), nrow=2, ncol=2) ## Then run the model: output <- model$run(21, d_time=0.1) ggplot2::autoplot(output) ## Change the beta_matrix to e.g. remove transmission from Group1->Group2: model$beta_matrix <- matrix(c(0,0,0.01,0), nrow=2, ncol=2) ## And run for more time points: output2 <- model$run(21, d_time=0.1) ggplot2::autoplot(output2) ## To combine time points: all_output <- dplyr::bind_rows(output, output2) summary(all_output)
These wrapper functions create a single-group model, run the model for the indicated time steps, and then return a data frame representing the model output. Note that starting values for S/E/I/R can be continuous or discrete for the deterministic models, but must be discrete numbers for the stochastic models.
sir_det( S = 99, I = 1, R = 0, beta = 0.25, gamma = 0.2, delta = 0.05, transmission_type = c("frequency", "density"), d_time = 1, max_time = 100 ) seir_det( S = 99, E = 0, I = 1, R = 0, numE = 3L, beta = 0.25, omega = 0.2, gamma = 0.2, delta = 0.05, vacc = 0, repl = 0, cull = 0, transmission_type = c("frequency", "density"), d_time = 1, max_time = 100 ) sir_stoc( S = 99, I = 1, R = 0, beta = 0.25, gamma = 0.2, delta = 0.05, transmission_type = c("frequency", "density"), d_time = 1, max_time = 100, iterations = 1L ) seir_stoc( S = 99, E = 0, I = 1, R = 0, numE = 3L, beta = 0.25, omega = 0.2, gamma = 0.2, delta = 0.05, vacc = 0, repl = 0, cull = 0, transmission_type = c("frequency", "density"), d_time = 1, max_time = 100, iterations = 1L )sir_det( S = 99, I = 1, R = 0, beta = 0.25, gamma = 0.2, delta = 0.05, transmission_type = c("frequency", "density"), d_time = 1, max_time = 100 ) seir_det( S = 99, E = 0, I = 1, R = 0, numE = 3L, beta = 0.25, omega = 0.2, gamma = 0.2, delta = 0.05, vacc = 0, repl = 0, cull = 0, transmission_type = c("frequency", "density"), d_time = 1, max_time = 100 ) sir_stoc( S = 99, I = 1, R = 0, beta = 0.25, gamma = 0.2, delta = 0.05, transmission_type = c("frequency", "density"), d_time = 1, max_time = 100, iterations = 1L ) seir_stoc( S = 99, E = 0, I = 1, R = 0, numE = 3L, beta = 0.25, omega = 0.2, gamma = 0.2, delta = 0.05, vacc = 0, repl = 0, cull = 0, transmission_type = c("frequency", "density"), d_time = 1, max_time = 100, iterations = 1L )
S |
the starting number of susceptible animals |
I |
the starting number of infected animals |
R |
the starting number of recovered animals |
beta |
the transmission rate parameter per unit time (must be positive) |
gamma |
the recovery rate parameter per unit time (must be positive) |
delta |
the reversion rate parameter per unit time (must be positive) |
transmission_type |
either frequency or density |
d_time |
the desired time step (delta time) |
max_time |
the desired maximum time point (must be greater than the time step) |
E |
the total starting number of exposed/latent animals (this will be distributed among the sub-compartments of E when the model is initialised) |
numE |
the number of sub-compartments within the exposed/latent state |
omega |
the latent progression rate parameter per unit time (must be positive) |
vacc |
the vaccination rate parameter per unit time (must be positive) |
repl |
the replacement rate parameter per unit time (must be positive) |
cull |
the targeted culling rate parameter per unit time (must be positive) |
iterations |
the number of iterations to run (stochastic models only) |
This is a wrapper function around the new/initialise methods for R6 and C++ implementations of within-group spread models used on the IPDMR course. For more information on available methods see the documentation for SEIRclass (the method signatures and outputs of the R6 and C++ implementations should be identical).
make_group( update_type = c("deterministic", "stochastic"), numE = 3L, numI = 1L, numR = 1L, group_name = NA_character_, implementation = c("C++", "R6") )make_group( update_type = c("deterministic", "stochastic"), numE = 3L, numI = 1L, numR = 1L, group_name = NA_character_, implementation = c("C++", "R6") )
update_type |
either stochastic or deterministic |
numE |
the number of sub-compartments desired for the E state (0 or more) |
numI |
the number of sub-compartments desired for the I state (1 or more) |
numR |
the number of sub-compartments desired for the R state (0 or more) |
group_name |
an optional name for the group |
implementation |
either C++ or R6 |
set.seed(2024-11-05) r6res <- make_group(update_type = "stochastic", implementation = "R6")$run(10, 0.1) set.seed(2024-11-05) cppres <- make_group(update_type = "stochastic", implementation = "C++")$run(10, 0.1) ## Should be identical: stopifnot(unlist(r6res) - unlist(cppres) == 0L)set.seed(2024-11-05) r6res <- make_group(update_type = "stochastic", implementation = "R6")$run(10, 0.1) set.seed(2024-11-05) cppres <- make_group(update_type = "stochastic", implementation = "C++")$run(10, 0.1) ## Should be identical: stopifnot(unlist(r6res) - unlist(cppres) == 0L)
These wrapper functions create a multi-group model, run the model for the indicated time steps, and then return a data frame representing the model output. Note that starting values for S/E/I/R can be continuous or discrete for the deterministic models, but must be discrete numbers for the stochastic models.
multi_seir_det( n_groups, beta_matrix, S = 99, E = 0, I = 1, R = 0, numE = 3L, beta = 0.25, omega = 0.2, gamma = 0.2, delta = 0.05, vacc = 0, repl = 0, cull = 0, transmission_within = c("frequency", "density"), transmission_between = c("density", "frequency"), d_time = 1, max_time = 100 ) multi_seir_stoc( n_groups, beta_matrix, S = 99, E = 0, I = 1, R = 0, numE = 3L, beta = 0.25, omega = 0.2, gamma = 0.2, delta = 0.05, vacc = 0, repl = 0, cull = 0, transmission_within = c("frequency", "density"), transmission_between = c("density", "frequency"), d_time = 1, max_time = 100, iterations = 1L )multi_seir_det( n_groups, beta_matrix, S = 99, E = 0, I = 1, R = 0, numE = 3L, beta = 0.25, omega = 0.2, gamma = 0.2, delta = 0.05, vacc = 0, repl = 0, cull = 0, transmission_within = c("frequency", "density"), transmission_between = c("density", "frequency"), d_time = 1, max_time = 100 ) multi_seir_stoc( n_groups, beta_matrix, S = 99, E = 0, I = 1, R = 0, numE = 3L, beta = 0.25, omega = 0.2, gamma = 0.2, delta = 0.05, vacc = 0, repl = 0, cull = 0, transmission_within = c("frequency", "density"), transmission_between = c("density", "frequency"), d_time = 1, max_time = 100, iterations = 1L )
n_groups |
the number of groups to generate |
beta_matrix |
a matrix representing between-group spread (must have nrow and ncol equal to n_groups) |
S |
the starting number of susceptible animals |
E |
the total starting number of exposed/latent animals (this will be distributed among the sub-compartments of E when the model is initialised) |
I |
the starting number of infected animals |
R |
the starting number of recovered animals |
numE |
the number of sub-compartments within the exposed/latent state |
beta |
the transmission rate parameter per unit time (must be positive) |
omega |
the latent progression rate parameter per unit time (must be positive) |
gamma |
the recovery rate parameter per unit time (must be positive) |
delta |
the reversion rate parameter per unit time (must be positive) |
vacc |
the vaccination rate parameter per unit time (must be positive) |
repl |
the replacement rate parameter per unit time (must be positive) |
cull |
the targeted culling rate parameter per unit time (must be positive) |
transmission_within |
the within-group transmission type (frequency or density) |
transmission_between |
the between-group transmission type (frequency or density) |
d_time |
the desired time step (delta time) |
max_time |
the desired maximum time point (must be greater than the time step) |
iterations |
the number of iterations to run (stochastic models only) |
output <- multi_seir_det(2, matrix(0,nrow=2,ncol=2)) ggplot2::autoplot(output) output <- multi_seir_stoc(2, matrix(0,nrow=2,ncol=2), iterations=3L) ggplot2::autoplot(output)output <- multi_seir_det(2, matrix(0,nrow=2,ncol=2)) ggplot2::autoplot(output) output <- multi_seir_stoc(2, matrix(0,nrow=2,ncol=2), iterations=3L) ggplot2::autoplot(output)
This is a within-group SEIR model class that can either be run on its own or embedded in a between-group model. Both R6 and C++ implementations are available - see make_group for a wrapper function that returns either implementation. However, the method signatures and outputs of the R6 and C++ implementations are (or at least should be!) identical, so this documentation applies to the C++ implementation as well as the R6 implementation.
Snumber of susceptible animals
Etotal number of exposed individuals - changing this value will cause the provided total to be distributed evenly (deterministic) or randomly (stochastic) between sub-compartments
Inumber of infected animals
Rnumber of recovered animals
betathe transmission rate parameter per unit time (must be positive)
omegathe latent progression rate parameter per unit time (must be positive)
gammathe recovery rate parameter per unit time (must be positive)
deltathe reversion rate parameter per unit time (must be positive)
vaccthe vaccination rate parameter per unit time (must be positive)
replthe replacement rate parameter per unit time (must be positive)
cullthe targeted culling rate parameter per unit time (must be positive)
timethe current time point of the model (read-only)
Nthe total number of animals in the group (read-only)
statea data frame representing the current state of the model (read-only)
trans_externalthe external transmission parameter
transmission_typeeither frequency or density
new()
Create a new within-group model (this would normally be called from make_group)
SEIRclass$new(
update_type = c("deterministic", "stochastic"),
numE = 3L,
group_name = NA_character_
)update_typeone of deterministic or stochastic
numEthe number of sub-compartments within the exposed/latent state
group_namean optional name for the group (will be included in the output state, if provided)
A new within-group model object
save()
Save the current state and parameter values for later retrieval using reset()
SEIRclass$save()
self, invisibly
reset()
Reset the current state and parameter values to their last saved state
SEIRclass$reset()
self, invisibly
update()
Update the state of the group for a single time point
SEIRclass$update(d_time)
d_timethe desired time step (delta time)
self, invisibly
run()
Update the state of the group for several time points
SEIRclass$run(add_time, d_time)
add_timethe additional time to add to the current time of the model
d_timethe desired time step (delta time)
a data frame of the model state at each (new) time point
print()
Print method giving an overview of the current state and parameter values
SEIRclass$print()
self, invisibly
clone()
The objects of this class are cloneable with this method.
SEIRclass$clone(deep = FALSE)
deepWhether to make a deep clone.
A simple continuous-time SI model
si_continuous( S = 9, I = 1, beta = 0.05, transmission_type = c("frequency", "density"), time_points = seq(0, 21, by = 0.1) )si_continuous( S = 9, I = 1, beta = 0.05, transmission_type = c("frequency", "density"), time_points = seq(0, 21, by = 0.1) )
S |
the starting number of susceptible animals (continuous number) |
I |
the starting number of infected animals (continuous number) |
beta |
the transmission rate per unit time (must be positive) |
transmission_type |
either frequency or density |
time_points |
a vector of time points to include in the data frame returned |
A simple discrete-time SI model
si_discrete( S = 9, I = 1, beta = 0.05, transmission_type = c("frequency", "density"), d_time = 1/24, max_time = 21 )si_discrete( S = 9, I = 1, beta = 0.05, transmission_type = c("frequency", "density"), d_time = 1/24, max_time = 21 )
S |
the starting number of susceptible animals (continuous number) |
I |
the starting number of infected animals (continuous number) |
beta |
the transmission rate per unit time (must be positive) |
transmission_type |
either frequency or density |
d_time |
the desired time step (delta time) |
max_time |
the desired maximum time point (must be greater than the time step) |
si_discrete(S=9, I=1, transmission_type="density") |> ggplot2::autoplot() si_discrete(S=9, I=1, transmission_type="frequency") |> ggplot2::autoplot()si_discrete(S=9, I=1, transmission_type="density") |> ggplot2::autoplot() si_discrete(S=9, I=1, transmission_type="frequency") |> ggplot2::autoplot()
This is a utility function that takes a desired mean and standard deviation for the distribution of transit times through a disease compartment, and suggests a suitable rate and number of sub-compartments that matches the provided mean and variance as closely as possible.
suggest_boxes(mean, sd)suggest_boxes(mean, sd)
mean |
the desired mean waiting time |
sd |
the desired standard deviation in waiting time |
the function returns a ggplot object illustrating the suggestions, and prints the suggested number of sub-compartments and rates to screen
# An average of 5 days, with sd of 2.5 days: suggest_boxes(mean=5, sd=3) # The number of sub-compartments is invariant to transformation of the rate # on different time scales, i.e. we get the same if specifying in hours: suggest_boxes(mean=5*24, sd=3*24)# An average of 5 days, with sd of 2.5 days: suggest_boxes(mean=5, sd=3) # The number of sub-compartments is invariant to transformation of the rate # on different time scales, i.e. we get the same if specifying in hours: suggest_boxes(mean=5*24, sd=3*24)