knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
library(EpiModel)
This vignette assumes familiarity with setting up and running network models in EpiModel. See the Network Modeling for Epidemics (NME) course materials and the EpiModel Gallery for background.
For working with nodal attributes and epidemic summary statistics, see the companion vignette Working with Custom Attributes and Summary Statistics in EpiModel. For network objects and edgelists, see Working with Network Objects in EpiModel.
In a model, parameters are the input variables that define aspects of system behavior. In a basic SIS (Susceptible-Infected-Susceptible) model, these might include an infection probability, an act rate, and a recovery rate. In simple models, each parameter is a single fixed value that remains constant throughout a simulation. Real-world applications, however, often require more flexibility: parameters that change over time (e.g., to represent an intervention rollout or behavioral change during an epidemic), parameters drawn from distributions (for sensitivity analysis or calibration), or large parameter sets managed via spreadsheets.
This vignette demonstrates how to implement:
data.frame, enabling spreadsheet-based parameter management.Scenarios allow you to define named sets of parameter values that can change at specific time steps during a simulation. This is the recommended approach for implementing time-varying parameters---for example, modeling an intervention that reduces transmission probability partway through an epidemic, or comparing multiple policy alternatives.
Our research paper implementing time-varying parameter scenarios for COVID-19 sexual distancing interventions was published in the Journal of Infectious Diseases. The code can be found here.
First, we set up a simple SIS model as the base case.
set.seed(10) nw <- network_initialize(n = 200) est <- netest(nw, formation = ~edges, target.stats = 60, coef.diss = dissolution_coefs(~offset(edges), 10, 0), # duration = 10, closed population verbose = FALSE ) param <- param.net(inf.prob = 0.9, rec.rate = 0.01, act.rate = 2) init <- init.net(i.num = 10) control <- control.net(type = "SIS", nsims = 1, nsteps = 250, verbose = FALSE)
Scenarios are defined using create_scenario_list, which takes a specially formatted data.frame as input and outputs a list of scenarios usable by EpiModel. We use tibble::tribble here for readability, but any data.frame constructor works.
library(dplyr) scenarios.df <- tribble( ~.scenario.id, ~.at, ~inf.prob, ~rec.rate, "base", 0, 0.9, 0.01, "initial_change", 0, 0.2, 0.01, "multiple_changes", 0, 0.1, 0.04, "multiple_changes", 100, 0.9, 0.01, "multiple_changes", 200, 0.1, 0.1 ) knitr::kable(scenarios.df)
The data.frame requires two columns:
.scenario.id: an identifier for the scenario..at: the time step at which the parameter changes should apply. Use 0 for changes that take effect before the simulation begins.Remaining columns are parameter names. They must start with a letter and contain only letters, numbers, or . (underscores are reserved for vector indexing, described below). If a cell is NA, EpiModel sets that parameter to NA---it does not mean "leave unchanged."
In this example, the "multiple_changes" scenario has three rows sharing the same .scenario.id, representing a single scenario with parameter changes at three different time points.
To convert from the data.frame to a usable list of scenarios:
scenarios.list <- create_scenario_list(scenarios.df) str(scenarios.list, max.level = 2)
We loop over all scenarios, using use_scenario to create a modified parameter object for each. Under the hood, use_scenario generates parameter updaters (described in the Advanced section below) that apply the specified changes at the designated time steps.
# List to hold simulation results d_list <- vector(mode = "list", length = length(scenarios.list)) names(d_list) <- names(scenarios.list) for (scenario in scenarios.list) { print(scenario$id) sc.param <- use_scenario(param, scenario) sim <- netsim(est, sc.param, init, control) d_sim <- as.data.frame(sim) d_sim[["scenario"]] <- scenario$id d_list[[scenario$id]] <- d_sim }
For the "multiple_changes" scenario, messages appear at time steps 100 and 200 indicating that parameters were modified during the simulation. The other scenarios are silent because their changes occur before the simulation begins (.at = 0).
Now we can plot the results. We use base R plotting here rather than plot.netsim() because we need to overlay results from separate netsim runs on a single axis.
plot(d_list$base$time, d_list$base$i.num, type = "l", col = 1, lwd = 2, ylim = c(0, 250), xlab = "Time Step", ylab = "Number Infected") lines(d_list$initial_change$time, d_list$initial_change$i.num, type = "l", col = 2, lwd = 2) lines(d_list$multiple_changes$time, d_list$multiple_changes$i.num, type = "l", col = 3, lwd = 2) abline(v = c(100, 200), lty = 2) legend("topleft", legend = names(d_list), col = 1:3, lwd = 2, cex = 0.9, bty = "n")
The "initial_change" scenario shows a slower epidemic than the base case because the infection probability is reduced from 0.9 to 0.2 from the start.
The "multiple_changes" scenario demonstrates three phases:
In full-scale modeling projects, parameters are often vectors. For example, hiv.test.rate might be a vector of length 3 containing weekly HIV screening probabilities for Black, Hispanic, and White persons.
To use vectors in scenarios, create a separate column for each element using the naming convention paramname_N, where N is the position in the vector:
scenarios.df <- tribble( ~.scenario.id, ~.at, ~hiv.test.rate_1, ~hiv.test.rate_2, ~hiv.test.rate_3, "base", 0, 0.001, 0.001, 0.001, "initial_change", 0, 0.002, 0.001, 0.002, "multiple_changes", 0, 0.002, 0.001, 0.002, "multiple_changes", 100, 0.004, 0.001, 0.004, "multiple_changes", 200, 0.008, 0.001, 0.008 ) knitr::kable(scenarios.df)
All elements of a vector parameter must be specified in the data.frame, even if some values do not change across scenarios (e.g., hiv.test.rate_2 is constant here but must still be included).
When working with many parameters, we recommend storing the scenarios.df as a CSV, RDS or Excel file for easier sharing and editing.
# CSV write.csv(scenarios.df, "scenarios.csv", row.names = FALSE) scenarios.df <- read.csv("scenarios.csv") # RDS saveRDS(scenarios.df, "scenarios.rds") scenarios.df <- readRDS("scenarios.rds")
The param.net function accepts a data.frame of parameters through the data.frame.params argument. This allows passing many parameters at once or working with a spreadsheet to track and update parameters before use in EpiModel.
The data.frame requires three columns:
param: The parameter name. For vector parameters (length > 1), append the position with an underscore (e.g., "p_1", "p_2").value: The parameter value (as a character string).type: The type of the parameter value. Only three values are accepted: "numeric", "logical", or "character".Additional columns (e.g., details, source) may be included for documentation but are ignored by EpiModel.
df_params <- tribble( ~param, ~value, ~type, "hiv.test.rate_1", "0.003", "numeric", "hiv.test.rate_2", "0.102", "numeric", "hiv.test.rate_3", "0.492", "numeric", "prep.require.lnt", "TRUE", "logical", "group_1", "first", "character", "group_2", "second", "character" ) knitr::kable(df_params)
Pass the table via data.frame.params. These parameters can be combined with named parameters for maximum flexibility. In case of conflict, named parameters take priority over those in the data.frame:
param <- param.net(data.frame.params = df_params, other.param = c(5, 10), act.rate = 1) param
Fixed parameters assume that each input value is known with certainty. In practice, we may want to explore how uncertainty in parameter values propagates to uncertainty in epidemic outcomes. EpiModel supports drawing parameter values from distributions, so that each simulation uses a different realization.
We demonstrate with 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), # duration = 10, closed population 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
Here we define four parameters: inf.prob (which will remain fixed), act.rate, dummy.param, and dummy.strat.param (which we will make random below). The dummy.strat.param parameter is a vector of length 2, which could represent a parameter stratified by subpopulation.
To draw parameters from distributions, use the random.params argument to param.net. There are two approaches.
Define a generator function for each random parameter:
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
The my.randoms list contains three elements:
act.rate uses the param_random function factory provided by EpiModel (see ?param_random). Each simulation samples one of the three values with equal probability.dummy.param is a zero-argument function that returns a single draw from a Beta(1, 2) distribution.dummy.strat.param is a zero-argument function that returns a vector of length 2, each element drawn from a normal distribution with different mean and standard deviation.Each element must be named after the parameter it fills and must be a function taking no arguments, returning a vector of the correct length for that parameter. When we print the parameter list before running the model, random parameters appear as function definitions since their values have not yet been realized.
control <- control.net(type = "SI", nsims = 3, nsteps = 5, verbose = FALSE) mod <- netsim(est, param, init, control) mod
After running 3 simulations, inf.prob remains under "Fixed Parameters" while the random parameters each have 3 realized values (one per simulation). The vector parameter dummy.strat.param shows <list> because each realization is itself a vector of length 2.
To inspect the realized values:
``` {r generators_inspect} all.params <- get_param_set(mod) all.params
These can be merged with epidemic output for analysis of parameter-outcome relationships: ```r epi <- as.data.frame(mod) left_join(epi, all.params)
Generator functions draw each parameter independently. When parameters need to be correlated---for example, when using Latin hypercube sampling or when one parameter is derived from another---use pre-defined parameter sets instead.
Define a data.frame where each row is a complete set of correlated parameter values:
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
Each row contains parameter values that will be used together in a single simulation. Vector parameters use the same paramname_N suffix convention as scenarios. This means underscores are reserved and cannot appear in parameter names themselves.
Save the parameter set in the my.randoms list under the 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
The inf.prob parameter remains fixed, act.rate remains independently random, and the two remaining parameters are drawn as correlated sets from the data.frame.
control <- control.net(type = "SI", nsims = 3, nsteps = 5, verbose = FALSE) mod <- netsim(est, param, init, control) mod
Verify that correlated sets are sampled together:
related.param
get_param_set(mod)
The scenarios system described above uses an internal updater module to implement parameter changes at designated time steps. This section describes the underlying updater mechanism directly. Use this when the scenario API is not flexible enough---for example, when you need relative (function-based) parameter changes or time-varying control settings.
An updater is a list with two required elements: at (the time step for the change) and param (a named list of new parameter values):
list( at = 10, param = list( inf.prob = 0.3, act.rate = 0.5 ) )
This updater sets inf.prob to 0.3 and act.rate to 0.5 at time step 10. Multiple updaters are combined in a list:
list.of.updaters <- list( list( at = 100, param = list( inf.prob = 0.3, act.rate = 0.3 ) ), list( at = 125, param = list( inf.prob = 0.01 ) ) )
Pass the updater list to param.net via the .param.updater.list argument:
param <- param.net( inf.prob = 0.1, act.rate = 0.1, .param.updater.list = list.of.updaters ) init <- init.net(i.num = 10) control <- control.net( type = "SI", nsims = 1, nsteps = 200, verbose = FALSE )
nw <- network_initialize(n = 100) est <- netest( nw, formation = ~edges, target.stats = 50, coef.diss = dissolution_coefs(~offset(edges), 10, 0), # duration = 10, closed population verbose = FALSE ) mod <- netsim(est, param, init, control)
The plot shows inflection points at the two parameter change points. At step 100, both inf.prob and act.rate increase from 0.1 to 0.3, accelerating transmission. At step 125, inf.prob drops to 0.01, slowing the epidemic.
plot(mod, mean.smooth = FALSE) abline(v = c(100, 125), lty = 2, col = "grey50")
Each updater can include an optional verbose element. When TRUE (the default), EpiModel prints a message describing the changes:
list( at = 10, param = list( inf.prob = 0.3, act.rate = 0.5 ), verbose = TRUE )
Instead of setting a parameter to a fixed new value, you can define a function that transforms the current value. This is useful for multiplicative changes or logit-scale adjustments.
list( at = 10, param = list( inf.prob = function(x) plogis(qlogis(x) + log(2)), act.rate = 0.5 ) )
In this updater, act.rate is set to the fixed value 0.5 as before. For inf.prob, instead of a value, we provide a function. At time step 10, EpiModel applies this function to the current value of inf.prob. If inf.prob was set to 0.1 by param.net, the function computes plogis(qlogis(0.1) + log(2)) = r round(plogis(qlogis(0.1) + log(2)), 4), which doubles the odds of transmission while keeping the result on the probability scale.
Control settings can also be changed during a simulation using the same updater mechanism. Each control updater uses a control element (instead of param):
list.of.updaters <- list( list( at = 100, control = list( resimulate.network = FALSE ) ), list( at = 125, control = list( verbose = FALSE ) ) )
This example turns off network resimulation at step 100 and disables model verbosity at step 125.
Pass control updaters to control.net via .control.updater.list, paralleling how parameter updaters are passed to param.net via .param.updater.list:
control <- control.net( type = "SI", nsims = 1, nsteps = 200, verbose = TRUE, .control.updater.list = list.of.updaters )
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.