knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
library(EpiModel)
This vignette discusses mechanisms usable inside EpiModel
network models with custom extension modules. More information about these may be found in the "New Network Models with EpiModel" section of the EpiModel Tutorials.
In a model, parameters are the input variables used to define aspects of the system behavior. In the basic built-in SIS (Susceptible-Infected-Susceptible) model, these parameters could be the act rate, the infection probability and the recovery rate. In simple models, each of these parameters are single fixed values that do not change over the course of a simulation. In more complex models, we may want more flexibility in model parameterization.
Therefore, in this vignette, we demonstrate how to implement:
First, we define a simple SI model.
nw <- network_initialize(n = 50) est <- netest( nw, formation = ~edges, target.stats = 25, coef.diss = dissolution_coefs(~offset(edges), 10, 0), verbose = FALSE ) param <- param.net( inf.prob = 0.3, act.rate = 0.5, dummy.param = 4, dummy.strat.param = c(0, 1) ) init <- init.net(i.num = 10) control <- control.net(type = "SI", nsims = 1, nsteps = 5, verbose = FALSE) mod <- netsim(est, param, init, control) mod
In the parameters, we set the value for inf.prob
and act.rate
as fixed, but we also define dummy.param
and dummy.strat.param
. These latter two parameters will serve to illustrate random parameters. dummy.strat.param
has two elements; this may be used, for example, to stratify a parameter by subpopulation.
The last line prints a summary of the model. In it we can see the value of the parameters under the "Fixed Parameters" section. Note the additional groups
parameter defined automatically by EpiModel
as part of the "SI" model definition.
To allow our parameters to be drawn from a distribution of random values, we use the random.params
argument to param.net
. There are two ways of defining which parameters are random, and the distribution of values for those random parameters to draw randomly and how to draw them.
The first option is to define a generator function for each parameter we want to treat as random.
my.randoms <- list( act.rate = param_random(c(0.25, 0.5, 0.75)), dummy.param = function() rbeta(1, 1, 2), dummy.strat.param = function() c( rnorm(1, 0.05, 0.01), rnorm(1, 0.15, 0.03) ) ) param <- param.net( inf.prob = 0.3, random.params = my.randoms ) param
Here we kept the inf.prob
parameter fixed at 0.3
and defined a list
object called my.randoms
containing 3 elements:
act.rate
uses the param_random
function factory defined by EpiModel
(see ?EpiModel::param_random
).dummy.param
is a function with no argument that returns a random value from a beta distribution.dummy.strat.param
is a function with no argument that returns 2 values sampled from normal distributions, each with different means and standard deviations.Each element is named after the parameter it will fill and MUST BE a function taking no argument and outputting a vector of the right size for the parameter: size 1 for act.rate
and dummy.param
; size 2 for dummy.strat.param
.
The rest of the model is run as before, although we increase the simulation count to three to demonstrate the parameter stochasticity.
control <- control.net(type = "SI", nsims = 3, nsteps = 5, verbose = FALSE) mod <- netsim(est, param, init, control) mod
After running 3 simulations we can see that 2 parameters are still displayed under the "Fixed Parameters" section: inf.prob
and groups
. The other ones are displayed under the "Random Parameters". act.rate
and dummy.param
now have 3 values associated with them, one per simulation. dummy.strat.param
have <list>
as value because each simulation has a vector of size 2.
We can inspect the values with the get_param_set
function:
``` {r generators_inspect} get_param_set(mod)
#### Parameter set The drawback of generator function approach above is that it cannot produce correlated parameters. For instance, we may want `dummy.param` and `dummy.strat.param` to be related to one another within a single simulation. We might also use another, external approach for generating parameter sets (e.g., Latin hypercube sampling of multiple parameters). For this, we need to pre-define a `data.frame` of the possible values: ```r n <- 5 related.param <- data.frame( dummy.param = rbeta(n, 1, 2) ) related.param$dummy.strat.param_1 <- related.param$dummy.param + rnorm(n) related.param$dummy.strat.param_2 <- related.param$dummy.param * 2 + rnorm(n) related.param
We now have a data.frame
with 5 rows and 3 columns. Each row contains a set of parameters values that will be used together in a model. This way we keep the relationship between each value.
The first column of the data.frame
is named dummy.param
, similar to the name of the parameter. For dummy.start.param
we need two columns as the parameter contains two values. To achieve this, we name the 2 columns dummy.start.param_1
and
dummy.strat.param_2
. The value after the underscore informs EpiModel
how it should combine these values. This in turn means that the underscore symbol is not allowed in proper parameter names.
Then we set up the rest of the parameters. related.param
is saved in the my.randoms
list under the special reserved name param_random_set
.
my.randoms <- list( act.rate = param_random(c(0.25, 0.5, 0.75)), param_random_set = related.param ) param <- param.net( inf.prob = 0.3, random.params = my.randoms ) param
Notice here that we combined a generator function for act.rate
and a set of
correlated parameters with param_random_set
.
control <- control.net(type = "SI", nsims = 3, nsteps = 5, verbose = FALSE) mod <- netsim(est, param, init, control) mod
The output is similar to what we saw with the generator functions.
When doing modeling work, the researcher may want to simulate interventions or events that vary over time. In (our 2021 paper)[https://academic.oup.com/jid/article/223/6/1019/6122459], we simulated the effects of COVID-19–Related behavior change and clinical service interruption on HIV and STI incidence. These included simulations with systemic events with a beginning and an end. To work with these situations, EpiModel offers an optional module to programmatically handle parameter changes across time steps.
This vignette assumes a basic understanding of EpiModel custom models. If the code does not make sense, please check the EpiModel Tutorials.
To define what parameters should change and when during the simulation, we need to define updaters. An updater is a list
with to named elements: at
,
the time step when the change will take place, and param
a named list of parameters to update.
list( at = 10, param = list( inf.prob = 0.3, act.rate = 0.5 ) )
This example defines an updater that will change the value of the inf.prob
parameter to 0.3
and the value to the act.rate
parameter to 0.5
. This
change will happen during the 10th time step.
As mentioned above, we usually want to define several changes at once. EpiModel accept a list
of updaters.
# Create a `list.of.updaters` list.of.updaters <- list( # this is one updater list( at = 100, param = list( inf.prob = 0.3, act.rate = 0.3 ) ), # this is another updater list( at = 125, param = list( inf.prob = 0.01 ) ) ) # The `list.of.updaters` goes into `param.net` under `param.updater.list` param <- param.net( inf.prob = 0.1, act.rate = 0.1, param.updater.list = list.of.updaters )
In this example, we define 2 updaters, one that occurs at time step 100 and the other one at time step 20. In practice, inf.prob
and act.rate
are set to 0.1
at the beginning by param.net
, at time step 100 they are both updated to 0.3
and at time step 125 act.rate
is reduced to 0.01
.
updater
ModuleThis functionality is part of an optional module provided by EpiModel
, updater.net
. This module must then be enabled in control.net
.
Below we set up a complete example with a closed population SI model using the parameters and updaters defined above. We then plot the size of the infected and susceptible population over time to see the effects of the updaters.
control <- control.net( type = NULL, # must be NULL as we use a custom module nsims = 1, nsteps = 200, verbose = FALSE, updater.FUN = updater.net, infection.FUN = infection.net ) nw <- network_initialize(n = 50) nw <- set_vertex_attribute(nw, "race", rbinom(50, 1, 0.5)) est <- netest( nw, formation = ~edges, target.stats = 25, coef.diss = dissolution_coefs(~offset(edges), 10, 0), verbose = FALSE ) init <- init.net(i.num = 10) mod <- netsim(est, param, init, control) plot(mod)
This is the minimal example to simplify the demonstration of these time-varying parameters.
Important: when enabling the updater.net
module, the researcher must pay particular attention to the order the modules are run. It is recommended to have
the updater.net
module run first so the changes happen at the beginning of the time step.
Some functionality of updater.net
were not described above and are explained in ?EpiModel::updater.net
.
When creating an updater, one can add an optional verbose
element to the list
. If TRUE
, updater.net
will output a message
describing what changes
where performed when it occurs during the simulation.
list( at = 10, param = list( inf.prob = 0.3, act.rate = 0.5 ), verbose = TRUE )
It is sometime useful to configure the changes with respect to the current value instead of a fixed new value. This is possible, as demonstrated below.
list( at = 10, param = list( inf.prob = function(x) plogis(qlogis(x) + log(2)), act.rate = 0.5 ) )
This updater will set the value of act.rate
to 0.5
as we saw earlier. But, for inf.prob
we put a function
(not a function call). In this case,
updater.net
will apply the function
to the current value of act.rate
. If we consider as in the previous example that act.rate
is set to 0.1
by
param.net
, then its new value will be obtained by adding an Odds Ratio of 2 to the original value plogis(qlogis(0.1) + log(2))
: 0.1818182
.
Similarly to time-varying parameters we can use time-varying controls. They work in the same way and use the same updater
module. The two differences are:
control
sub-list instead of a param
sub-listcontrol$control.updater.list
# Create a `list.of.updaters` list.of.updaters <- list( # this is one updater list( at = 100, control = list( resimulate.network = FALSE ) ), # this is another updater list( at = 125, control = list( verbose = FALSE ) ) ) # The `list.of.updaters` goes into `control.net` under `control.updater.list` control <- control.net( type = NULL, # must be NULL as we use a custom module nsims = 1, nsteps = 200, verbose = TRUE, resimulate.network = TRUE, updater.FUN = updater.net, infection.FUN = infection.net, control.updater.list = list.of.updaters )
This example set 2 updaters, one that will turn off network resimulation at timestep 100 and another toggling of the verbosity at step 125.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.