vignettes/L-How_To_Add_Model.Rmd
L-How_To_Add_Model.RmdThis vignette discusses how to add a new L module to
ramp.library.
As background, we strongly recommend reading the
ramp.xds vignette about the L
Interface.
Each L 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
L module.
Do this:
Download the markdown file that created this document “X-How_To_Add_Lodel.Rmd.”
Set your working directory to the location of the downloaded copy of “X-How_To_Add_Lodel.Rmd”
Execute these commands:
library(knitr)
knitr::purl("L-How_To_Add_Model.Rmd",
output = "Lmodule-newLname.R",
documentation = 0,
quiet=TRUE)You should see a new file called Lmodule-newLname.R
in the working directory. You can now move and rename it.
Open the file and do a global search for newLname
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 L_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 L with the your first variable
(e.g. mosquito population density), and replace
F_L with the formula that computes the derivatives. Do the
same with Y and add lines for each variable.
Lodify derivs = c(dL, 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 [dLdt] for the newLname model
#' @inheritParams ramp.xds::dLdt
#' @return a [numeric] vector
#' @export
dLdt.newLname <- function(t, y, xds_obj, i) {
# do not change this
eta <- xds_obj$terms$eta[[i]]
# do not change this
# attach the variables by name
with(get_L_vars(y, xds_obj, i),{
# do not change this
# attach the parameters by name
with(xds_obj$L_obj[[i]], {
# Do not change this
dL <- F_L(eta, L)
...
# Change this
derivs = c(dL, ...)
# return the derivatives
return(derivs)
})
})
}
#' @title Compute the derivatives for parasite infection dynamics in human population strata
#' @description Implements [UpdateLt] for the newLname model
#' @inheritParams ramp.xds::UpdateLt
#' @return a [numeric] vector
#' @export
Update_Lt.newLname<- function(t, y, xds_obj, s) {
# do not change this
eta <- xds_obj$terms$eta[[i]]
# do not change this
# attach the variables by name
with(get_L_vars(y, xds_obj, i),{
# do not change this
# attach the parameters by name
with(xds_obj$L_obj[[i]], {
# Change this
Lt <- L + Lambda + F_L(L)
...
# Change this
states = c(Lt, ...)
# return the derivatives
return(states)
})
})
}Two functions update bionomic parameter values at each time step, called in sequence before the dynamics are computed.
LBionomics.newLname computes the current bionomic
parameter values by calling the forcing function for each parameter. It
also resets all effect sizes to 1.
Replace p1, p2, … with the names of the
bionomic parameters in your model. Add a call to a forcing function for
each one.
For each parameter p1, add a matching effect size
reset: es_p1 <- rep(1, nHabitats).
#' @title Mosquito bionomics for `newLname` (**L**)
#' @description Implements [LBionomics] for the newLname model
#' @inheritParams ramp.xds::LBionomics
#' @return an **`xds`** object
#' @export
LBionomics.newLname <- function(t, y, xds_obj, s) {
with(xds_obj$L_obj[[s]], {
# Change this: call the forcing function for each bionomic parameter
xds_obj$L_obj[[s]]$p1_t <- F_p1(t, xds_obj, s)
xds_obj$L_obj[[s]]$p2_t <- F_p2(t, xds_obj, s)
# Change this: reset one effect size for each bionomic parameter
xds_obj$L_obj[[s]]$es_p1 <- rep(1, nHabitats)
xds_obj$L_obj[[s]]$es_p2 <- rep(1, nHabitats)
return(xds_obj)
})}LEffectSizes.newLname applies vector control effect
sizes to the bionomic parameters set by LBionomics. For
each parameter, the effective value is the baseline times the effect
size.
p1 with baseline p1_t
and effect size es_p1, add a line:
xds_obj$L_obj[[s]]$p1 <- p1_t * es_p1.
#' @title Apply effect sizes for `newLname` (**L**)
#' @description Implements [LEffectSizes] for the newLname model
#' @inheritParams ramp.xds::LEffectSizes
#' @return an **`xds`** object
#' @export
LEffectSizes.newLname <- function(t, y, xds_obj, s) {
with(xds_obj$L_obj[[s]], {
# Change this: multiply each baseline parameter by its effect size
xds_obj$L_obj[[s]]$p1 <- p1_t * es_p1
xds_obj$L_obj[[s]]$p2 <- p2_t * es_p2
return(xds_obj)
})}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_L.newLname. If you do want your module to have
that functionality, add the effects of mass treatment to the parameters
that define state transitions in dLdt or
Update_Lt and set mda=TRUE and
msat=TRUE in the skill set.
#' @title Lake parameters for newLname 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_L_obj_newLname = function(nPatches, options=list(),
p1=1, p2=2, p3=3){
with(options,{
L_obj = list()
class(L_obj) <- c("newLname")
# Change this
L_obj$p1 = checkIt(p1, nPatches)
L_obj$p2 = checkIt(p2, nPatches)
L_obj$p3 = checkIt(p3, nPatches)
# Don't change this
L_obj <- setup_K_obj(nPatches, L_obj)
return(L_obj)
})}If you did a global search and replace on newLname, this
function will not need to be modified
#' @title Setup L_obj.newLname
#' @description Implements [setup_L_obj] for the newLname model
#' @inheritParams ramp.xds::setup_L_obj
#' @return a [list] vector
#' @export
setup_L_obj.newLname = function(newLname, xds_obj, i, options=list()){
xds_obj$L_obj[[i]] = make_L_obj_newLname(xds_obj$nPatches, options)
return(xds_obj)
}In the following, each variable in your model gets its own pair of
lines. Replace L in L_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_L_ix] for the newLname model.
#' @inheritParams ramp.xds::setup_L_ix
#' @return none
#' @importFrom utils tail
#' @export
setup_L_ix.newLname <- function(xds_obj, i) {with(xds_obj,{
L_ix <- seq(from = max_ix+1, length.out=xds_obj$nPatches[i])
max_ix <- tail(L_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$L_obj[[i]]$ix = list(L_ix=L_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 L in L and L_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$L_obj`
#' @inheritParams ramp.xds::get_L_vars
#' @return a [list]
#' @export
get_L_vars.newLname <- function(y, xds_obj, i) {
with(xds_obj$L_obj[[i]]$ix,{
H = y[H_ix]
L = y[L_ix]
Y = y[Y_ix]
...
X0 = H-L-Y
return(list(H=H, L=L, Y=Y,...,X0=X0))})
}In the following, each variable in your model gets its own pair of lines.
Leave H as-is
Replace L in L and L_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 newLname model
#' @description Implements [parse_L_orbits] for the newLname model
#' @inheritParams ramp.xds::parse_L_orbits
#' @return a named list
#' @export
parse_L_orbits.newLname <- function(outputs, xds_obj, i) {
with(xds_obj$L_obj[[i]]$ix,{
H <- outputs[,H_ix]
L <- outputs[,L_ix]
Y <- outputs[,Y_ix]
X0 <- H-L-Y
vars <- list(H=H, L=L, 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 *newLname*
#'
#' @inheritParams ramp.xds::setup_L_inits
#'
#' @return a **`ramp.xds`** object
#'
#' @export
setup_L_inits.newLname = function(xds_obj, H, i, options=list()){
xds_obj$L_obj[[i]]$inits = make_L_inits_newLname(xds_obj$nPatches[i], H, options)
return(xds_obj)
}
#' @title Lake 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 L0 the initial value for L
#' @param Y0 the initial value for L
#' @return a [list]
#' @export
make_L_inits_newLname = function(nPatches, H, options = list(), L0=NULL, Y0=1){with(options,{
stopifnot(is.numeric(L0))
stopifnot(is.numeric(Y0))
L = checkIt(L0, nPatches)
Y = checkIt(Y0, nPatches)
return(list(L=L, Y=Y, ...))
})}
#' @title Return the parameters as a list
#'
#' @inheritParams ramp.xds::change_L_inits
#'
#' @return an **`xds`** object
#' @export
change_L_inits.SIS <- function(xds_obj, i=1, options=list()) {
with(get_L_inits(xds_obj, i), with(options,{
xds_obj = change_H(H, xds_obj, i)
xds_obj$Xinits[[i]]$L = L
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_emerge] for the SIS model.
#' @inheritParams ramp.xds::F_emerge
#' @return a [numeric] vector of length `nPatches`
#' @export
F_emerge.newLname <- function(y, xds_obj, i) {
with(get_L_vars(y, xds_obj, i),
with(xds_obj$L_obj[[i]], {
emerge = ...
return(emerge)
}))}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_L] for the SIS model.
#' @inheritParams ramp.xds::F_L
#' @return a [numeric] vector of length `nPatches`
#' @export
F_L.newLname <- function(y, xds_obj, i) {
with(get_L_vars(y, xds_obj, i),
with(xds_obj$L_obj[[i]], {
L = ...
return(L)
}))}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 L 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(L, nPatches, clrs=c("darkblue","darkred"), llty=1){
with(L,{
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 newLname model
#'
#' @param L a list with the outputs of parse_deout_X_newLname
#' @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_newLname = function(L, nPatches, clrs=c("darkblue","darkred"), llty=1){
with(L,{
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){
UseLethod("xds_plot_X", xds_obj$L_obj[[i]])
}
#' Plot the density of infected individuals for the newLname 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$L[[i]],
plot(time, 0*time, type = "n", ylim = c(0, max(H)),
ylab = "# Infected", xlab = "Time"))
xds_lines_X_SIS(vars$L[[i]], xds_obj$Hpar[[i]]$nPatches, clrs, llty)
}
#' Plot the density of infected individuals for the newLname model
#'
#' @inheritParams ramp.xds::xds_plot_X
#' @export
xds_plot_X.newLname = 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$L[[i]],
plot(time, 0*time, type = "n", ylim = c(0, max(H)),
ylab = "# Infected", xlab = "Time"))
xds_lines_X_newLname(vars$L[[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 Lalaria Transmission 2
This vignette was written to help add new modules to
ramp.library.
Reiner RC Jr, Perkins TA, Barker CL, Niu T, Chaves LF, Ellis AL, 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 JL, Citron DT, Lbabazi 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↩︎