vignettes/MY-How_To_Add_Model.Rmd
MY-How_To_Add_Model.RmdThis vignette discusses how to add a new MY module
to ramp.library.
As background, we strongly recommend reading the
ramp.xds vignette about the MY
Interface.
Each MY module implements a unique model family
describing human/host infection dynamics and immunity that can be used
by ramp.xds.
Before writing a new module, please check and see if the model family is
already part of the SimBA library.
If you need help, please contact Professor David L Smith
We have written this vignette to produce a new template file using
knitr::purl. The template can be modified to build a new
MY module.
Do this:
Download the markdown file that created this document “X-How_To_Add_Model.Rmd.”
Set your working directory to the location of the downloaded copy of “X-How_To_Add_Model.Rmd”
Execute these commands:
library(knitr)
knitr::purl("X-How_To_Add_Model.Rmd",
output = "MYmodule-newMYname.R",
documentation = 0,
quiet=TRUE)You should see a new file called
MYmodule-newMYname.R in the working directory. You can now
move and rename it.
Open the file and do a global search for newMYname
and replace it with the a short string that will define the model as its
Xname.
Edit the functions following the directions below.
We recommend that new models conform to a set of style guidelines:
All the names of the variables by name using the indices attached
to MY_obj
All documentation is written in roxygen2
Functions should use with(*,{...}) to make the code
easier to read.
A vignette should accompany each model, including citations to the relevant publications, and proposed tests of the code.
The functions should easy to relate to the equations defining a model, following notation defined in an accompanying vignette.
The sections that need to be modified are marked
# Change this
Replace M with the your first variable
(e.g. mosquito population density), and replace
F_M with the formula that computes the derivatives. Do the
same with Y and add lines for each variable.
Modify derivs = c(dM, dY, ...) to return the LHS of
the derivative equations you just wrote.
#' @title Compute the derivatives for parasite infection dynamics in human population strata
#' @description Implements [dMYdt] for the newMYname model
#' @inheritParams ramp.xds::dMYdt
#' @return a [numeric] vector
#' @export
dMYdt.newMYname <- function(t, y, xds_obj, i) {
# do not change this
Lambda <- xds_obj$terms$Lambda[[i]]
kappa <- xds_obj$terms$kappa[[i]]
# do not change this
# attach the variables by name
with(get_MY_vars(y, xds_obj, i),{
# do not change this
# attach the parameters by name
with(xds_obj$MY_obj[[i]], {
# Do not change this
dM <- F_M(Lambda, M)
dY <- F_Y(kappa, M, Y)
...
# Change this
derivs = c(dM, dY, ...)
# return the derivatives
return(derivs)
})
})
}
#' @title Compute the derivatives for parasite infection dynamics in human population strata
#' @description Implements [UpdateMYt] for the newMYname model
#' @inheritParams ramp.xds::UpdateMYt
#' @return a [numeric] vector
#' @export
Update_MYt.newMYname<- function(t, y, xds_obj, s) {
Lambda <- xds_obj$terms$Lambda[[i]]
kappa <- xds_obj$terms$kappa[[i]]
# do not change this
# attach the variables by name
with(get_MY_vars(y, xds_obj, i),{
# do not change this
# attach the parameters by name
with(xds_obj$MY_obj[[i]], {
# Change this
Mt <- M + Lambda + F_M(M)
Yt <- Y + F_Y(kappa, M, Y)
...
# Change this
states = c(Mt, Yt, ...)
# return the derivatives
return(states)
})
})
}In the following, replace p1 with the first
parameter, p2 with the second, and so on.
Unless you know what you’re doing, leave the “Ports for human / host demography” unchanged.
If you don’t want to add ports for mass treatment, delete the
line with mda and msat and leave
mda=FALSE and msat=FALSE in
skill_set_MY.newMYname. If you do want your module to have
that functionality, add the effects of mass treatment to the parameters
that define state transitions in dMYdt or
Update_MYt and set mda=TRUE and
msat=TRUE in the skill set.
#' @title Make parameters for newMYname human model, with defaults
#' @param nPatches is the number of population strata
#' @param options a [list] that could overwrite defaults
#' @param p1 the first parameter
#' @param p2 the second parameter
#' @param p3 the third parameter
#' @return a [list]
#' @export
make_MY_obj_newMYname = function(nPatches, options=list(),
p1=1, p2=2, p3=3){
with(options,{
MY_obj = list()
class(MY_obj) <- c("newMYname")
# Change this
MY_obj$p1 = checkIt(p1, nPatches)
MY_obj$p2 = checkIt(p2, nPatches)
MY_obj$p3 = checkIt(p3, nPatches)
# Don't change this
MY_obj <- setup_K_obj(nPatches, MY_obj)
return(MY_obj)
})}If you did a global search and replace on newMYname,
this function will not need to be modified
#' @title Setup MY_obj.newMYname
#' @description Implements [setup_MY_obj] for the newMYname model
#' @inheritParams ramp.xds::setup_MY_obj
#' @return a [list] vector
#' @export
setup_MY_obj.newMYname = function(newMYname, xds_obj, i, options=list()){
xds_obj$MY_obj[[i]] = make_MY_obj_newMYname(xds_obj$nPatches, options)
return(xds_obj)
}In the following, each variable in your model gets its own pair of
lines. Replace M in M_ix with the first
variable name, Y in Y_ix with the second, and
add lines for each variable by copying and pasting and editing the pair
of lines. Be sure the replacement gets applied to all four
instances.
#' @title Add indices for human population to parameter list
#' @description Implements [setup_MY_ix] for the newMYname model.
#' @inheritParams ramp.xds::setup_MY_ix
#' @return none
#' @importFrom utils tail
#' @export
setup_MY_ix.newMYname <- function(xds_obj, i) {with(xds_obj,{
M_ix <- seq(from = max_ix+1, length.out=xds_obj$nPatches[i])
max_ix <- tail(M_ix, 1)
Y_ix <- seq(from = max_ix+1, length.out=xds_obj$nPatches[i])
max_ix <- tail(Y_ix, 1)
...
xds_obj$max_ix = max_ix
xds_obj$MY_obj[[i]]$ix = list(M_ix=M_ix, Y_ix=Y_ix, ...)
return(xds_obj)
})}In the following, each variable in your model gets its own pair of lines.
Leave H as-is
Replace M in M and M_ix
with the first variable name.
Replace Y in Y and Y_ix
with the second variable name.
Copy the form for every other variable
The variable X0 is the class that gets newborns. If
you use it or if you want to have it available, compute it and add it to
the returned list.
#' @title Return the variables as a list
#' @description This method dispatches on the type of `xds_obj$MY_obj`
#' @inheritParams ramp.xds::get_MY_vars
#' @return a [list]
#' @export
get_MY_vars.newMYname <- function(y, xds_obj, i) {
with(xds_obj$MY_obj[[i]]$ix,{
H = y[H_ix]
M = y[M_ix]
Y = y[Y_ix]
...
X0 = H-M-Y
return(list(H=H, M=M, Y=Y,...,X0=X0))})
}In the following, each variable in your model gets its own pair of lines.
Leave H as-is
Replace M in M and M_ix
with the first variable name.
Replace Y in Y and Y_ix
with the second variable name.
Copy the form for every other variable
The variable X0 is the class that gets newborns. If
you use it or if you want to have it available, compute it and add it to
the returned list.
#' @title parse the output of deSolve and return variables for the newMYname model
#' @description Implements [parse_MY_orbits] for the newMYname model
#' @inheritParams ramp.xds::parse_MY_orbits
#' @return a named list
#' @export
parse_MY_orbits.newMYname <- function(outputs, xds_obj, i) {
with(xds_obj$MY_obj[[i]]$ix,{
H <- outputs[,H_ix]
M <- outputs[,M_ix]
Y <- outputs[,Y_ix]
X0 <- H-M-Y
vars <- list(H=H, M=M, Y=Y, X0=X0)
return(vars)
})}If you did the global search and replace, you won’t need to change this function.
#' @title Setup initial values for *newMYname*
#'
#' @inheritParams ramp.xds::setup_MY_inits
#'
#' @return a **`ramp.xds`** object
#'
#' @export
setup_MY_inits.newMYname = function(xds_obj, H, i, options=list()){
xds_obj$MY_obj[[i]]$inits = make_MY_inits_newMYname(xds_obj$nPatches[i], H, options)
return(xds_obj)
}
#' @title Make initial values for the Xname human model, with defaults
#' @param nPatches the number of strata in the model
#' @param options a [list] to overwrite defaults
#' @param M0 the initial value for M
#' @param Y0 the initial value for M
#' @return a [list]
#' @export
make_MY_inits_newMYname = function(nPatches, H, options = list(), M0=NULL, Y0=1){with(options,{
stopifnot(is.numeric(M0))
stopifnot(is.numeric(Y0))
M = checkIt(M0, nPatches)
Y = checkIt(Y0, nPatches)
return(list(M=M, Y=Y, ...))
})}
#' @title Return the parameters as a list
#'
#' @inheritParams ramp.xds::change_MY_inits
#'
#' @return an **`xds`** object
#' @export
change_MY_inits.SIS <- function(xds_obj, i=1, options=list()) {
with(get_MY_inits(xds_obj, i), with(options,{
xds_obj = change_H(H, xds_obj, i)
xds_obj$Xinits[[i]]$M = M
xds_obj$Xinits[[i]]$Y = Y
...
return(xds_obj)
}))}In the following, replace … with a formula to compute the effective density of infectious humans: it is population density by strata weighted by the probability of transmitting.
#' @title Size of effective infectious human population
#' @description Implements [F_fqZ] for the SIS model.
#' @inheritParams ramp.xds::F_fqZ
#' @return a [numeric] vector of length `nPatches`
#' @export
F_fqZ.newMYname <- function(y, xds_obj, i) {
with(get_MY_vars(y, xds_obj, i),
with(xds_obj$MY_obj[[i]], {
fqZ = ...
return(fqZ)
}))}In the following, replace … with a formula to compute the effective density of infectious humans: it is population density by strata weighted by the probability of transmitting.
#' @title Size of effective infectious human population
#' @description Implements [F_fqM] for the SIS model.
#' @inheritParams ramp.xds::F_fqM
#' @return a [numeric] vector of length `nPatches`
#' @export
F_fqM.newMYname <- function(y, xds_obj, i) {
with(get_MY_vars(y, xds_obj, i),
with(xds_obj$MY_obj[[i]], {
fqM = ...
return(fqM)
}))}
#' @title Infection blocking pre-erythrocytic immunity
#' @description Implements [F_eggs] for the newMYname model.
#' @inheritParams ramp.xds::F_eggs
#' @return a [numeric] vector of length `nPatches`
#' @export
F_eggs.newMYname <- function(y, xds_obj, i) {
with(xds_obj$MY_obj[[i]],{
########################
# retrieve or compute it
########################
eggs = ...
########################
# return it
########################
return(eggs)
})
}To make ramp.xds easy to use, we provide some functions
that parse the outputs and compute some standard metrics.
#' Add lines for the density of infected individuals for the SIS model
#'
#' @param MY a list with the outputs of parse_deout_X_SIS
#' @param nPatches the number of population strata
#' @param clrs a vector of colors
#' @param llty an integer (or integers) to set the `lty` for plotting
#'
#' @export
xds_lines_X_SIS = function(MY, nPatches, clrs=c("darkblue","darkred"), llty=1){
with(MY,{
if(nPatches==1) {
lines(time, S, col=clrs[1], lty = llty[1])
lines(time, I, col=clrs[2], lty = llty[1])
}
if(nPatches>1){
if (length(clrs)==2) clrs=matrix(clrs, 2, nPatches)
if (length(llty)==1) llty=rep(llty, nPatches)
for(i in 1:nPatches){
lines(time, S[,i], col=clrs[1,i], lty = llty[i])
lines(time, I[,i], col=clrs[2,i], lty = llty[i])
}
}
})}
#' Add lines for the density of infected individuals for the newMYname model
#'
#' @param MY a list with the outputs of parse_deout_X_newMYname
#' @param nPatches the number of population strata
#' @param clrs a vector of colors
#' @param llty an integer (or integers) to set the `lty` for plotting
#'
#' @export
xds_lines_X_newMYname = function(MY, nPatches, clrs=c("darkblue","darkred"), llty=1){
with(MY,{
if(nPatches==1) {
lines(time, S, col=clrs[1], lty = llty[1])
lines(time, I, col=clrs[2], lty = llty[1])
}
if(nPatches>1){
if (length(clrs)==2) clrs=matrix(clrs, 2, nPatches)
if (length(llty)==1) llty=rep(llty, nPatches)
for(i in 1:nPatches){
lines(time, S[,i], col=clrs[1,i], lty = llty[i])
lines(time, I[,i], col=clrs[2,i], lty = llty[i])
}
}
})}We provide a basic plotting function. The function
xds_lines_X_* does not have an S3 version.
Instead, it is called by the S3 function
xds_plot_X
#' Basic plotting for epidemiological models
#'
#' @param xds_obj a list that defines an `ramp.xds` model (*e.g.*, generated by `xde_setup()`)
#' @param i the host species index
#' @param clrs a vector of colors
#' @param llty an integer (or integers) to set the `lty` for plotting
#' @param stable a logical: set to FALSE for `orbits` and TRUE for `stable_orbits`
#' @param add_axes a logical: plot axes only if TRUE
#'
#' @export
xds_plot_X = function(xds_obj, i=1, clrs="black", llty=1, stable=FALSE, add_axes=TRUE){
UseMethod("xds_plot_X", xds_obj$MY_obj[[i]])
}
#' Plot the density of infected individuals for the newMYname model
#'
#' @inheritParams ramp.xds::xds_plot_X
#' @export
xds_plot_X.SIS = function(xds_obj, i=1, clrs=c("darkblue","darkred"), llty=1, stable=FALSE, add_axes=TRUE){
vars=with(xds_obj$outputs,if(stable==TRUE){stable_orbits}else{orbits})
if(add_axes==TRUE)
with(vars$MY[[i]],
plot(time, 0*time, type = "n", ylim = c(0, max(H)),
ylab = "# Infected", xlab = "Time"))
xds_lines_X_SIS(vars$MY[[i]], xds_obj$Hpar[[i]]$nPatches, clrs, llty)
}
#' Plot the density of infected individuals for the newMYname model
#'
#' @inheritParams ramp.xds::xds_plot_X
#' @export
xds_plot_X.newMYname = function(xds_obj, i=1, clrs=c("darkblue","darkred"), llty=1, stable=FALSE, add_axes=TRUE){
vars=with(xds_obj$outputs,if(stable==TRUE){stable_orbits}else{orbits})
if(add_axes==TRUE)
with(vars$MY[[i]],
plot(time, 0*time, type = "n", ylim = c(0, max(H)),
ylab = "# Infected", xlab = "Time"))
xds_lines_X_newMYname(vars$MY[[i]], xds_obj$Hpar[[i]]$nPatches, clrs, llty)
}ramp.xds has a modular design. The
modular design makes it possible to break down published models to serve
as dynamical components in new models for malaria, a kind of
plug-and-play functionality for two kinds of systems:
systems of ordinary and delay differential equations models for malaria and other mosquito-transmitted pathogens.
stochastic and deterministic discrete time systems for malaria and other mosquito-transmitted pathogens.
This package – ramp.library – is a collection of models
to implement a large number of dynamical model families and other
algorithms taken from the literature describing malaria and other
mosquito-transmitted pathogens (see Reiner, et al. 2013)1. Those
models have been implemented here as code that has undergone rigorous
testing and that can be used to build models for simulation-based
analytics.
The mathematics supporting this framework were explained in Spatial Dynamics of Malaria Transmission 2
This vignette was written to help add new modules to
ramp.library.
Reiner RC Jr, Perkins TA, Barker CM, Niu T, Chaves LF, Ellis AM, et al. A systematic review of mathematical models of mosquito-borne pathogen transmission: 1970-2010. J R Soc Interface. 2013;10: 20120921.↩︎
Wu SL, Henry JM, Citron DT, Mbabazi Ssebuliba D, Nakakawa Nsumba J, Sánchez C. HM, et al. (2023) Spatial dynamics of malaria transmission. PLoS Comput Biol 19(6): e1010684. https://doi.org/10.1371/journal.pcbi.1010684↩︎