Skip to contents

Simulates stochastic network epidemic models for infectious disease. This is typically the final step in the network modeling pipeline, after network estimation with netest and model diagnostics with netdx.

Usage

netsim(x, param, init, control)

Arguments

x

If control$start == 1, either a fitted network model object of class netest, a clique layer of class netclique, an observed network layer of class netcensus, or a list of such objects. If control$start > 1, an object of class netsim. When multiple networks are used (multi-layer models), pass a list with one netest, netclique, or netcensus object per network layer, in any order; the node sets (including network size and nodal attributes) are assumed to be the same for all networks.

param

Model parameters, as an object of class param.net. Includes transmission probability (inf.prob), act rate (act.rate), recovery rate (rec.rate), and demographic rates for models with vital dynamics. Custom parameters for extended models may also be passed through param.net.

init

Initial conditions, as an object of class init.net. Specifies the initial number of infected nodes (i.num) and, for SIR models, recovered nodes (r.num). For two-group models, the corresponding .g2 parameters are also required.

control

Control settings, as an object of class control.net. Key settings include type (disease model: "SI", "SIR", or "SIS"), nsteps (number of time steps), nsims (number of simulations), tergmLite (lightweight mode for performance), and resimulate.network (required for models with vital dynamics). For extended models, custom module functions are also passed here.

Value

A list of class netsim 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.

  • stats: a list containing two sublists, nwstats for any network statistics saved in the simulation, and transmat for the transmission matrix saved in the simulation. See control.net() for further details.

  • network: a list of lists of networkDynamic or networkLite objects, with one list of objects for each model simulation.

If control$raw.output == TRUE: A list of the raw (pre-processed) netsim_dat objects, for use in simulation continuation.

The epi data can be extracted as a data frame with as.data.frame.netsim(), with options for per-simulation values (out = "vals"), means (out = "mean"), standard deviations (out = "sd"), or quantiles (out = "qnt").

Details

Stochastic network models explicitly represent phenomena within and across edges (pairs of nodes that remain connected) over time. This enables edges to have duration, allowing for repeated transmission-related acts within the same dyad, specification of edge formation and dissolution rates, control over the temporal sequencing of multiple edges, and specification of network-level features. A detailed description of these models, along with examples, is found in the Network Modeling for Epidemics course materials.

The netsim function performs modeling of both the base model types and original models. 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).

Original models may be parameterized by writing new process modules that either take the place of existing modules (for example, disease recovery), or supplement the set of existing processes with a new one contained in a new module. This functionality is documented in the Extending EpiModel section of the Network Modeling for Epidemics course materials. The list of modules within netsim available for modification is listed in modules.net().

Module Pipeline

At each time step, netsim executes a series of modules in sequence. For base models, the default pipeline is:

  1. Network resimulation (resim_nets): updates the network structure (if resimulate.network = TRUE).

  2. Infection (infection.net): simulates disease transmission across discordant edges (where one partner is susceptible and the other is infected).

  3. Recovery (recovery.net): simulates recovery from infection (SIR and SIS models only).

  4. Departures (departures.net): simulates node exits from the network (if vital dynamics are enabled).

  5. Arrivals (arrivals.net): simulates new node entries into the network (if vital dynamics are enabled).

  6. Prevalence (prevalence.net): calculates summary statistics.

See modules.net() for full details on each module.

Performance and tergmLite

Setting tergmLite = TRUE in control.net uses a lightweight network representation (networkLite) that is substantially faster than the full networkDynamic representation, often by a factor of 20–50x. This is recommended for large networks or when running many simulations. The trade-off is that full networkDynamic objects are not available for post-hoc analysis; use the cumulative edgelist (cumulative.edgelist = TRUE) instead to track partnership histories.

The resimulate.network control must be set to TRUE when demographic processes (arrivals and departures) change the network composition over time. Without it, the network structure evolves independently of the epidemic and demographic dynamics.

Partnership Histories (Cumulative Edgelist)

For contact tracing, partnership-duration analysis, or reachability analysis, enable cumulative-edgelist tracking via control.net(cumulative.edgelist = TRUE, save.cumulative.edgelist = TRUE). This records every edge formed during the simulation along with its start and stop time steps, and attaches the result to the returned object as sim$cumulative.edgelist. Under tergmLite = TRUE this is the recommended substitute for the full networkDynamic history.

The helper family (get_cumulative_edgelist, get_cumulative_edgelists_df, get_partners, get_cumulative_degree, get_forward_reachable, get_backward_reachable) operates on the tracked data. See vignette("network-objects", package = "EpiModel") for the full lifecycle.

Multi-Network Models

For models with multiple overlapping network layers (e.g., sexual and needle-sharing networks), pass a list of netest objects to the x argument, one per network layer. Each layer has its own formation/dissolution dynamics but shares the same node set. See the multilayer documentation for specifying controls that vary by layer, and the Multi-Layer Networks chapter of the Network Modeling for Epidemics course materials for a worked example building from independent to cross-dependent layers.

A layer made of groups that are cliques by definition, such as households, classrooms, wards, or cabins, may instead be built with netclique and placed in the same list. This is an addition to netest, not a replacement for it. An ERGM remains the natural model for a network whose ties depend on nodal and dyadic predictors (degree, mixing by attribute, clustering), whether those ties turn over quickly, slowly, or not at all; the model is the compact description of that structure. In a clique layer every pair of group members is a contact and there is no tie-formation process to estimate. An ERGM can reproduce such cliques only through a nodematch term targeted at its maximum, where no finite coefficient exists, so the fit is unstable and rarely completes every group. A clique layer is skipped by the network resimulation and the edges correction, its edges are removed with departing nodes, and arriving nodes are placed on it under the rule chosen in netclique.

A fully observed network, a static network or a networkDynamic with edge spells from sensors, proximity loggers, or contact tracing, may be placed in the list with netcensus. Here estimation is skipped because the whole network is in hand rather than a sample to generalize from; the simulation reads the edges active at each time step from the observed object, and the node set is fixed, so vital dynamics are refused.

The built-in infection modules read every layer alike; to give each layer its own transmission probability or act rate, pass inf.prob and act.rate as multilayer objects in param.net. The transmission matrix records the layer of every transmission in its network column.

Restarting and Checkpointing

Simulations can be checkpointed and restarted if interrupted. Set .checkpoint.steps and .checkpoint.dir in control.net to enable automatic checkpointing. To restart a simulation from a prior netsim output, pass the netsim object as x and set control$start to one greater than the final time step of the prior simulation. See the Checkpointing Simulations section of control.net for full details. The epi trackers and network statistics of the output hold every time step from 1 to nsteps, those of x followed by the new ones, and as.data.frame.netsim and plot.netsim label them 1 to nsteps.

When restarting, every simulation restarts from the first simulation of x, with a warning if x holds more than one. make_restart_point makes a lightweight x holding the single simulation to restart from (get_sims selects one from a netsim object). The parameters are taken from param alone, which must hold all the parameters of the model: the ones saved in x are not used.

References

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.

See also

Estimate the network model with netest and diagnose model fit with netdx before running simulations. Extract model results with as.data.frame.netsim(). Summarize the time-specific model results with summary.netsim(). Plot the model results with plot.netsim(). Extract the network with get_network, the transmission matrix with get_transmat, and derive new epi statistics with mutate_epi.

Examples

if (FALSE) { # \dontrun{
## Example 1: SI Model without Network Feedback
# 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 = 100, nsims = 5, verbose.int = 0)
mod1 <- netsim(est1, param, init, control)

# Print, plot, and summarize the results
mod1
plot(mod1)
summary(mod1, at = 50)

## Example 2: SIR Model with Network Feedback (Two-Group)
nw <- network_initialize(n = 100)
nw <- set_vertex_attribute(nw, "group", rep(1:2, each = 50))
formation <- ~edges
target.stats <- 50

# Recalculate dissolution coefficient with departure rate
coef.diss <- dissolution_coefs(dissolution = ~offset(edges), duration = 20,
                               d.rate = 0.0021)

# Reestimate the model with new coefficient
est2 <- netest(nw, formation, target.stats, coef.diss)

# Set parameters to include demographic rates
param <- param.net(inf.prob = 0.3, inf.prob.g2 = 0.15,
                   rec.rate = 0.02, rec.rate.g2 = 0.02,
                   a.rate = 0.002, a.rate.g2 = NA,
                   ds.rate = 0.001, ds.rate.g2 = 0.001,
                   di.rate = 0.001, di.rate.g2 = 0.001,
                   dr.rate = 0.001, dr.rate.g2 = 0.001)
init <- init.net(i.num = 10, i.num.g2 = 10,
                 r.num = 0, r.num.g2 = 0)
control <- control.net(type = "SIR", nsteps = 100, nsims = 5,
                       resimulate.network = TRUE, tergmLite = TRUE)

# Simulate the model with new network fit
mod2 <- netsim(est2, param, init, control)

# Print, plot, and summarize the results
mod2
plot(mod2)
summary(mod2, at = 40)

## Example 3: Post-Simulation Analysis
# Extract epi data as a data frame
as.data.frame(mod1)
as.data.frame(mod1, out = "mean")

# Extract the transmission matrix from simulation 1
get_transmat(mod1, sim = 1)

# Derive new epi statistics and plot them
mod1 <- mutate_epi(mod1, prev = i.num / num)
plot(mod1, y = "prev")
} # }