Package 'EpiModel'

Description

Tools for simulating mathematical models of infectious disease dynamics. Epidemic model classes include deterministic compartmental models, stochastic individual-contact models, and stochastic network models. Network models use the robust statistical methods of exponential-family random graph models (ERGMs) from the Statnet suite of software packages in R. Standard templates for epidemic modeling include SI, SIR, and SIS disease types. EpiModel features an API for extending these templates to address novel scientific research aims. Full methods for EpiModel are detailed in Jenness et al. (2018, doi:10.18637/jss.v084.i08).

Details

EpiModel provides tools for building, simulating, and analyzing mathematical models of infectious disease dynamics. Models may be built for pedagogical purposes using the built-in disease types, or for research purposes using the extension API to define arbitrary disease states and processes.

The distinguishing capability of EpiModel is stochastic network modeling, in which partnerships form, persist, and dissolve according to a statistical model estimated from network data. Simpler compartmental and individual-based classes are also provided, both for teaching and for comparison against the network results.

Model Classes and Infectious Disease Types

EpiModel provides functionality for three classes of epidemic models:

• Deterministic Compartmental Models (DCMs): continuous-time models solved with ordinary differential equations, under an assumption of random mixing within the population. Sensitivity analyses across parameter values or initial conditions are specified by passing vectors to param.dcm() or init.dcm().

• Stochastic Individual Contact Models (ICMs): discrete-time, individual-based analogs of the DCM class, which add random variation to each component of the transmission system, from infection to recovery to vital dynamics (arrivals and departures). Contacts are randomly mixed, with no persistent partnership structure.

• Stochastic Network Models: discrete-time, individual-based models in which the contact structure is represented explicitly as a dynamic network. Partnership formation and dissolution are simulated from a temporal exponential-family random graph model (TERGM), using the statistical framework implemented in the Statnet suite of R packages. Disease then spreads across the edges of that evolving network. This class represents features that random mixing cannot, including degree distributions, partnership duration, concurrency, assortative mixing, and feedback between network structure and epidemic dynamics.

EpiModel supports three built-in infectious disease types, available across all three model classes, in one-group and two-group configurations and with or without demography:

• Susceptible-Infectious (SI): a two-state disease with life-long infection and no recovery. HIV is one example, although for that case it is common to model infection stages as separate compartments.

• Susceptible-Infectious-Recovered (SIR): a three-state disease with life-long immunity following infection. Measles is one example, although modern models for measles also require consideration of vaccination patterns in the population.

• Susceptible-Infectious-Susceptible (SIS): a two-state disease in which one may transition back and forth between the susceptible and infected states throughout life. Examples include bacterial sexually transmitted infections such as gonorrhea.

These built-in types are starting points. The extension API described below allows them to be replaced with arbitrarily complex disease systems.

Model Parameterization and Simulation

Each model class uses three setup functions to collect the epidemic parameters, initial conditions, and control settings:

• param.dcm(), param.icm(), and param.net() set the epidemic parameters, such as the rate of contacts or acts between agents, the probability of transmission per act, and the recovery and demographic rates for models that include those transitions.

• init.dcm(), init.icm(), and init.net() set the initial conditions, which for built-in types are the numbers or, where applicable, the specific agents infected or recovered at the outset.

• control.dcm(), control.icm(), and control.net() set the remaining simulation controls, including the disease type, the number of time steps, and the number of simulations or runs.

With the model parameterized, the simulation functions are:

• dcm() for deterministic compartmental models.

• icm() for stochastic individual contact models.

• Network models are simulated in a three-step process:

  1. netest() estimates the statistical model for the network structure itself, that is, how partnerships form and dissolve over time. This function wraps the estimation routines in the ergm and tergm packages. The typical input is a set of target statistics summarizing the formation process together with a mean partnership duration, an approach known as egocentric inference because those targets may be calculated from an egocentric sample of the population. Dissolution coefficients are calculated with dissolution_coefs().

  2. netdx() runs diagnostics on the fitted model by simulating the network forward in time without disease, to check that the simulated network recovers the formation and dissolution targets. The net-controls page covers the estimation and simulation settings to adjust when it does not.

  3. netsim() simulates the stochastic network epidemic model, given the fitted network model from netest() along with the parameters, initial conditions, and controls defined above.

Network models may also use multiple network layers over a shared node set, for example separate sexual and needle-sharing contact networks. Each layer is estimated separately and the fits are passed together to netsim(), with multilayer() specifying any controls that vary by layer.

Extending the Built-In Models

Research applications typically require disease states or processes beyond the built-in types. The extension API supports this for the network class: setting type = NULL in control.net() and supplying module functions (for example, infection.FUN = my_infection_module) replaces or supplements the built-in modules. A module is a function with the signature ⁠function(dat, at)⁠ that reads from and returns the netsim_dat object holding all simulation state.

Modules should use the accessor functions rather than manipulating that object directly. These include get_attr() and set_attr() for nodal attributes, get_param() for parameters, get_epi() and set_epi() for epidemic summary statistics, append_core_attr() and append_attr() for new arrivals, discord_edgelist() for serodiscordant pairs, and set_transmat() for recording transmission events.

Extension support is specific to the network class. ICMs support the built-in SI, SIR, and SIS types only; users needing custom individual-based dynamics should use netsim(). DCMs accept an original system of differential equations through the new.mod argument to control.dcm().

Working with Simulation Output

All three model classes return objects with print, summary, plot, and as.data.frame methods, so the same analysis workflow applies throughout. Simulation output may be converted to a long data frame with as.data.frame.dcm(), as.data.frame.icm(), or as.data.frame.netsim(), extended with derived statistics using mutate_epi(), subset with get_sims(), and combined across runs with merge.icm() or merge.netsim().

Network models additionally record the simulated networks, retrieved with get_network(); the network statistics, retrieved with get_nwstats(); and the transmission tree, retrieved with get_transmat() and convertible to a phylogeny with as.phylo.transmat(). Partnership histories may be tracked over time with the cumulative edgelist and queried with get_partners(), get_forward_reachable(), and get_backward_reachable().

Learning Resources

The package vignettes cover getting started (Intro) and then working with model parameters, network objects, and nodal attributes and summary statistics in extension models.

Beyond the package documentation:

• The Network Modeling for Epidemics course materials provide a full treatment of network model theory, estimation, diagnostics, and epidemic simulation, followed by the extension API and advanced topics such as multi-layer networks and epidemics over fully observed networks, at https://epimodel.github.io/sismid/.

• The EpiModel Gallery provides worked extension models, from adding a single compartment through multi-stage disease models with interventions and cost-effectiveness analysis, at https://epimodel.github.io/EpiModel-Gallery/.

• The EpiModel website collects these resources and a suggested learning pathway, at https://www.epimodel.org/.

Author(s)

Maintainer: Samuel Jenness samuel.m.jenness@emory.edu

Authors:

• Samuel Jenness samuel.m.jenness@emory.edu

• Steven M. Goodreau goodreau@uw.edu

• Martina Morris morrism@uw.edu

• Adrien Le Guillou contact@aleguillou.org

• Chad Klumb cklumb@gmail.com

Other contributors:

• Skye Bender-deMoll skyebend@uw.edu [contributor]

References

The EpiModel website is at https://www.epimodel.org/, and the source code is at https://github.com/EpiModel/EpiModel. Bug reports and feature requests are welcome.

Our primary methods paper on EpiModel is published in the Journal of Statistical Software. If you use EpiModel for any research or teaching purposes, please cite this reference:

Jenness SM, Goodreau SM, and Morris M. EpiModel: An R Package for Mathematical Modeling of Infectious Disease over Networks. Journal of Statistical Software. 2018; 84(8): 1-47. doi:10.18637/jss.v084.i08.

EpiModel is the foundation for a set of extension packages maintained by our group. EpiModelHIV (https://github.com/EpiModel/EpiModelHIV) models HIV and bacterial sexually transmitted infections, and EpiModelCOVID (https://github.com/EpiModel/EpiModelCOVID) models SARS-CoV-2. Simulation on high-performance computing clusters is supported by EpiModelHPC (https://github.com/EpiModel/EpiModelHPC).

See Also

Useful links:

• https://www.epimodel.org/

• https://epimodel.github.io/EpiModel/

• Report bugs at https://github.com/EpiModel/EpiModel/issues/

Description

This function performs a simple operation of updating the edgelist attribute n that tracks the total network size implicit in an edgelist representation of the network.

Usage

add_vertices(el, nv)

Arguments

Details

This function is used in EpiModel modules to add vertices (nodes) to the edgelist object to account for entries into the population (e.g., births and in-migration).

Value

Returns the matrix of current edges, el, with the population size attribute updated based on the number of new vertices specified in nv.

Examples

library("EpiModel")
nw <- network_initialize(100)
formation <- ~edges
target.stats <- 50
coef.diss <- dissolution_coefs(dissolution = ~offset(edges), duration = 20)
x <- netest(nw, formation, target.stats, coef.diss, verbose = FALSE)

param <- param.net(inf.prob = 0.3)
init <- init.net(i.num = 10)
control <- control.net(type = "SI", nsteps = 100, nsims = 5,
                       tergmLite = TRUE, resimulate.network = TRUE)

# Edgelist representation after initialization
crosscheck.net(x, param, init, control)
dat <- initialize.net(x, param, init, control, s = 1)
el <- get_edgelist(dat, network = 1)

# Check current network size
attributes(el)$n

# Add 10 vertices
el <- add_vertices(el, 10)

# Check new network size
attributes(el)$n

Description

Apportions a vector of values given a specified frequency distribution of those values such that the length of the output is robust to rounding and other instabilities.

Usage

apportion_lr(vector.length, values, proportions, shuffled = FALSE)

Arguments

Value

A vector of length vector.length containing the apportioned values from values.

Examples

## Not run: 
## Example 1: Without rounding
apportioned_vec_1 <- apportion_lr(4, c(1, 2, 3, 4, 5),
                                     c(0.25, 0, 0.25, 0.25, 0.25))

## Example 2: With rounding
apportioned_vec_2 <- apportion_lr(5, c(1, 2, 3, 4, 5),
                                     c(0.21, 0, 0.29, 0.25, 0.25))

## End(Not run)

Description

Arrive New Nodes to the netsim_dat Object

Usage

arrive_nodes(dat, nArrivals)

Arguments

Details

nArrivals new nodes are added to the network data stored on the netsim_dat object. If tergmLite is FALSE, these nodes are activated from the current timestep onward. Attributes for the new nodes must be set separately. On each clique layer (see netclique()), the new nodes are then placed under the layer's arrivals rule, which may set the layer's grouping attribute for them and add their edges. A model with an observed network layer (see netcensus()) has a fixed node set, and both this function and depart_nodes() stop if asked to change it.

Note that this function only supports arriving new nodes; returning to an active state nodes that were previously active in the network is not supported.

Value

the updated netsim_dat object with nArrivals new nodes added

Description

Converts an object that carries timed-edge information into a cumulative_edgelist: a data.frame with head, tail, start, stop columns matching the format used elsewhere in EpiModel. The typical use is feeding a networkDynamic object (from a non-tergmLite simulation, an external source, or a post-hoc construction) into the reachability tools.

Usage

as_cumulative_edgelist(x)

Arguments

Details

The networkDynamic method subtracts 1 from each spell's terminus so that the resulting start/stop range is inclusive on both ends, matching the convention used by get_cumulative_edgelist() and the reachability functions.

Value

A cumulative_edgelist object: a data.frame with at least the columns head, tail, start, stop. Edges are active from time start to time stop inclusive; NA in stop means the edge was never observed to dissolve.

See Also

get_forward_reachable() / get_backward_reachable() for typical downstream uses; dedup_cumulative_edgelist() when the source contains overlapping spells for the same pair of nodes.

Other cumulative_edgelist: dedup_cumulative_edgelist(), get_cumulative_degree(), get_cumulative_edgelist(), get_cumulative_edgelists_df(), get_partners(), reachable-nodes, update_cumulative_edgelist()

Examples

## Not run: 
# Convert a saved networkDynamic from a non-tergmLite netsim run
nd <- get_network(sim, sim = 1, network = 1)
el_cuml <- as_cumulative_edgelist(nd)
fwd <- get_forward_reachable(el_cuml, from_step = 1, to_step = 100)

## End(Not run)

Description

Convert an Edgelist into a Tibble

Usage

as_tibble_edgelist(el)

Arguments

Value

The edgelist in tibble form with two columns named head and tail.

Description

This function extracts a model run from an object of class dcm into a data frame using the generic as.data.frame function.

Usage

## S3 method for class 'dcm'
as.data.frame(x, row.names = NULL, optional = FALSE, run = NULL, ...)

Arguments

Details

Model output from dcm() simulations are available as a data frame with this helper function. The output data frame will include columns for time, the size of each compartment, the overall population size (the sum of compartment sizes), and the size of each flow.

For models with multiple runs (i.e., varying parameters - see example below), the default with the run parameter not specified will output all runs in a single stacked data frame.

Value

A data frame containing the data from x.

Examples

## Example 1: One-group SIS model with varying act.rate
param <- param.dcm(inf.prob = 0.2, act.rate = seq(0.05, 0.5, 0.05),
                   rec.rate = 1/50)
init <- init.dcm(s.num = 500, i.num = 1)
control <- control.dcm(type = "SIS", nsteps = 10)
mod1 <- dcm(param, init, control)
as.data.frame(mod1)
as.data.frame(mod1, run = 1)
as.data.frame(mod1, run = 10)

## Example 2: Two-group SIR model with vital dynamics
param <- param.dcm(inf.prob = 0.2, inf.prob.g2 = 0.1,
                   act.rate = 3, balance = "g1",
                   rec.rate = 1/50, rec.rate.g2 = 1/50,
                   a.rate = 1/100, a.rate.g2 = NA,
                   ds.rate = 1/100, ds.rate.g2 = 1/100,
                   di.rate = 1/90, di.rate.g2 = 1/90,
                   dr.rate = 1/100, dr.rate.g2 = 1/100)
init <- init.dcm(s.num = 500, i.num = 1, r.num = 0,
                 s.num.g2 = 500, i.num.g2 = 1, r.num.g2 = 0)
control <- control.dcm(type = "SIR", nsteps = 10)
mod2 <- dcm(param, init, control)
as.data.frame(mod2)

Description

This function extracts model simulations for objects of classes icm and netsim into a data frame using the generic as.data.frame function.

Usage

## S3 method for class 'icm'
as.data.frame(
  x,
  row.names = NULL,
  optional = FALSE,
  out = "vals",
  sim = NULL,
  qval = NULL,
  repair = "drop",
  ...
)

## S3 method for class 'netsim'
as.data.frame(
  x,
  row.names = NULL,
  optional = FALSE,
  out = "vals",
  sim = NULL,
  repair = "drop",
  ...
)

Arguments

Details

These methods work for both icm and netsim class models. The available output includes time-specific means, standard deviations, quantiles, and simulation values (compartment and flow sizes) from these stochastic model classes. Means, standard deviations, and quantiles are calculated by taking the row summary (i.e., each row of data corresponds to a time step) across all simulations in the model output.

Value

A data frame containing the data from x.

Examples

## Stochastic ICM SIS model
param <- param.icm(inf.prob = 0.8, act.rate = 2, rec.rate = 0.1)
init <- init.icm(s.num = 500, i.num = 1)
control <- control.icm(type = "SIS", nsteps = 10,
                       nsims = 3, verbose = FALSE)
mod <- icm(param, init, control)

# Default output all simulation runs, default to all in stacked data.frame
as.data.frame(mod)
as.data.frame(mod, sim = 2)

# Time-specific means across simulations
as.data.frame(mod, out = "mean")

# Time-specific standard deviations across simulations
as.data.frame(mod, out = "sd")

# Time-specific quantile values across simulations
as.data.frame(mod, out = "qnt", qval = 0.25)
as.data.frame(mod, out = "qnt", qval = 0.75)

## Not run: 
## Stochastic SI Network Model
nw <- network_initialize(n = 100)
formation <- ~edges
target.stats <- 50
coef.diss <- dissolution_coefs(dissolution = ~offset(edges), duration = 20)
est <- netest(nw, formation, target.stats, coef.diss, verbose = FALSE)

param <- param.net(inf.prob = 0.5)
init <- init.net(i.num = 10)
control <- control.net(type = "SI", nsteps = 10, nsims = 3, verbose = FALSE)
mod <- netsim(est, param, init, control)

# Same data extraction methods as with ICMs
as.data.frame(mod)
as.data.frame(mod, sim = 2)
as.data.frame(mod, out = "mean")
as.data.frame(mod, out = "sd")
as.data.frame(mod, out = "qnt", qval = 0.25)
as.data.frame(mod, out = "qnt", qval = 0.75)

## End(Not run)

Description

This function extracts timed edgelists for objects of class netdx into a data frame using the generic as.data.frame function.

Usage

## S3 method for class 'netdx'
as.data.frame(x, row.names = NULL, optional = FALSE, sim = NULL, ...)

Arguments

Value

A data frame containing the data from x.

Examples

# Initialize and parameterize the network model
nw <- network_initialize(n = 100)
formation <- ~edges
target.stats <- 50
coef.diss <- dissolution_coefs(dissolution = ~offset(edges), duration = 20)

# Model estimation
est <- netest(nw, formation, target.stats, coef.diss, verbose = FALSE)

# Simulate the network with netdx
dx <- netdx(est, nsims = 3, nsteps = 10, keep.tedgelist = TRUE,
            verbose = FALSE)

# Extract data from the first simulation
as.data.frame(dx, sim = 1)

# Extract data from all simulations
as.data.frame(dx)

Description

This methods ensures that the data.frame is correctly formatted as an epi.data.frame

Usage

as.epi.data.frame(df)

Arguments

Description

Converts a transmission matrix from the get_transmat function into a network::network class object.

Usage

## S3 method for class 'transmat'
as.network(x, ...)

Arguments

Details

When converting from a transmat to a network object, this functions copies the edge attributes within the transmission matrix ('at', 'infDur', 'transProb', 'actRate', and 'finalProb') into edge attributes on the network.

Value

A network::network object.

Description

Converts a transmission matrix from the get_transmat function into a phylo class object.

Usage

## S3 method for class 'transmat'
as.phylo(x, vertex.exit.times = NULL, ...)

Arguments

Details

Converts a transmat() object containing information about the history of a simulated infection into a ape::phylo object representation suitable for plotting as a tree with ape::plot.phylo(). Each infection event becomes a 'node' (horizontal branch) in the resulting phylo tree, and each network vertex becomes a 'tip' of the tree. The infection events are labeled with the vertex ID of the infector to make it possible to trace the path of infection.

The infection timing information is included to position the phylo-nodes, with the lines to the tips drawn to the max time value +1 (unless vertex.exit.times are passed in it effectively assumes all vertices are active until the end of the simulation).

If the transmat contains multiple infection seeds (there are multiple trees with separate root nodes), this function will return a list of class multiPhylo, each element of which is a phylo object. See ape::read.tree().

Value

A phylo class object.

Examples

set.seed(13)

# Fit a random mixing TERGM with mean degree of 1 and mean edge
# duration of 20 time steps
nw <- network_initialize(n = 100)
formation <- ~edges
target.stats <- 50
coef.diss <- dissolution_coefs(dissolution = ~offset(edges), duration = 20)
est <- netest(nw, formation, target.stats, coef.diss, verbose = FALSE)

# Parameterize the epidemic model as SI with one infected seed
param <- param.net(inf.prob = 0.5)
init <- init.net(i.num = 1)
control <- control.net(type = "SI", nsteps = 40, nsims = 1, verbose = FALSE)

# Simulate the model
mod1 <- netsim(est, param, init, control)

# Extract the transmission matrix
tm <- get_transmat(mod1)
head(tm, 15)

# Convert to phylo object and plot
tmPhylo <- as.phylo.transmat(tm)
oldpar <- par(mar = c(1,1,1,1))
plot(tmPhylo, show.node.label = TRUE,
              root.edge = TRUE,
              cex = 0.75)
par(oldpar)

Description

Assigns nodes to groups of a target size distribution, subject to a composition rule that keeps dependent members (such as children) in groups that contain an anchor member (such as an adult). This is the reverse of sample_groups(): the nodal attributes already exist, from a census age distribution for example, and the group ids are built to match them.

Usage

assign_groups(
  size.dist,
  role = NULL,
  anchor = NULL,
  dependent = NULL,
  n = NULL
)

Arguments

Details

The algorithm has four steps:

1. Group sizes are drawn from size.dist until they sum to at least n; the last group is trimmed so that the sizes sum to exactly n.

2. One anchor is placed in each group, in random order, until either the groups or the anchors run out. When there are fewer anchors than groups, the groups without an anchor cannot receive dependents.

3. Dependents are placed in the open slots of anchored groups, each open slot being equally likely, so that a group receives dependents in proportion to its remaining size. When the dependents outnumber those slots, the surplus is added to random anchored groups, which then exceed their drawn size; this is the only case in which the realized size distribution departs from size.dist, and a message reports it.

4. All remaining nodes (further anchors, and roles that are neither) fill the remaining open slots in the same way.

Without role, or without anchor, nodes are assigned to the drawn group sizes at random and no composition rule applies.

The rule guarantees cross-generational contact for every dependent. It does not otherwise control group composition: the number of dependents per group, the pairing of anchors, and whether a third role (older adults, for example) lives alone or with others all follow from the size distribution and the population shares of the roles. When the composition itself is the quantity to control, build the population from a table of group types with sample_groups() instead.

Value

An integer vector of group ids, one per node, with ids running from 1 to the number of groups.

See Also

sample_groups() and netclique().

Examples

# Ages from a census-like distribution, then households around them
set.seed(1)
n <- 1000
age <- sample(c("child", "adult", "elderly"), n, replace = TRUE,
              prob = c(0.22, 0.60, 0.18))
size.dist <- c("1" = 0.28, "2" = 0.34, "3" = 0.15, "4" = 0.14, "5" = 0.09)
hh_id <- assign_groups(size.dist, role = age,
                       anchor = c("adult", "elderly"), dependent = "child")

table(tabulate(hh_id))                          # realized size distribution
has_anchor <- tapply(age %in% c("adult", "elderly"), hh_id, any)
all(has_anchor[as.character(hh_id[age == "child"])])  # every child has one

nw <- network_initialize(n)
nw <- set_vertex_attribute(nw, "age", age)
nw <- set_vertex_attribute(nw, "hh_id", hh_id)
est_hh <- netclique(nw, group.attr = "hh_id")
print(est_hh, by = "age")

Description

Marks transmission probabilities, act rates, or recovery rates as varying with the time since infection, for the inf.prob, inf.prob.g2, act.rate, rec.rate, and rec.rate.g2 arguments of param.net(). The values are given either by stage, with the length of each stage in time steps, or one per time step since infection. The built-in netsim() modules apply the value at each infected node's current duration of infection.

Usage

by_infection_duration(x, durations = NULL)

Arguments

Details

A duration of infection is the number of time steps since the infected node was infected, counted from 1: a node infected at time step 10 is in its first step of infection at time steps 10 and 11, its second at time step 12, and so on. For transmission, the duration is that of the infected partner; for recovery, that of the recovering node.

The two forms describe the same thing. by_infection_duration(c(acute = 0.25, chronic = 0.02), durations = c(10, Inf)) gives 0.25 in the first ten time steps of an infection and 0.02 from the eleventh on, which can represent an acute stage of high infectiousness; by_infection_duration(c(rep(0.25, 10), 0.02)) is the same parameter given per time step. The stage form suits parameters reported by stage of infection; the per-step form suits a profile computed elsewhere, such as an infectiousness curve.

This is variation over the course of each infection, not over calendar time: to change a parameter at a given time step of the simulation, use scenarios or parameter updaters (see vignette("model-parameters", package = "EpiModel")).

In multi-layer models, an entry of a multilayer() parameter may itself be a by_infection_duration object: inf.prob = multilayer(by_infection_duration(c(0.5, 0.1)), 0.05). Arithmetic on the object, such as x * 2, scales the values and keeps the stages.

With the built-in modules, a plain vector passed to these arguments stops netsim() with an error, since a plain vector carries no record of what its positions mean; custom modules keep their own reading of plain vectors, such as values by group or by layer. Parameter tables cannot represent these objects yet, so param.net_to_table() stops on them.

Value

An object of class by_infection_duration: the numeric vector of values x, with the stage durations in its durations attribute. Without durations, every stage lasts one time step but the last.

See Also

param.net(), multilayer().

Examples

# An acute stage of ten time steps with a higher transmission probability
inf.prob <- by_infection_duration(c(acute = 0.25, chronic = 0.02),
                                  durations = c(10, Inf))
inf.prob

# Recovery impossible in the first 20 time steps of an infection, then
# certain
by_infection_duration(c(0, 1), durations = c(20, Inf))

# One value per time step, here an infectiousness profile
by_infection_duration(round(dgamma(1:15, shape = 3, rate = 0.5), 3))

param <- param.net(inf.prob = inf.prob, act.rate = 1)
param


nw <- network_initialize(n = 100)
est <- netest(nw, formation = ~edges, target.stats = 50,
              coef.diss = dissolution_coefs(~offset(edges), 10),
              verbose = FALSE)
sim <- netsim(est, param, init.net(i.num = 10),
              control.net(type = "SI", nsteps = 25, nsims = 1,
                          verbose = FALSE))
tm <- get_transmat(sim)
table(tm$infDur, tm$transProb)

Description

Checks for consistency in the implied network statistics of a two-group network in which the group size and group-specific degree distributions are specified.

Usage

check_degdist_bal(num.g1, num.g2, deg.dist.g1, deg.dist.g2)

Arguments

Details

This function outputs the number of nodes of degree 0 to g, where g is the length of a fractional degree distribution vector, given that vector and the size of the group. This utility is used to check for balance in implied degree given that fractional distribution within two-group network simulations, in which the degree-constrained counts must be equal across groups.

Examples

# An unbalanced distribution
check_degdist_bal(num.g1 = 500, num.g2 = 500,
                  deg.dist.g2 = c(0.40, 0.55, 0.03, 0.02),
                  deg.dist.g1 = c(0.48, 0.41, 0.08, 0.03))

# A balanced distribution
check_degdist_bal(num.g1 = 500, num.g2 = 500,
                  deg.dist.g1 = c(0.40, 0.55, 0.04, 0.01),
                  deg.dist.g2 = c(0.48, 0.41, 0.08, 0.03))

Description

Creates a new color-named temporally-extended attribute (TEA) variable in a networkDynamic object containing a disease status TEA in numeric format.

Usage

color_tea(
  nd,
  old.var = "testatus",
  old.sus = "s",
  old.inf = "i",
  old.rec = "r",
  new.var = "ndtvcol",
  new.sus = NULL,
  new.inf = NULL,
  new.rec = NULL,
  verbose = TRUE
)

Arguments

Details

The ndtv package (https://cran.r-project.org/package=ndtv) produces animated visuals for dynamic networks with evolving edge structures and nodal attributes. Nodal attribute dynamics in ndtv movies require a temporally extended attribute (TEA) containing a standard R color for each node at each time step. By default, the EpiModel package uses TEAs to store disease status history in network model simulations run in netsim(). But that status TEA is in numeric format (0, 1, 2). The color_tea function transforms those numeric values of that disease status TEA into a TEA with color values in order to visualize status changes in ndtv.

The convention in plot.netsim() is to color the susceptible nodes as blue, infected nodes as red, and recovered nodes as green. Alternate colors may be specified using the new.sus, new.inf, and new.rec parameters, respectively.

Using the color_tea function with a netsim object requires that TEAs for disease status be used and that the networkDynamic object be saved in the output: tergmListe must be set to FALSE in control.net().

Value

The updated object of class networkDynamic.

See Also

netsim() and the ndtv package documentation.

Description

Plots a compartment flow diagram for deterministic compartmental models, stochastic individual contact models, and stochastic network models.

Usage

comp_plot(x, at, digits, ...)

## S3 method for class 'netsim'
comp_plot(x, at = 1, digits = 3, ...)

## S3 method for class 'icm'
comp_plot(x, at = 1, digits = 3, ...)

## S3 method for class 'dcm'
comp_plot(x, at = 1, digits = 3, run = 1, ...)

Arguments

Details

The comp_plot function provides a visual summary of an epidemic model at a specific time step. The information contained in comp_plot is the same as in the summary functions for a model, but presented graphically as a compartment flow diagram.

For dcm class plots, specify the model run number if the model contains multiple runs, as in a sensitivity analysis. For icm and netsim class plots, the run argument is not used; the plots show the means and standard deviations across simulations at the specified time step.

These plots are currently limited to one-group models for each of the three model classes. That functionality may be expanded in future software releases.

Examples

## Example 1: DCM SIR model with varying act.rate
param <- param.dcm(inf.prob = 0.2, act.rate = 5:7,
                   rec.rate = 1/3, a.rate = 1/90, ds.rate = 1/100,
                   di.rate = 1/35, dr.rate = 1/100)
init <- init.dcm(s.num = 1000, i.num = 1, r.num = 0)
control <- control.dcm(type = "SIR", nsteps = 25, verbose = FALSE)
mod1 <- dcm(param, init, control)
comp_plot(mod1, at = 25, run = 3)

## Example 2: ICM SIR model with 3 simulations
param <- param.icm(inf.prob = 0.2, act.rate = 3, rec.rate = 1/50,
                   a.rate = 1/100, ds.rate = 1/100,
                   di.rate = 1/90, dr.rate = 1/100)
init <- init.icm(s.num = 500, i.num = 1, r.num = 0)
control <- control.icm(type = "SIR", nsteps = 25,
                       nsims = 3, verbose = FALSE)
mod2 <- icm(param, init, control)
comp_plot(mod2, at = 25, digits = 1)

Description

Sets the controls for deterministic compartmental models simulated with dcm().

Usage

control.dcm(
  type,
  nsteps,
  dt = 1,
  odemethod = "lsoda",
  dede = FALSE,
  new.mod = NULL,
  sens.param = TRUE,
  print.mod = FALSE,
  verbose = FALSE,
  ode.args = NULL,
  solver = NULL,
  ...
)

Arguments

Details

control.dcm sets the required control settings for any deterministic compartmental models solved with the dcm() function. Controls are required for both base model types and original models. For all base models, the type argument is a necessary parameter and it has no default.

Value

An EpiModel object of class control.dcm.

ODE Solver

As of EpiModel 2.5.5, the default ODE solver was changed from "rk4" (fixed-step Runge-Kutta 4) to "lsoda" (adaptive step-size). The "lsoda" method automatically adjusts its internal step size, providing numerical stability across a wide range of parameter values. The fixed-step "rk4" solver could produce numerical instability (negative compartment sizes, NaN values) for models with high transmission rates or stiff dynamics. Users who require the previous behavior can set odemethod = "rk4" explicitly.

New Model Functions

The form of the model function for base models may be displayed with the print.mod argument set to TRUE. In this case, the model will not be run. These model forms may be used as templates to write original model functions.

These new models may be input and solved with dcm() using the new.mod argument, which requires as input a model function.

Custom Solvers

By default, dcm() solves the model with one call to deSolve::ode over all time steps, passing odemethod and any ode.args. A model whose rates change abruptly (an intervention that starts on a given day, for example) is solved more accurately by restarting the solver at each change. The solver argument allows this, or any other solver, for base and new models alike. It is called with the initial state y, the solution times times, the model function func, the parameters parms, and the control settings control, and must return a matrix or data frame with time in the first column and one named column per state variable, as ode does. Any further columns it returns, such as derived outputs, are kept as model outputs. Columns whose names contain .flow are converted from cumulative counts to counts per time step after the solver returns, as for the default solver.

See Also

Use param.dcm() to specify model parameters and init.dcm() to specify the initial conditions. Run the parameterized model with dcm().

Description

Sets the controls for stochastic individual contact models simulated with icm().

Usage

control.icm(type, nsteps, nsims = 1, verbose = FALSE, verbose.int = 0)

Arguments

Details

control.icm sets the required control settings for any stochastic individual contact model solved with the icm() function. ICM simulations use the built-in SI, SIR, and SIS disease types only. The type argument is required and has no default. For custom or extension epidemic models, use the network model class via control.net() instead.

Value

An EpiModel object of class control.icm.

See Also

Use param.icm() to specify model parameters and init.icm() to specify the initial conditions. Run the parameterized model with icm().

Description

Sets the controls for stochastic network models simulated with netsim.

Usage

control.net(
  type,
  nsteps,
  start = 1,
  nsims = 1,
  ncores = 1,
  resimulate.network = FALSE,
  tergmLite = FALSE,
  edges.correct.attr = NULL,
  cumulative.edgelist = FALSE,
  truncate.el.cuml = 0,
  tracked.attributes = NULL,
  tracked.attributes.once = NULL,
  attr.rules,
  epi.by,
  initialize.FUN = initialize.net,
  resim_nets.FUN = resim_nets,
  summary_nets.FUN = summary_nets,
  infection.FUN = NULL,
  recovery.FUN = NULL,
  departures.FUN = NULL,
  arrivals.FUN = NULL,
  nwupdate.FUN = nwupdate.net,
  prevalence.FUN = prevalence.net,
  verbose.FUN = verbose.net,
  module.order = NULL,
  save.nwstats = TRUE,
  nwstats.formula = "formation",
  save.transmat = TRUE,
  save.network,
  save.run = FALSE,
  save.cumulative.edgelist = FALSE,
  save.other,
  verbose = TRUE,
  verbose.int = 1,
  skip.check = FALSE,
  raw.output = FALSE,
  future.use.plan = FALSE,
  tergmLite.track.duration = FALSE,
  set.control.ergm = control.simulate.formula(MCMC.burnin = 2e+05),
  set.control.tergm = control.simulate.formula.tergm(MCMC.maxchanges =
    .Machine$integer.max),
  save.diss.stats = TRUE,
  dat.updates = NULL,
  ...
)

Arguments

Details

control.net sets the required control settings for any network model solved with the netsim function. Controls are required for both base model types and when passing original process modules. For an overview of control settings for base models, consult the Network Modeling for Epidemics course materials For all base models, the type argument is a necessary parameter and it has no default.

Value

An EpiModel object of class control.net.

The attr.rules Argument

The attr.rules parameter is used to specify the rules for how nodal attribute values for incoming nodes should be set. These rules are only necessary for models in which there are incoming nodes (i.e., arrivals). There are three rules available for each attribute value:

• current: new nodes will be assigned this attribute in proportion to the distribution of that attribute among existing nodes at that current time step.

• t1: new nodes will be assigned this attribute in proportion to the distribution of that attribute among nodes at time 1 (that is, the proportions set in the original network for netest).

• Value: all new nodes will be assigned this specific value, with no variation. For example, the rules list attr.rules = list(race = "t1", sex = "current", status = "s") specifies how the race, sex, and status attributes should be set for incoming nodes. By default, the rule is "current" for all attributes except status, in which case it is "s" (that is, all incoming nodes are susceptible).

Checkpointing Simulations

netsim has a built-in checkpoint system to prevent losing computation work if the function is interrupted (SIGINT, power loss, time limit exceeded on a computation cluster). When enabled, each simulation will be saved every .checkpoint.steps time steps. Then, if a checkpoint enabled simulation is launched again with netsim, it will restart at the last checkpoint available in the saved data.

To enable the checkpoint capabilities of netsim, two control arguments have to be set: .checkpoint.steps, which is a positive number of time steps to be run between each file save; and .checkpoint.dir, which is the path to a directory to save the checkpointed data. If .checkpoint.dir directory does not exist, netsim will attempt to create it on the first checkpoint save. With these two controls defined, one can simply re-run netsim with the same arguments to restart a set of simulations that were interrupted.

Simulations are checkpointed individually: for example, if 3 simulations are run on a single core, the first 2 are finished, then the interruption occurs during the third, netsim will only restart the third one from the last checkpoint.

A .checkpoint.compress argument sets the compress argument of saveRDS used to save the checkpointed data. It defaults to TRUE, which is gzip compression and usually works well on netsim objects. Any value accepted by the compress argument of saveRDS may be used instead, such as "xz" or "bzip2" for smaller but slower files, or FALSE for no compression.

By default, if netsim reaches the end of all simulations, the checkpoint data directory and its content are removed before returning the netsim object. The .checkpoint.keep argument can be set to TRUE to prevent this removal to inspect the raw simulation objects.

New Modules

Base network models use a set of module functions that specify how the individual nodes in the network are subjected to infection, recovery, demographics, and other processes. Core modules are those listed in the .FUN arguments. For each module, there is a default function used in the simulation. The default infection module, for example, is contained in the infection.net function.

For original models, one may substitute replacement module functions for any of the default functions. New modules may be added to the workflow at each time step by passing a module function via the ... argument. Consult the Extending EpiModel section of the Network Modeling for Epidemics course materials. One may remove existing modules, such as arrivals.FUN, from the workflow by setting the parameter value for that argument to NULL.

End Horizon

netsim implements an "End Horizon" mechanism, where a set of modules are removed from the simulation at a specific time step. This is enabled through the end.horizon parameter to control.net.

This parameter must receive a list with fields at, the time step at which the end horizon occurs, and modules, a character vector with the names of the modules to remove. For example, list(at = 208, modules = c("arrivals.FUN", "infections.FUN")).

See Also

Use param.net to specify model parameters and init.net to specify the initial conditions. Run the parameterized model with netsim.

Description

This helper function populates a netsim_dat main list object with the minimal required elements. All parameters are optional. When none are given the resulting object is only a shell list of class netsim_dat with the different named elements defined as empty lists.

Usage

create_dat_object(
  param = list(),
  init = list(),
  control = list(),
  run = list()
)

Arguments

Value

A netsim_dat main list object.

Description

An EpiModel scenario allows one or multiple set of parameters to be applied to a model a predefined timesteps. They are usually used by a researcher who wants to model counterfactuals using a pre calibrated model.

Usage

create_scenario_list(scenarios.df)

Arguments

Value

a list of EpiModel scenarios

scenarios.df

The scenarios.df is a data.frame of values to be used as parameters.

It must contain a ".at" column, specifying when the changes should occur. It requires the "updater" module of EpiModel. See, vignette. If the ".at" value of a row is less than two, the changes will be applied to the parameter list iteself. The second mandatory column is ".scenario.id". It is used to distinguish the different scenarios. If multiple rows share the same ".scenario.id", the resulting scenario will contain one updater per row. This permits modifying parameters at multiple points in time. (e.g. an intervention limited in time).

The other column names must correspond either to: the name of one parameter if this parameter is of size 1 or the name of the parameter with "_1", "_2", "N" with the second part being the position of the value for a parameter of size > 1. This means that the parameter names cannot contain any underscore "". (e.g "a.rate", "d.rate_1", "d.rate_2")

Vector parameters

A column with a position suffix sets that position of the current value of the parameter and leaves the other positions as they are. For example, with prep.start.rate = c(0.0055, 0.0042, 0.0066) in param, a scenario with only a prep.start.rate_1 column of 0.5 gives prep.start.rate = c(0.5, 0.0042, 0.0066), both when it is applied before the simulation by use_scenario() and when it is applied by an updater during the simulation. Before EpiModel 2.7.0, the parameter was replaced by a vector built from the scenario columns alone (0.5 here), so every position had to be given.

A column without a suffix (prep.start.rate) still replaces the whole parameter with its single value. Positions past the end of the current vector extend it when they leave no gap: with a vector of length 3, prep.start.rate_4 alone extends it to length 4, while prep.start.rate_5 alone is an error. A scenario cannot shorten a vector. Only plain vectors can be set by position: a parameter given as a multilayer() or by_infection_duration() object, or as a list, is an error. Use a parameter updater that sets the whole value instead.

Description

Solves deterministic compartmental epidemic models for infectious disease.

Usage

dcm(param, init, control)

Arguments

Details

The dcm function uses the ordinary differential equation solver in the deSolve package to model disease as a deterministic compartmental system. The parameterization for these models follows the standard approach in EpiModel, with epidemic parameters, initial conditions, and control settings.

The dcm function performs modeling of both base model types and original models with new structures. Base model types include one-group and two-group models with disease types for Susceptible-Infected (SI), Susceptible-Infected-Recovered (SIR), and Susceptible-Infected-Susceptible (SIS). Both base and original models require the param, init, and control inputs.

Solver settings such as tolerances are passed to deSolve through the ode.args argument of control.dcm(), and the solver argument replaces the solver altogether, for example to restart it at abrupt changes in a model's rates.

Value

A list of class dcm with the following elements:

• param: the epidemic parameters passed into the model through param, with additional parameters added as necessary.

• control: the control settings passed into the model through control, with additional controls added as necessary.

• epi: a list of data frames, one for each epidemiological output from the model. Outputs for base models always include the size of each compartment, as well as flows in, out of, and between compartments.

References

Soetaert K, Petzoldt T, Setzer W. Solving Differential Equations in R: Package deSolve. Journal of Statistical Software. 2010; 33(9): 1-25. doi:10.18637/jss.v033.i09.

See Also

Extract the model results with as.data.frame.dcm(). Summarize the time-specific model results with summary.dcm(). Plot the model results with plot.dcm(). Plot a compartment flow diagram with comp_plot().

Examples

## Example 1: SI Model (One-Group)
# Set parameters
param <- param.dcm(inf.prob = 0.2, act.rate = 0.25)
init <- init.dcm(s.num = 500, i.num = 1)
control <- control.dcm(type = "SI", nsteps = 500)
mod1 <- dcm(param, init, control)
mod1
plot(mod1)

## Example 2: SIR Model with Vital Dynamics (One-Group)
param <- param.dcm(inf.prob = 0.2, act.rate = 5,
                   rec.rate = 1/3, a.rate = 1/90, ds.rate = 1/100,
                   di.rate = 1/35, dr.rate = 1/100)
init <- init.dcm(s.num = 500, i.num = 1, r.num = 0)
control <- control.dcm(type = "SIR", nsteps = 500)
mod2 <- dcm(param, init, control)
mod2
plot(mod2)

## Example 3: SIS Model with act.rate Sensitivity Parameter
param <- param.dcm(inf.prob = 0.2, act.rate = seq(0.1, 0.5, 0.1),
                   rec.rate = 1/50)
init <- init.dcm(s.num = 500, i.num = 1)
control <- control.dcm(type = "SIS", nsteps = 500)
mod3 <- dcm(param, init, control)
mod3
plot(mod3)

## Example 4: SI Model with Vital Dynamics (Two-Group)
param <- param.dcm(inf.prob = 0.4,  inf.prob.g2 = 0.1,
                   act.rate = 0.25, balance = "g1",
                   a.rate = 1/100, a.rate.g2 = NA,
                   ds.rate = 1/100, ds.rate.g2 = 1/100,
                   di.rate = 1/50, di.rate.g2 = 1/50)
init <- init.dcm(s.num = 500, i.num = 1,
                 s.num.g2 = 500, i.num.g2 = 0)
control <- control.dcm(type = "SI", nsteps = 500)
mod4 <- dcm(param, init, control)
mod4
plot(mod4)

## Example 5: SI Model with Varying Initial Conditions
param <- param.dcm(inf.prob = 0.2, act.rate = 0.25)
init <- init.dcm(s.num = 500, i.num = c(1, 5, 25))
control <- control.dcm(type = "SI", nsteps = 500)
mod5 <- dcm(param, init, control)
mod5
plot(mod5)

Description

Merges overlapping spells for the same ⁠(head, tail)⁠ pair into a single row. Useful when concatenating cumulative edgelists across multiple simulations or sources, where the same partnership can appear in more than one segment. The reachability functions assume non-overlapping input.

Usage

dedup_cumulative_edgelist(el)

Arguments

Value

A cumulative edgelist with no overlapping edges for the same ⁠(head, tail)⁠ pair.

See Also

as_cumulative_edgelist() for one source of duplicates; get_forward_reachable() / get_backward_reachable(), which expect non-overlapping input.

Other cumulative_edgelist: as_cumulative_edgelist(), get_cumulative_degree(), get_cumulative_edgelist(), get_cumulative_edgelists_df(), get_partners(), reachable-nodes, update_cumulative_edgelist()

Examples

## Not run: 
el_cuml <- dedup_cumulative_edgelist(el_raw)

## End(Not run)

Description

Given a current two-column matrix of edges and a vector of vertex IDs, this function removes any rows of the edgelist in which the IDs are present.

Usage

delete_edges(el, vid)

Arguments

Value

Returns an updated edgelist object, with any edges including the specified vertices removed.

Description

Given a current two-column matrix of edges and a vector of IDs to delete from the matrix, this function first removes any rows of the edgelist in which the IDs are present and then permutes downward the index of IDs on the edgelist that were numerically larger than the IDs deleted.

Usage

delete_vertices(el, vid)

Arguments

Details

This function is used in EpiModel modules to remove vertices (nodes) from the edgelist object to account for exits from the population (e.g., deaths and out-migration).

Value

Returns an updated edgelist object, el, with the edges of deleted vertices removed from the edgelist and the ID numbers of the remaining edges permuted downward.

Examples

library("EpiModel")
set.seed(12345)
nw <- network_initialize(100)
formation <- ~edges
target.stats <- 50
coef.diss <- dissolution_coefs(dissolution = ~offset(edges), duration = 20)
x <- netest(nw, formation, target.stats, coef.diss, verbose = FALSE)

param <- param.net(inf.prob = 0.3)
init <- init.net(i.num = 10)
control <- control.net(type = "SI", nsteps = 100, nsims = 5,
                       tergmLite = TRUE, resimulate.network = TRUE)

# Set seed for reproducibility
set.seed(123456)

# Edgelist representation after initialization
crosscheck.net(x, param, init, control)
dat <- initialize.net(x, param, init, control, s = 1)
el <- get_edgelist(dat, network = 1)

# Current edges
head(el, 20)

# Remove nodes 1 and 2
nodes.to.delete <- 1:2
el <- delete_vertices(el, nodes.to.delete)

# Newly permuted edges
head(el, 20)

Description

Depart Nodes from the netsim_dat Object

Usage

depart_nodes(dat, departures)

Arguments

Details

If tergmLite is FALSE, the vertex ids departures are deactivated (from the current timestep onward) in each networkDynamic stored in dat$run$nw. If tergmLite is TRUE, the vertex ids departures are deleted from dat$run$el, dat$run$attr, and dat$run$net_attr.

Value

the updated netsim_dat object with the nodes in departures departed

Description

Calculates dissolution coefficients, given a dissolution model and average edge duration, to pass as offsets to an ERGM/TERGM model fit in netest.

Usage

dissolution_coefs(dissolution, duration, d.rate = 0)

Arguments

Details

This function performs two calculations for dissolution coefficients used in a network model estimated with netest():

1. Transformation: the mean durations of edges in a network are mathematically transformed to logit coefficients.

2. Adjustment: in a dynamic network simulation in an open population (in which there are departures), it is further necessary to adjust these coefficients; this upward adjustment accounts for departure as a competing risk to edge dissolution.

The current dissolution models supported by this function and in network model estimation in netest() are as follows:

• ~offset(edges): a homogeneous dissolution model in which the edge duration is the same for all partnerships. This requires specifying one duration value.

• ~offset(edges) + offset(nodematch("<attr>")): a heterogeneous model in which the edge duration varies by whether the nodes in the dyad have similar values of a specified attribute. The duration vector should now contain two values: the first is the mean edge duration of non-matched dyads, and the second is the duration of the matched dyads.

• ~offset(edges) + offset(nodemix("<attr>")): a heterogeneous model that extends the nodematch model to include non-binary attributes for homophily. The duration vector should first contain the base value, then the values for every other possible combination in the term.

Value

A list of class disscoef with the following elements:

• dissolution: right-hand sided STERGM dissolution formula passed in the function call.

• duration: mean edge durations passed into the function.

• coef.crude: mean durations transformed into logit coefficients.

• coef.adj: crude coefficients adjusted for the risk of departure on edge persistence, if the d.rate argument is supplied.

• coef.form.corr: corrections to be subtracted from formation coefficients.

• d.rate: the departure rate.

• diss.model.type: the form of the dissolution model; options include edgesonly, nodematch, and nodemix.

Examples

## Homogeneous dissolution model with no departures
dissolution_coefs(dissolution = ~offset(edges), duration = 25)

## Homogeneous dissolution model with departures
dissolution_coefs(dissolution = ~offset(edges), duration = 25,
                  d.rate = 0.001)

## Heterogeneous dissolution model in which same-race edges have
## shorter duration compared to mixed-race edges, with no departures
dissolution_coefs(dissolution = ~offset(edges) + offset(nodematch("race")),
                  duration = c(20, 10))

## Heterogeneous dissolution model in which same-race edges have
## shorter duration compared to mixed-race edges, with departures
dissolution_coefs(dissolution = ~offset(edges) + offset(nodematch("race")),
                  duration = c(20, 10), d.rate = 0.001)

## Not run: 
## Extended example for differential homophily by age group
# Set up the network with nodes categorized into 5 age groups
nw <- network_initialize(n = 1000)
age.grp <- sample(1:5, 1000, TRUE)
nw <- set_vertex_attribute(nw, "age.grp", age.grp)

# durations = non-matched, age.grp1 & age.grp1, age.grp2 & age.grp2, ...
# TERGM will include differential homophily by age group with nodematch term
# Target stats for the formation model are overall edges, and then the number
# matched within age.grp 1, age.grp 2, ..., age.grp 5
form <- ~edges + nodematch("age.grp", diff = TRUE)
target.stats <- c(450, 100, 125, 40, 80, 100)

# Target stats for the dissolution model are duration of non-matched edges,
# then duration of edges matched within age.grp 1, age.grp 2, ..., age.grp 5
durs <- c(60, 30, 80, 100, 125, 160)
diss <- dissolution_coefs(~offset(edges) +
                            offset(nodematch("age.grp", diff = TRUE)),
                          duration = durs)

# Fit the TERGM
fit <- netest(nw, form, target.stats, diss)

# Full diagnostics to evaluate model fit
dx <- netdx(fit, nsims = 10, ncores = 4, nsteps = 300)
print(dx)

# Simulate one long time series to examine timed edgelist
dx <- netdx(fit, nsims = 1, nsteps = 5000, keep.tedgelist = TRUE)

# Extract timed-edgelist
te <- as.data.frame(dx)
head(te)

# Limit to non-censored edges
te <- te[which(te$onset.censored == FALSE & te$terminus.censored == FALSE),
         c("head", "tail", "duration")]
head(te)

# Look up the age group of head and tail nodes
te$ag.head <- age.grp[te$head]
te$ag.tail <- age.grp[te$tail]
head(te)

# Recover average edge durations for age-group pairing
mean(te$duration[te$ag.head != te$ag.tail])
mean(te$duration[te$ag.head == 1 & te$ag.tail == 1])
mean(te$duration[te$ag.head == 2 & te$ag.tail == 2])
mean(te$duration[te$ag.head == 3 & te$ag.tail == 3])
mean(te$duration[te$ag.head == 4 & te$ag.tail == 4])
mean(te$duration[te$ag.head == 5 & te$ag.tail == 5])
durs

## End(Not run)

Description

Outputs a table of the number and percent of edges that are left-censored, right-censored, both-censored, or uncensored for a networkDynamic object.

Usage

edgelist_censor(el)

Arguments

Details

Given a STERGM simulation over a specified number of time steps, the edges within that simulation may be left-censored (started before the first step), right-censored (continued after the last step), right and left-censored, or uncensored. The amount of censoring will increase when the average edge duration approaches the length of the simulation.

Value

A 4 x 2 table containing the number and percent of edges in el that are left-censored, right-censored, both-censored, or uncensored.

Examples

# Initialize and parameterize network model
nw <- network_initialize(n = 100)
formation <- ~edges
target.stats <- 50
coef.diss <- dissolution_coefs(dissolution = ~offset(edges), duration = 20)

# Model estimation
est <- netest(nw, formation, target.stats, coef.diss, verbose = FALSE)

# Simulate the network and extract a timed edgelist
dx <- netdx(est, nsims = 1, nsteps = 100, keep.tedgelist = TRUE,
      verbose = FALSE)
el <- as.data.frame(dx)

# Calculate censoring
edgelist_censor(el)

Description

Launches interactive Shiny applications for exploring deterministic compartmental models (DCMs) and stochastic individual contact models (ICMs) of infectious disease transmission.

Usage

epiweb(class, ...)

Arguments

Details

epiweb launches a web-based graphical interface for configuring, running, and analyzing epidemic models built on EpiModel. Each application provides interactive controls for model parameters, real-time (or on-demand) visualization with plotly, narrative model summaries, and downloadable data tables. Both applications include a Guide tab with a comprehensive user manual covering the underlying theory, parameter definitions, and instructions for interpreting output.

DCM Application (class = "dcm")

The DCM app solves systems of ordinary differential equations for one-group SI, SIR, and SIS models. Models run instantly and the interface is fully reactive (no run button). Features include:

• Scenario presets: Flu-like (SIR), STI-like (SIS), and Measles-like (SIR) configurations that set all parameters to epidemiologically plausible values.

• Interventions: An optional mid-epidemic intervention that reduces the transmission probability by a user-specified efficacy, shown as a vertical dashed line on the plot.

• Sensitivity analysis: Vary any epidemic parameter across a range of values to compare trajectories side by side with a blue-to-red color ramp.

• Vital dynamics: Optional births and deaths to model endemic equilibria over longer time horizons.

• Summary tab: Displays the basic reproduction number (R₀), peak timing, cumulative infections, attack rate, and contextual interpretation.

• Data tab: Searchable table of model output with CSV download.

ICM Application (class = "icm")

The ICM app runs agent-based stochastic microsimulations for one-group SI, SIR, and SIS models. Because ICMs are computationally more expensive, a Run Model button controls when simulations execute. Features include:

• Stochastic visualization: Multiple simulations are summarized with mean trajectory lines and interquartile range (IQR) ribbons showing the middle 50\ individual trajectories.

• Scenario presets: The same Flu-like, STI-like, and Measles-like configurations as the DCM app, adapted for smaller populations suitable for agent-based simulation.

• Interventions: Same intervention system as DCM, natively supported by param.icm.

• Vital dynamics: Optional births and deaths, as in DCM.

• Summary tab: Displays R₀, epidemic timeline from the mean trajectory, and a stochastic variation section showing the range and standard deviation of peak infections and final prevalence across simulations.

• Data tab: View simulation output as means, standard deviations, or individual simulation values, with CSV download.

Required Packages

Both applications require the shiny, bslib, DT, and plotly packages. Install any missing packages with install.packages(c("shiny", "bslib", "DT", "plotly")).

See Also

dcm and param.dcm for the DCM modeling API; icm and param.icm for the ICM modeling API.

Examples

## Not run: 
## Deterministic compartmental models
epiweb(class = "dcm")

## Stochastic individual contact models
epiweb(class = "icm")

## End(Not run)

Description

Plots quantile bands given a data.frame with stochastic model results from icm() or netsim().

Usage

geom_bands(mapping, lower = 0.25, upper = 0.75, alpha = 0.25, ...)

Arguments

Details

This is a wrapper around ggplot::stat_summary with a ribbon geom as aesthetic output.

Examples

param <- param.icm(inf.prob = 0.2, act.rate = 0.25)
init <- init.icm(s.num = 500, i.num = 1)
control <- control.icm(type = "SI", nsteps = 250, nsims = 5)
mod1 <- icm(param, init, control)
df <- as.data.frame(mod1)
df.mean <- as.data.frame(mod1, out = "mean")

library(ggplot2)
ggplot() +
   geom_line(data = df, mapping = aes(time, i.num, group = sim),
   alpha = 0.25, lwd = 0.25, color = "firebrick") +
   geom_bands(data = df, mapping = aes(time, i.num),
              lower = 0.1, upper = 0.9, fill = "firebrick") +
   geom_line(data = df.mean, mapping = aes(time, i.num)) +
   theme_minimal()

Description

Returns an adjacency list from an edge list

Usage

get_adj_list(el, n_nodes)

Arguments

Details

The adjacency list is a list of length n_nodes. The entry for each node is a integer vector containing the index of all the nodes connected to it. This layout makes it directly subsetable in O(1) at the expanse of memory usage. To get all connections to the nodes 10 and 15 : ⁠unlist(adj_list[c(10, 15)]⁠

Value

An adjacency list for the network

Description

Retrieve, create, and edit the nodal attributes of the main netsim_dat object of network models, from inside module functions. Each nodal attribute is a vector with one value per node, indexed by the positional identifiers of the nodes (posit_ids).

• get_attr() returns one attribute, for all nodes or for the nodes in posit_ids. get_attr_list() returns several attributes as a named list.

• set_attr() replaces the values of an attribute, for all nodes or for the nodes in posit_ids, and creates the attribute if it does not exist. add_attr() creates an attribute set to NA for all nodes.

• append_core_attr() adds n.new nodes and sets their core attributes. append_attr() then adds the values of each other attribute for the same nodes.

• remove_node_attr() removes the nodes in posit_ids from every attribute.

• get_core_attributes() lists the core attributes with their types.

See net-accessor for an overview of the accessors.

Usage

get_attr(dat, item, posit_ids = NULL, override.null.error = FALSE)

get_attr_list(dat, item = NULL)

set_attr(dat, item, value, posit_ids = NULL, override.length.check = FALSE)

add_attr(dat, item)

append_attr(dat, item, value, n.new)

append_core_attr(dat, at, n.new)

remove_node_attr(dat, posit_ids)

get_core_attributes()

Arguments

Value

For get_attr(), the attribute vector, or its values for the nodes in posit_ids; for get_attr_list(), a named list of attribute vectors; for get_core_attributes(), a named list with a value of the correct type for each core attribute. The other functions return the updated netsim_dat object.

Core Attributes

The append_core_attr function initializes the attributes necessary for EpiModel to work (the four core attributes are: "active", "unique_id", "entrTime", and "exitTime"). These attributes are used in the initialization phase of the simulation, to create the nodes (see initialize.net()); and also used when adding nodes during the simulation (see arrivals.net()).

One Value Per Node

Every nodal attribute holds exactly one value per node, and netsim() checks this at the end of each time step. When set_attr and append_attr create an attribute, the nodes they do not write to are filled with NA to preserve that: set_attr with posit_ids fills the unselected nodes, and append_attr fills the nodes that predate the ones being appended. This means append_core_attr must be called before append_attr when adding new nodes, as it is what makes those nodes known to the rest of the object.

Intended Usage

These accessors operate on the live netsim_dat object that EpiModel constructs internally when a simulation starts and passes through each module at every time step. The expected calling context is inside a user-supplied module function (e.g., a custom infection.FUN, recovery.FUN, or arrivals/departures module) registered with control.net().

They are not a general-purpose tool for editing the user-facing input or output objects of netsim():

• Calling them on a param.net(), init.net(), or control.net() object before a run will appear to succeed (those objects are also list-like) but produces an object that is not a valid simulation input.

• Calling them on the object returned by netsim() modifies only the saved input slot, not anything that would change the result of a re-run. Note also that netsim() strips the ⁠*.FUN⁠ entries from its returned control slot, so feeding that slot back into a new simulation will fail.

For modifications outside a running simulation, use the following instead:

• Editing parameters before a run: update_params() for a param.net object, or direct list assignment (e.g., p$inf.prob <- 0.5).

• Editing init or control before a run: direct list assignment (e.g., ctrl$nsteps <- 1000), or rebuild with a fresh call to the init.net() / control.net() constructor.

• Scheduled mid-simulation changes: the .param.updater.list argument to param.net() and the .control.updater.list argument to control.net(). See the "model-parameters" vignette for the underlying updater module and the scenario API built on top of it.

Mutability

The set_, append_, and add_ functions DO NOT modify the netsim_dat object in place. The result must be assigned back to dat in order to be registered: ⁠dat <- set_*(dat, item, value)⁠.

See Also

net-accessor for an overview of the accessors, and get_epi() and get_param() for the other groups; arrive_nodes() and depart_nodes() for adding and removing nodes from the networks; get_unique_ids() and get_posit_ids() for converting node identifiers.

Examples

dat <- create_dat_object(control = list(nsteps = 150))
dat <- append_core_attr(dat, at = 1, n.new = 10)

dat <- add_attr(dat, "age")
dat <- set_attr(dat, "age", runif(10, 15, 65))
dat <- set_attr(dat, "status", rep(c("s", "i"), c(8, 2)))
get_attr_list(dat)
get_attr_list(dat, c("age", "active"))
get_attr(dat, "status")
get_attr(dat, "status", posit_ids = c(1, 9))

# Edit selected nodes
dat <- set_attr(dat, "status", "i", posit_ids = 3)
get_attr(dat, "status")

# Add 5 nodes: core attributes first, then the others
dat <- append_core_attr(dat, at = 2, n.new = 5)
dat <- append_attr(dat, "status", "s", 5)
dat <- append_attr(dat, "age", 15, 5)
get_attr(dat, "age")
get_attr(dat, "entrTime")

Description

Reconstructs the state of all tracked attributes at a specific time step from the delta-compressed attribute history. Returns a tibble with one row per active node and one column per tracked attribute.

Usage

get_attr_at(sim, at, sim_num = 1, compute_age = FALSE)

Arguments

Value

A tibble with columns unique_id, active, and one column per tracked attribute, containing the most recent value at or before at for each active node.

See Also

record_attr_history(), get_attr_history()

Examples

nw <- network_initialize(n = 100)
est <- netest(nw, formation = ~edges, target.stats = 50,
              coef.diss = dissolution_coefs(~offset(edges), 10),
              verbose = FALSE)
control <- control.net(type = "SI", nsteps = 20, nsims = 1,
                       tracked.attributes = "status", verbose = FALSE)
sim <- netsim(est, param.net(inf.prob = 0.3, act.rate = 1),
              init.net(i.num = 10), control)

# The status of every node at time step 10
at10 <- get_attr_at(sim, at = 10)
head(at10)
table(at10$status)

Description

Extract the Attributes History from Network Simulations

Usage

get_attr_history(sims)

Arguments

Value

A list of data.frames, one for each "measure" recorded in the simulation by the record_attr_history function.

Examples

nw <- network_initialize(n = 100)
est <- netest(nw, formation = ~edges, target.stats = 50,
              coef.diss = dissolution_coefs(~offset(edges), 10),
              verbose = FALSE)
control <- control.net(type = "SI", nsteps = 20, nsims = 1,
                       tracked.attributes = "status", verbose = FALSE)
sim <- netsim(est, param.net(inf.prob = 0.3, act.rate = 1),
              init.net(i.num = 10), control)

# One data frame per tracked attribute, with the changes at each step
hist <- get_attr_history(sim)
names(hist)
head(hist$status)

Description

Returns all the node connected directly or indirectly to a set of nodes

Usage

get_connected_nodes(adj_list, nodes)

Arguments

Value

A vector of nodes indexes that are connected together with the ones provided in the nodes argument. The nodes themselves are not listed in this output

Description

Counts the number of distinct partners each of a set of index nodes has accumulated over the tracked history. A thin wrapper around get_partners() that collapses the partner list to a count per index.

Usage

get_cumulative_degree(
  dat,
  index_posit_ids,
  networks = NULL,
  truncate = Inf,
  only.active.nodes = FALSE
)

Arguments

Value

A data.frame with 2 columns:

• index_pid: the positional ID (see get_posit_ids()) of the indexes.

• degree: the cumulative degree of the index.

Cumulative Degree

The cumulative degree of a node is the number of distinct partners connected to it during the tracked time window. The window is whatever history the cumulative edgelist currently contains (controlled at simulation time by control$truncate.el.cuml and at call time by the truncate argument).

See Also

get_partners() for the full partner list backing this count.

Other cumulative_edgelist: as_cumulative_edgelist(), dedup_cumulative_edgelist(), get_cumulative_edgelist(), get_cumulative_edgelists_df(), get_partners(), reachable-nodes, update_cumulative_edgelist()

Examples

## Not run: 
get_cumulative_degree(dat, index_posit_ids = 1:50, truncate = 52)

## End(Not run)

Description

Returns the cumulative edgelist for a given network: a record of every edge EpiModel has tracked, with the time steps at which each edge formed and (if no longer active) dissolved. This is the canonical mechanism for querying partnership histories during or after a simulation, and is particularly important under tergmLite = TRUE, where the full networkDynamic history is not retained.

Usage

get_cumulative_edgelist(dat, network)

Arguments

Details

Cumulative-edgelist tracking is opt-in. It must be enabled via control.net(cumulative.edgelist = TRUE) (see control.net()); calling this function on a simulation that did not enable tracking raises an error.

Inside a custom module (where dat is the live netsim_dat object), call this function directly. After a netsim() run, the saved edgelists are attached to the returned object as sim$cumulative.edgelist[[s]] (one element per simulation, populated when save.cumulative.edgelist = TRUE); read those directly rather than calling this function on the processed output.

For the current (active-only, no history) edgelist, see get_edgelist() and get_edgelists_df().

Value

A tibble with four columns:

• head: the unique ID (see get_unique_ids()) of the head node.

• tail: the unique ID of the tail node.

• start: the time step at which the edge formed.

• stop: the time step at which the edge dissolved, or NA if the edge is still active. Edges are active over ⁠[start, stop]⁠ inclusive.

See Also

control.net() for the controlling flags (cumulative.edgelist, truncate.el.cuml, save.cumulative.edgelist). vignette("network-objects", package = "EpiModel") walks through the full lifecycle.

Other cumulative_edgelist: as_cumulative_edgelist(), dedup_cumulative_edgelist(), get_cumulative_degree(), get_cumulative_edgelists_df(), get_partners(), reachable-nodes, update_cumulative_edgelist()

Examples

## Not run: 
# Inside a custom module (dat is the live netsim_dat object):
el <- get_cumulative_edgelist(dat, network = 1)

# Post-simulation on a processed netsim object, read the saved slot:
sim$cumulative.edgelist[[1]]

## End(Not run)

Description

Combines the cumulative edgelists from one or more network layers into a single data.frame, adding a network column to identify the layer. Like get_cumulative_edgelist(), this requires control.net(cumulative.edgelist = TRUE).

Usage

get_cumulative_edgelists_df(dat, networks = NULL)

Arguments

Value

A data.frame with 5 columns:

• head: the unique ID (see get_unique_ids()) of the head node.

• tail: the unique ID of the tail node.

• start: the time step at which the edge formed.

• stop: the time step at which the edge dissolved, or NA if still active.

• network: the numerical index of the network on which the partnership lives.

Note: column names head/tail here match the single-network get_cumulative_edgelist() output. The index/partner naming used by get_partners() is reserved for queries about specific index nodes.

See Also

get_cumulative_edgelist() (single network), control.net() (cumulative.edgelist, save.cumulative.edgelist). vignette("network-objects", package = "EpiModel") walks through the full lifecycle.

Other cumulative_edgelist: as_cumulative_edgelist(), dedup_cumulative_edgelist(), get_cumulative_degree(), get_cumulative_edgelist(), get_partners(), reachable-nodes, update_cumulative_edgelist()

Description

Return the Current Timestep

Usage

get_current_timestep(dat)

Arguments

Value

The current timestep.

Description

A fast method for querying the current degree of all individuals within a network.

Usage

get_degree(x)

Arguments

Details

Individual-level data on the current degree of nodes within a network is often useful for summary statistics. Given a network class object, net, one way to look up the current degree is to get a summary of the ERGM term, sociality, as in: summary(net ~ sociality(nodes = NULL)). But that is computationally inefficient for a number of reasons. This function provides a fast method for generating the vector of degrees using a query of the edgelist. It is even faster if the parameter x is already transformed into an edgelist.

Value

A vector of length equal to the total network size, containing the current degree of each node in the network.

Examples

nw <- network_initialize(n = 500)

set.seed(1)
fit <- ergm(nw ~ edges, target.stats = 250)
sim <- simulate(fit)

# Slow ERGM-based method
ergm.method <- unname(summary(sim ~ sociality(nodes = NULL)))
ergm.method

# Fast tabulate method with network object
deg.net <- get_degree(sim)
deg.net

# Even faster if network already transformed into an edgelist
el <- as.edgelist(sim)
deg.el <- get_degree(el)
deg.el

identical(as.integer(ergm.method), deg.net, deg.el)

Description

Calculates the cumulative or momentary degree distribution from the timed edgelist of a dynamic netdx() simulation.

Usage

get_degree_dist(
  x,
  type = c("cumulative", "momentary"),
  sims = NULL,
  window = NULL,
  unique.partners = TRUE
)

Arguments

Details

The momentary degree distribution is the distribution plotted by plot.netdx() when degree(k) terms are included in the nwstats.formula of netdx(): it counts the partners a node has at one moment in time. The cumulative degree distribution counts the partners a node accumulates across the whole observation window, so a node with one partner at a time and a node with several simultaneous partners may have the same cumulative degree.

Partnerships are counted if they are active at any point in the window: a partnership contributes to the cumulative degree of both of its nodes when its spell overlaps window, and to the momentary degree of both nodes at each time step of the window in which it is active. Partnerships that are ongoing when the simulation ends are counted for the steps they are observed.

Because the cumulative degree grows with the length of the observation window, window should be set to a substantively meaningful period (for example, one year of time steps) when comparing across models or against data.

Value

A data.frame with one row per simulation and degree value, and the columns:

• sim: the simulation number.

• degree: the degree value, from 0 to the maximum observed across the selected simulations.

• count: the number of nodes with that degree. For type = "momentary", this is the mean number of nodes per time step, and so need not be a whole number.

• prop: count divided by the network size.

See Also

plot.netdx() with type = "cumldeg" to plot these distributions, and netdx() to generate the input object.

Examples

nw <- network_initialize(n = 100)
est <- netest(nw, formation = ~edges, target.stats = 50,
              coef.diss = dissolution_coefs(~offset(edges), duration = 10),
              verbose = FALSE)
dx <- netdx(est, nsims = 2, nsteps = 50, keep.tedgelist = TRUE,
            verbose = FALSE)

# Cumulative degree over the full simulation
cd <- get_degree_dist(dx, type = "cumulative")
head(cd)

# Mean cumulative degree per simulation
tapply(cd$degree * cd$prop, cd$sim, sum)

# Momentary degree over the last 25 steps
get_degree_dist(dx, type = "momentary", window = c(25, 50))

Description

This function outputs an edgelist from the specified network, selecting the method depending on the stored network type.

Usage

get_edgelist(dat, network)

Arguments

Value

An edgelist in matrix form with two columns. Each column contains the posit_ids (see get_posit_ids) of the nodes in each edge.

Description

Get the Edgelist(s) from the Specified Network(s)

Usage

get_edgelists_df(dat, networks = NULL)

Arguments

Value

A data.frame with the following columns:

• head: Positional ID of the head node.

• tail: Positional ID of the tail node.

• network: The numerical index of the network on which the edge is located.

Description

Retrieve, create, and edit the epidemic outputs of the main netsim_dat object of network models, from inside module functions. Each epidemic output is a vector with one value per time step, which becomes a column of the data frame returned by as.data.frame.netsim() after the simulation.

• get_epi() returns one output, for all time steps or for the time steps in at. get_epi_list() returns several outputs as a named list.

• set_epi() records value at one time step, the current one by default, and creates the output if it does not exist. add_epi() creates an output set to NA at every time step.

Unlike the other set_ functions, set_epi() takes a time step at, which comes last and must be passed by name.

See net-accessor for an overview of the accessors.

Usage

get_epi(dat, item, at = NULL, override.null.error = FALSE)

get_epi_list(dat, item = NULL)

set_epi(dat, item, value, ..., at = NULL)

add_epi(dat, item)

Arguments

Value

For get_epi(), the output vector, or its values at the time steps in at; for get_epi_list(), a named list of output vectors. The other functions return the updated netsim_dat object.

Intended Usage

These accessors operate on the live netsim_dat object that EpiModel constructs internally when a simulation starts and passes through each module at every time step. The expected calling context is inside a user-supplied module function (e.g., a custom infection.FUN, recovery.FUN, or arrivals/departures module) registered with control.net().

They are not a general-purpose tool for editing the user-facing input or output objects of netsim():

• Calling them on a param.net(), init.net(), or control.net() object before a run will appear to succeed (those objects are also list-like) but produces an object that is not a valid simulation input.

• Calling them on the object returned by netsim() modifies only the saved input slot, not anything that would change the result of a re-run. Note also that netsim() strips the ⁠*.FUN⁠ entries from its returned control slot, so feeding that slot back into a new simulation will fail.

For modifications outside a running simulation, use the following instead:

• Editing parameters before a run: update_params() for a param.net object, or direct list assignment (e.g., p$inf.prob <- 0.5).

• Editing init or control before a run: direct list assignment (e.g., ctrl$nsteps <- 1000), or rebuild with a fresh call to the init.net() / control.net() constructor.

• Scheduled mid-simulation changes: the .param.updater.list argument to param.net() and the .control.updater.list argument to control.net(). See the "model-parameters" vignette for the underlying updater module and the scenario API built on top of it.

Mutability

The set_, append_, and add_ functions DO NOT modify the netsim_dat object in place. The result must be assigned back to dat in order to be registered: ⁠dat <- set_*(dat, item, value)⁠.

See Also

net-accessor for an overview of the accessors, and get_attr() and get_param() for the other groups; epi_trackers() for outputs computed by tracker functions at each time step; mutate_epi() for outputs derived after the simulation.

Examples

dat <- create_dat_object(control = list(nsteps = 5))
dat <- set_current_timestep(dat, 1)

dat <- add_epi(dat, "i.num")
dat <- set_epi(dat, "i.num", 10)
dat <- set_epi(dat, "s.num", 90)
get_epi_list(dat)

# Record at a time step other than the current one
dat <- set_epi(dat, "i.num", 12, at = 2)
get_epi_list(dat, c("i.num", "s.num"))
get_epi(dat, "i.num")
get_epi(dat, "i.num", at = c(1, 2))

Description

Given a formation formula for a network model, outputs a character vector of vertex attributes to be used in netsim() simulations.

Usage

get_formula_term_attr(form, nw)

Arguments

Value

A character vector of vertex attributes.

Description

Extracts the network object from either a network epidemic model object generated with netsim, a network diagnostic simulation generated with netdx, or a netsim_dat object used internally in netsim. For netdx or netsim with tergmLite == FALSE, the extracted network object is a networkDynamic, which can be collapsed down to a static network object with the collapse and at arguments. For netsim with tergmLite == TRUE, the extracted network object is the final networkLite, the collapse argument should be FALSE, and the at argument should be missing. For netsim_dat, the collapse and at arguments are not supported, and the network object is either the current networkLite (if tergmLite == TRUE) or the current networkDynamic (if tergmLite == FALSE).

Usage

get_network(x, ...)

## S3 method for class 'netdx'
get_network(x, sim = 1, collapse = FALSE, at = NULL, ...)

## S3 method for class 'netsim'
get_network(x, sim = 1, network = 1, collapse = FALSE, at = NULL, ...)

## S3 method for class 'netsim_dat'
get_network(x, network = 1L, ...)

Arguments

Details

This function requires that the network object is saved during the network simulation while running either netsim() or netdx(). For the former, that is specified by setting the save.network parameter in control.net() to TRUE. For the latter, that is specified with the keep.tnetwork parameter directly in netdx().

Value

For netdx or netsim with tergmLite == FALSE, a networkDynamic object (if collapse = FALSE) or a static network object (if collapse = TRUE). For netsim with tergmLite == TRUE or netsim_dat with tergmLite == TRUE, a networkLite object. For netsim_dat with tergmLite == FALSE, a networkDynamic object.

Examples

# Set up network and TERGM formula
nw <- network_initialize(n = 100)
nw <- set_vertex_attribute(nw, "group", rep(1:2, each = 50))
formation <- ~edges
target.stats <- 50
coef.diss <- dissolution_coefs(dissolution = ~offset(edges), duration = 20)

# Estimate the model
est <- netest(nw, formation, target.stats, coef.diss)

# Run diagnostics, saving the networkDynamic objects
dx <- netdx(est, nsteps = 10, nsims = 3, keep.tnetwork = TRUE,
            verbose = FALSE)

# Extract the network for simulation 2 from dx object
get_network(dx, sim = 2)

# Extract and collapse the network from simulation 1 at time step 5
get_network(dx, collapse = TRUE, at = 5)

# Parameterize the epidemic model, and simulate it
param <- param.net(inf.prob = 0.3, inf.prob.g2 = 0.15)
init <- init.net(i.num = 10, i.num.g2 = 10)
control <- control.net(type = "SI", nsteps = 10, nsims = 3, verbose = FALSE)
mod <- netsim(est, param, init, control)

# Extract the network for simulation 2 from mod object
get_network(mod, sim = 2)

## Extract and collapse the network from simulation 1 at time step 5
get_network(mod, collapse = TRUE, at = 5)

Description

Gets all network attributes except "mnext" from its network argument.

Usage

get_network_attributes(x)

Arguments

Details

This function is used in EpiModel workflows to copy relevant network attributes from the network object to the netsim_dat object when initializing netsim runs.

Value

Returns the named list of network attributes.

Examples

nw <- network_initialize(100)
get_network_attributes(nw)

Description

Given a simulated network, outputs a character vector of vertex attributes to be used in netsim() simulations.

Usage

get_network_term_attr(nw)

Arguments

Value

A character vector of vertex attributes.

Description

Extracts network statistics from a network epidemic model simulated with netsim or a network diagnostics object simulated with netdx. Statistics can be returned either as a single data frame or as a list of matrices (one matrix for each simulation).

Usage

get_nwstats(x, sim = NULL, network = 1, mode = c("data.frame", "list"))

Arguments

Value

A data frame or list of matrices containing the network statistics.

Examples

# Two-group Bernoulli random graph TERGM
nw <- network_initialize(n = 100)
nw <- set_vertex_attribute(nw, "group", rep(1:2, each = 50))
formation <- ~edges
target.stats <- 50
coef.diss <- dissolution_coefs(dissolution = ~offset(edges), duration = 20)
est <- netest(nw, formation, target.stats, coef.diss, verbose = FALSE)

dx <- netdx(est, nsim = 3, nsteps = 10, verbose = FALSE,
            nwstats.formula = ~edges + isolates)
get_nwstats(dx)
get_nwstats(dx, sim = 1)

# SI epidemic model
param <- param.net(inf.prob = 0.3, inf.prob.g2 = 0.15)
init <- init.net(i.num = 10, i.num.g2 = 10)
control <- control.net(type = "SI", nsteps = 10, nsims = 3,
                       nwstats.formula = ~edges + meandeg + degree(0:5),
                       verbose = FALSE)
mod <- netsim(est, param, init, control)

# Extract the network statistics from all or sets of simulations
get_nwstats(mod)
get_nwstats(mod, sim = 2)
get_nwstats(mod, sim = c(1, 3))

# On the fly summary stats
summary(get_nwstats(mod))
colMeans(get_nwstats(mod))

Description

Retrieve, create, and edit the parameters, control settings, and initial conditions of the main netsim_dat object of network models, from inside module functions. They start as the param.net(), control.net(), and init.net() inputs of netsim(), and the three share one set of accessors:

• get_param(), get_control(), and get_init() return one element. get_param_list(), get_control_list(), and get_init_list() return several elements as a named list.

• set_param(), set_control(), and set_init() replace the value of an element, and create the element if it does not exist. add_param(), add_control(), and add_init() create an element set to NA.

• get_network_control() returns the entry for one network layer of a control setting given as a multilayer() object.

Unlike set_attr() and set_epi(), these set_ functions take neither nodes nor a time step: the new value is used from then on, until it is set again.

See net-accessor for an overview of the accessors.

Usage

get_param(dat, item, override.null.error = FALSE)

get_param_list(dat, item = NULL)

set_param(dat, item, value)

add_param(dat, item)

get_control(dat, item, override.null.error = FALSE)

get_control_list(dat, item = NULL)

get_network_control(dat, network, item, override.null.error = FALSE)

set_control(dat, item, value)

add_control(dat, item)

get_init(dat, item, override.null.error = FALSE)

get_init_list(dat, item = NULL)

set_init(dat, item, value)

add_init(dat, item)

Arguments

Value

For get_param(), get_control(), and get_init(), the value of the element; for get_network_control(), its entry for layer network; for the ⁠get_*_list()⁠ functions, a named list of elements. The other functions return the updated netsim_dat object.

Intended Usage

These accessors operate on the live netsim_dat object that EpiModel constructs internally when a simulation starts and passes through each module at every time step. The expected calling context is inside a user-supplied module function (e.g., a custom infection.FUN, recovery.FUN, or arrivals/departures module) registered with control.net().

They are not a general-purpose tool for editing the user-facing input or output objects of netsim():

• Calling them on a param.net(), init.net(), or control.net() object before a run will appear to succeed (those objects are also list-like) but produces an object that is not a valid simulation input.

• Calling them on the object returned by netsim() modifies only the saved input slot, not anything that would change the result of a re-run. Note also that netsim() strips the ⁠*.FUN⁠ entries from its returned control slot, so feeding that slot back into a new simulation will fail.

For modifications outside a running simulation, use the following instead:

• Editing parameters before a run: update_params() for a param.net object, or direct list assignment (e.g., p$inf.prob <- 0.5).

• Editing init or control before a run: direct list assignment (e.g., ctrl$nsteps <- 1000), or rebuild with a fresh call to the init.net() / control.net() constructor.

• Scheduled mid-simulation changes: the .param.updater.list argument to param.net() and the .control.updater.list argument to control.net(). See the "model-parameters" vignette for the underlying updater module and the scenario API built on top of it.

Mutability

The set_, append_, and add_ functions DO NOT modify the netsim_dat object in place. The result must be assigned back to dat in order to be registered: ⁠dat <- set_*(dat, item, value)⁠.

See Also

net-accessor for an overview of the accessors, and get_attr() and get_epi() for the other groups; update_params() for editing a param.net object outside a simulation.

Examples

dat <- create_dat_object(control = list(nsteps = 150))

dat <- add_param(dat, "x")
get_param(dat, "x")
dat <- set_param(dat, "x", 0.4)
dat <- set_param(dat, "y", 0.8)
get_param_list(dat)
get_param(dat, "x")

# The control and init accessors work the same way
dat <- add_control(dat, "x")
dat <- set_control(dat, "x", 0.4)
get_control_list(dat, c("nsteps", "x"))
get_control(dat, "nsteps")

dat <- add_init(dat, "x")
dat <- set_init(dat, "x", 0.4)
get_init_list(dat)
get_init(dat, "x")

# A control setting with one entry per network layer
dat <- set_control(dat, "z", multilayer(0.1, 0.2))
get_network_control(dat, network = 2, "z")

Description

Pulls every cumulative-edgelist row touching one of the supplied "index" nodes, returning each ⁠(index, partner)⁠ pair together with the partnership start/stop times and the network layer. This is the building block for contact tracing over the simulated network history.

Usage

get_partners(
  dat,
  index_posit_ids,
  networks = NULL,
  truncate = Inf,
  only.active.nodes = FALSE
)

Arguments

Details

Indexes are passed as positional IDs but the output uses unique IDs, because partners may include nodes that have already departed (and thus no longer have a positional ID). Use get_unique_ids() and get_posit_ids() to convert between the two systems.

The truncate argument here filters by edge age: only edges whose stop step is within truncate time steps of the current step are kept (active edges are always included). It operates on whatever history is already stored in the cumulative edgelist; it cannot recover edges that were dropped earlier by update_cumulative_edgelist() or by control$truncate.el.cuml.

Value

A data.frame with 5 columns:

• index: the unique IDs of the indexes.

• partner: the unique IDs of the partners/contacts.

• start: the time step at which the edge started.

• stop: the time step in which the edge stopped; if ongoing, then NA is returned.

• network: the numerical index for the network on which the partnership/contact is located.

See Also

get_cumulative_degree() for a partner count per index; get_cumulative_edgelist() / get_cumulative_edgelists_df() for the underlying edgelist; get_unique_ids() / get_posit_ids() for ID conversion.

Other cumulative_edgelist: as_cumulative_edgelist(), dedup_cumulative_edgelist(), get_cumulative_degree(), get_cumulative_edgelist(), get_cumulative_edgelists_df(), reachable-nodes, update_cumulative_edgelist()

Examples

## Not run: 
# Contacts of nodes 1..10 across all networks within the last 30 steps,
# excluding partners who have since departed:
get_partners(dat,
             index_posit_ids = 1:10,
             truncate = 30,
             only.active.nodes = TRUE)

## End(Not run)

Description

Subsets the entire netsim object to a subset of simulations, essentially functioning like a reverse of merge.

Usage

get_sims(x, sims, var = NULL)

Arguments

Value

An updated object of class netsim containing only the simulations specified in sims and the variables specified in var.

Examples

# Network model estimation
nw <- network_initialize(n = 100)
formation <- ~edges
target.stats <- 50
coef.diss <- dissolution_coefs(dissolution = ~offset(edges), duration = 20)
est1 <- netest(nw, formation, target.stats, coef.diss, verbose = FALSE)

# Epidemic model
param <- param.net(inf.prob = 0.3)
init <- init.net(i.num = 10)
control <- control.net(type = "SI", nsteps = 10, nsims = 3, verbose.int = 0)
mod1 <- netsim(est1, param, init, control)

# Get sim 2
s.g2 <- get_sims(mod1, sims = 2)

# Get sims 2 and 3 and keep only a subset of variables
s.g2.small <- get_sims(mod1, sims = 2:3, var = c("i.num", "si.flow"))

# Extract the mean simulation for the variable i.num
sim.mean <- get_sims(mod1, sims = "mean", var = "i.num")

Description

Return an adjacency list of subnets

Usage

get_subnet_adj_list(adj_list)

Arguments

Details

A graph with 4 components: 1, 2, 3, 4, and 5 and 6, 7, 8 would yield a list like so: 1: 2, 3, 4 2: 1 3: 1 4: 1 5: numeric(0) 6: 7, 8 7: 6, 8: 6

This format speeds up the construction of reachable sets on dense networks

Value

An adjacency list where only the first node of a subnet contains the subnet and all other contain only the first node

Description

Gets a vertex attribute from an object of class network. This function simplifies the related function in the network package.

Usage

get_vertex_attribute(x, attrname)

Arguments