21 SIS Epidemic Across a Dynamic Network
We build and simulate a network-based epidemic model in EpiModel. EpiModel uses temporal exponential-family random graph models (TERGMs) to estimate and simulate the whole (sociocentric) network of a population, based on individual-level, dyad-level, and network-level patterns of density, degree, assortativity, and other features influencing edge formation and dissolution. This is a multi-step process: it begins with estimation of a temporal ERGM and continues with simulation of a dynamic network and epidemic processes on top of that network.
In this tutorial, we work through a model of a Susceptible-Infected-Susceptible (SIS) epidemic. One example of an SIS disease would be a bacterial sexually transmitted infection such as Gonorrhea, in which individuals may acquire infection from sexual contact with an infected partner, and then recover from infection either through natural clearance or through antibiotic treatment.
Network modeling is not limited to STIs. The same formation and dissolution machinery applies to directly transmitted infections, such as a respiratory pathogen spreading over household, workplace, or school contact networks, where the edges represent the relevant person-to-person contacts rather than sexual partnerships.
We will use a simplifying assumption of a closed population, in which there are no entries or exits from the network; this may be justified by the short time span over which the epidemic will be simulated.
Download the R script to follow along with this tutorial here.
21.1 Network Model Estimation
To get started, load the EpiModel library.
The first step in our network model is to specify a network structure, including features like size and nodal attributes. The network_initialize function creates an object of class network. Below we show an example of initializing a network of 500 nodes, with no edges between them at the start. Edges represent sexual partnerships (mutual person-to-person contact), so this is an undirected network.
The sizes of the networks represented in this workshop are smaller than what might be used for a research-level model, mostly for computational efficiency. Larger network sizes over longer time intervals are typically used for research purposes.
21.1.1 Model Parameterization
This example will start simple, with a formula that represents the network density and the level of concurrency (overlapping sexual partnerships) in the population. This is a dyad-dependent ERGM, since the probability of edge formation between any two nodes depends on the existence of edges between those nodes and other nodes. The concurrent term is defined as the number of nodes with at least two partners at any time. Following the notation of the tergm package, we specify this using a right-hand side (RHS) formula. In addition to concurrency, we will use a constraint on the degree distribution. This will cap the degree of any person at 3, with no nodes allowed to have 4 or more ongoing partnerships. This type of constraint could reflect a truncated sampling scheme for partnerships within a survey (e.g., respondents only asked about their 3 most recent partners), or a model assumption about limits of human activity.
Target statistics will be the input mechanism for formation model terms. The edges term will be a function of mean degree, or the average number of ongoing partnerships. With an arbitrarily specified mean degree of 0.7, the corresponding target statistic is 175: \(edges = mean \ degree \times \frac{N}{2}\).
We will also specify that 22% of individuals exhibit concurrency, so the corresponding target statistic is 110: \(concurrent = 0.22 \times N = 0.22 \times 500\). An edges-only model implies an approximately Poisson degree distribution, so at a mean degree of 0.7 the probability of two or more partners is 1 - dpois(0, 0.7) - dpois(1, 0.7), about 16%; the assumed 22% therefore sits slightly above chance. The target statistic for the number of individuals with a momentary degree of 4 or more is 0, reflecting our assumed constraint. The three values in target.stats correspond, in order, to the three formation terms: edges = 175, concurrent = 110, and degrange(from = 4) = 0.
The dissolution model is parameterized from a mean partnership duration estimated from cross-sectional egocentric data. It differs from the formation model in two respects.
First, the dissolution model is not estimated in the ERGM. It is instead passed in as a fixed coefficient, conditional on which the formation model is estimated. The dissolution model terms are calculated analytically using the dissolution_coefs function, the output of which is passed into the netest model estimation function.
Second, whereas formation models may be arbitrarily complex, dissolution models are limited to a set of dyad-independent models. Dyad-independent means that a partnership’s probability of dissolution does not depend on the other partnerships in the network. These supported models are listed in the dissolution_coefs function help page.
The model we will use is an edges-only model, implying a homogeneous probability of dissolution for all partnerships in the network. The average duration of these partnerships will be specified at 50 time steps, which will be days in our model.
A time step is arbitrary, representing whatever unit we choose when we set durations and rates; we treat a step as a day here, while other Module 4 tutorials use weeks, purely for illustration.
Dissolution Coefficients
=======================
Dissolution Model: ~offset(edges)
Target Statistics: 50
Crude Coefficient: 3.89182
Mortality/Exit Rate: 0
Adjusted Coefficient: 3.89182The output from this function indicates both an adjusted and crude coefficient. They are equivalent here because the adjustment applies only in an open population with entries and exits, where the coefficient must also account for partnerships that end when a partner leaves the population (through death or departure) rather than through dissolution itself. Upcoming workshop material will showcase when they differ as a result of exits from the network.
21.1.2 Model Estimation and Diagnostics
In EpiModel, network model estimation is performed with the netest function, which is a wrapper around the estimation functions in the ergm and tergm packages. The function arguments are as follows:
function (nw, formation, target.stats, coef.diss, constraints = NULL,
coef.form = NULL, edapprox = TRUE, set.control.ergm = control.ergm(),
set.control.tergm = control.tergm(MCMC.maxchanges = .Machine$integer.max),
set.control.ergm.ego = NULL, verbose = FALSE, nested.edapprox = TRUE,
...)
NULLThe four arguments that must be specified with each function call are:
nw: an initialized empty network.formation: a RHS formation formula.target.stats: target statistics for the formation model.coef.diss: output object fromdissolution_coefs, containing the dissolution coefficients.
Other arguments that may be helpful to understand when getting started are:
constraints: this is another way of inputting model constraints (seehelp("ergm")).coef.form: sets the coefficient values of any offset terms in the formation model (those that are not explicitly estimated but fixed).edapprox: selects the dynamic estimation method. IfTRUE(the default), uses the approximation method; ifFALSE, the direct method.- Direct method: uses the functionality of the
tergmpackage to estimate the separable formation and dissolution models for the network. This is often not used because of computational time. - Approximation method: uses
ergmestimation for a cross-sectional network (the prevalence of edges) with an analytic adjustment of the edges coefficient to account for dissolution (i.e., transformation from prevalence to incidence). This approximation method may introduce bias into estimation in certain cases (high density and short durations) but these are typically not a concern for the low density cases in epidemiologically relevant networks.
- Direct method: uses the functionality of the
21.1.2.1 Estimation
Because we have a dyad-dependent model, MCMC will be used to estimate the coefficients of the model given the target statistics. We leave edapprox at its default, so this fit uses the edges dissolution approximation (edapprox = TRUE), which is the method the diagnostics below assess.
21.1.2.2 Diagnostics
There are two forms of model diagnostics for a dynamic ERGM fit with netest: static and dynamic diagnostics. When the approximation method has been used, static diagnostics check the fit of the cross-sectional model to target statistics. Dynamic diagnostics check the fit of the model adjusted to account for edge dissolution.
When running a dynamic network simulation, it is good to start with the dynamic diagnostics, and if there are fit problems, work back to the static diagnostics to determine if the problem is due to the cross-sectional fit itself or with the dynamic adjustment (i.e., the approximation method). A proper fitting ERGM using the approximation method does not guarantee well-performing dynamic simulations.
Here we will examine dynamic diagnostics only. These are run with the netdx function, which simulates from the model fit object returned by netest. One must specify the number of simulations from the dynamic model and the number of time steps per simulation. Choice of both simulation parameters depends on the stochasticity in the model, which is a function of network size, model complexity, and other factors. The nwstats.formula contains the network statistics to monitor in the diagnostics: it may contain statistics in the formation model and also others. By default, it is the formation model. Finally, we are keeping the “timed edgelist” (a record of every partnership with its start and end time) with keep.tedgelist.
We have also built parallelization into the EpiModel simulation functions, so it is also possible to run multiple simulations at the same time using your computer’s multi-core design. You can find the number of cores in your system with:
Then you can run the multi-core simulations by specifying ncores (EpiModel will prevent you from specifying more cores than you have available).
Printing the object will show the object structure and diagnostics. Both formation and duration diagnostics show a good fit relative to their targets. For the formation diagnostics, the mean statistics are the mean of the cross sectional statistics at each time step across all simulations. The Pct Diff column shows the relative difference between the mean and targets. There are two forms of dissolution diagnostics. The edge duration row shows the mean duration of partnerships across the simulations; calculating this involves some imputation due to the length censoring at the start of the simulation. The next row shows the percent of current edges dissolving at each time step; this can be less intuitive than duration, but it does not require the imputation. The percentage of edges dissolving is the inverse of the expected duration: if the duration is 50 days, then we expect that 1/50 (or 2%) to dissolve each day.
EpiModel Network Diagnostics
=======================
Diagnostic Method: Dynamic
Simulations: 10
Time Steps per Sim: 1000
Formation Diagnostics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 175 176.351 0.772 1.174 1.151 5.120 13.484
meandeg NA 0.705 NA 0.005 NA 0.020 0.054
degree0 NA 270.493 NA 1.348 NA 5.120 14.077
degree1 NA 119.841 NA 0.452 NA 1.644 9.220
degree2 NA 96.137 NA 0.671 NA 3.457 11.033
degree3 NA 13.529 NA 0.225 NA 1.038 4.029
degree4 NA 0.000 NA NaN NA 0.000 0.000
concurrent 110 109.666 -0.304 0.834 -0.401 4.318 12.378
Duration Diagnostics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 50 50.658 1.315 0.325 2.021 1.211 3.838
Dissolution Diagnostics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 0.02 0.02 -0.571 0 -1.092 0 0.011Plotting the diagnostics object will show the time series of the target statistics against any targets. The other options used here specify to smooth the mean lines, give them a thicker line width, and plot each statistic in a separate panel. The black dashed lines show the value of the target statistics for any terms in the model. Similar to the numeric summaries, the plots show a good fit over the time series.
The simulated network statistics from diagnostic object may be extracted into a data.frame with get_nwstats.
time sim edges meandeg degree0 degree1 degree2 degree3 degree4 concurrent
1 1 1 186 0.744 265 112 109 14 0 123
2 2 1 183 0.732 268 113 104 15 0 119
3 3 1 179 0.716 270 115 102 13 0 115
4 4 1 180 0.720 271 112 103 14 0 117
5 5 1 182 0.728 268 116 100 16 0 116
6 6 1 177 0.708 274 113 98 15 0 113
7 7 1 182 0.728 272 107 106 15 0 121
8 8 1 180 0.720 273 107 107 13 0 120
9 9 1 174 0.696 276 110 104 10 0 114
10 10 1 176 0.704 271 118 99 12 0 111
11 11 1 177 0.708 272 114 102 12 0 114
12 12 1 177 0.708 274 110 104 12 0 116
13 13 1 184 0.736 269 107 111 13 0 124
14 14 1 178 0.712 274 106 110 10 0 120
15 15 1 178 0.712 273 109 107 11 0 118
16 16 1 181 0.724 273 103 113 11 0 124
17 17 1 180 0.720 274 102 114 10 0 124
18 18 1 175 0.700 278 106 104 12 0 116
19 19 1 175 0.700 276 108 106 10 0 116
20 20 1 178 0.712 274 106 110 10 0 120The dissolution model fit may also be assessed with plots by specifying either the duration or dissolution type, as defined above. The duration diagnostic is based on the average age of edges at each time step, up to that time step. An imputation algorithm is used for left-censored edges (i.e., those that exist at t1); you can turn off this imputation to see the effects of censoring with duration.imputed = FALSE. Both metrics show a good fit of the dissolution model to the target duration of 50 time steps.
By inspecting the timed edgelist, we can see the burn-in period directly with censoring of onset times. The as.data.frame function is used to extract this edgelist object.
onset terminus tail head onset.censored terminus.censored duration edge.id
1 0 1 3 55 FALSE FALSE 1 1
2 0 156 3 484 FALSE FALSE 156 2
3 0 252 7 404 FALSE FALSE 252 3
4 0 73 7 420 FALSE FALSE 73 4
5 0 2 9 245 FALSE FALSE 2 5
6 0 54 9 298 FALSE FALSE 54 6
7 0 41 10 102 FALSE FALSE 41 7
8 0 117 13 59 FALSE FALSE 117 8
9 0 21 13 168 FALSE FALSE 21 9
10 0 91 15 57 FALSE FALSE 91 10
11 0 19 18 267 FALSE FALSE 19 11
12 0 92 18 326 FALSE FALSE 92 12
13 0 78 18 380 FALSE FALSE 78 13
14 0 115 25 366 FALSE FALSE 115 14
15 0 9 27 386 FALSE FALSE 9 15
16 0 198 29 336 FALSE FALSE 198 16
17 0 24 35 196 FALSE FALSE 24 17
18 0 87 35 211 FALSE FALSE 87 18
19 0 7 38 400 FALSE FALSE 7 19
20 0 48 40 193 FALSE FALSE 48 20If the model diagnostics had suggested a poor fit, then additional diagnostics and fitting would be necessary, especially the cross-sectional diagnostics (setting dynamic to FALSE in netdx). Note that the number of simulations may be very large here and there are no time steps specified because each simulation is a cross-sectional network.
The plots now represent individual simulations from an MCMC chain, rather than time steps.
This lack of temporality is now evident when looking at the raw data.
sim edges concurrent deg4+
1 1 173 103 0
2 2 161 102 0
3 3 177 120 0
4 4 169 105 0
5 5 169 106 0
6 6 173 115 0
7 7 177 115 0
8 8 179 112 0
9 9 177 111 0
10 10 173 119 0
11 11 172 112 0
12 12 188 126 0
13 13 172 103 0
14 14 190 128 0
15 15 175 107 0
16 16 170 101 0
17 17 166 99 0
18 18 168 107 0
19 19 160 100 0
20 20 189 124 0If the cross-sectional model fits well but the dynamic model does not, then a full STERGM estimation may be necessary (using edapprox = FALSE). If the cross-sectional model does not fit well, different control parameters for the ERGM estimation may be necessary (see the help file for netdx for instructions).
21.2 Epidemic Simulation
EpiModel simulates disease epidemics over dynamic networks by integrating dynamic model simulations with the simulation of other epidemiological processes such as disease transmission and recovery. Like the network model simulations, these processes are also simulated stochastically so that the range of potential outcomes under the model specifications is estimated.
The specification of epidemiological processes to model may be arbitrarily complex, but EpiModel includes a number of “built-in” model types within the software. Additional components will be programmed and plugged into the simulation API (just like any epidemic model); we will introduce this later, and cover this in depth in our advanced workshop. Here, we will start simple with an SIS epidemic using this built-in functionality. This is just the starting point to what you can do in EpiModel!
21.2.1 Epidemic Model Parameters
Our SIS model will rely on three parameters. The act rate is the number of sexual acts that occur within a partnership each time unit. The overall frequency of acts per person per unit time is a function of the number of ongoing partnerships (mean degree) and this act rate parameter. The infection probability is the risk of transmission given contact with an infected person. The recovery rate for an SIS epidemic is the speed at which infected individuals become susceptible again. For a bacterial STI like gonorrhea, this may be a function of biological attributes like sex or use of therapeutic agents like antibiotics.
EpiModel uses three helper functions to input epidemic parameters, initial conditions, and other control settings for the epidemic model. Each function holds a distinct part of the specification: param.net holds the epidemic-process parameters (transmission probability, act rate, recovery rate), init.net holds the initial conditions (how many people start infected), and control.net holds the simulation settings and engine (model type, number of simulations, number of time steps, and any added modules). First, we use the param.net function to input the per-act transmission probability in inf.prob and the number of acts per partnership per unit time in act.rate. The recovery rate implies that the average duration of disease is 10 days (1/rec.rate).
For initial conditions in this model, we only need to specify the number of infected individuals at the outset of the epidemic. The remaining individuals in the network will be classified as disease susceptible.
The control settings specify the structural elements of the model. These include the disease type, number of simulations, and number of time steps per simulation. nsims is the number of stochastic replicates; run enough of them to characterize the run-to-run variability in outcomes (small or noisy models need more). nsteps is the number of time steps; set it long enough to cover the time horizon of interest, here for the SIS epidemic to reach its endemic equilibrium. (Here again we could use the model multi-core functionality by specifying an ncores value, but these models run so quickly that it’s not necessary.)
21.2.2 Simulating the Epidemic Model
Once the model has been parameterized, simulating the model is straightforward. One must pass the fitted network model object from netest along with the parameters, initial conditions, and control settings to the netsim function. With a no-feedback model like this (i.e., there are no vital dynamics parameters), the full dynamic network time series is simulated at the start of each epidemic simulation, and then the epidemiological processes are simulated over that structure.
Printing the model output lists the inputs and outputs of the model. The output includes the sizes of the compartments (s.num is the number susceptible and i.num is the number infected) and flows (si.flow is the number of infections and is.flow is the number of recoveries). Methods for extracting this output are discussed below.
EpiModel Simulation
=======================
Model class: netsim
Simulation Summary
-----------------------
Model type: SIS
No. simulations: 5
No. time steps: 500
No. NW groups: 1
Fixed Parameters
---------------------------
inf.prob = 0.4
act.rate = 2
rec.rate = 0.1
groups = 1
Model Output
-----------------------
Variables: s.num i.num num si.flow is.flow
Networks: sim1 ... sim5
Transmissions: sim1 ... sim5
Formation Statistics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 175 177.926 1.672 1.964 1.489 7.442 12.812
concurrent 110 110.332 0.302 1.301 0.255 6.308 11.594
deg4+ 0 0.000 NaN NaN NaN 0.000 0.000
Duration Statistics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 50 49.396 -1.209 0.482 -1.254 1.386 3.404
Dissolution Statistics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 0.02 0.02 0.285 0 0.269 0 0.01121.2.3 Model Analysis
Now that the model has been simulated, the next step is to analyze the data. This includes plotting the epidemiological output, the networks over time, and extracting other raw data.
21.2.3.1 Epidemic Plots
Plotting the output from the epidemic model using the default arguments will display the size of the compartments in the model across simulations. The means across simulations at each time step are plotted with lines, and the polygon band shows the inter-quartile range across simulations.
Graphical elements may be toggled on and off. The popfrac argument specifies whether to use the absolute size of compartments versus proportions.
Code
Whereas the default will print the compartment sizes, other elements of the simulation may be plotted by name with the y argument. Here we plot both flow sizes using smoothed means, which converge at model equilibrium by the end of the time series.
21.2.3.2 Network Plots
Another available plot type is a network plot to visualize the individual nodes and edges at a specific time point. Network plots are output by setting the type parameter to "network". To plot the disease infection status on the nodes, use the col.status argument: blue indicates susceptible and red infected. It is necessary to specify both a time step and a simulation number to plot these networks.
21.2.3.3 Time-Specific Model Summaries
The summary function with the output of netsim will show the model statistics at a specific time step. Here we output the statistics at the final time step, where roughly two-thirds of the population are infected.
EpiModel Summary
=======================
Model class: netsim
Simulation Details
-----------------------
Model type: SIS
No. simulations: 5
No. time steps: 500
No. NW groups: 1
Model Statistics
------------------------------
Time: 500
------------------------------
mean sd pct
Suscept. 326.8 16.208 0.654
Infect. 173.2 16.208 0.346
Total 500.0 0.000 1.000
S -> I 19.4 2.702 NA
I -> S 19.6 1.949 NA
------------------------------ 21.2.3.4 Data Extraction
The as.data.frame function may be used to extract the model output into a data frame object for easy analysis outside of the built-in EpiModel functions. The function default will output the raw data for all simulations for each time step.
sim time s.num i.num num si.flow is.flow
1 1 1 490 10 500 NA NA
2 1 2 489 11 500 6 5
3 1 3 485 15 500 5 1
4 1 4 482 18 500 3 0
5 1 5 481 19 500 3 2
6 1 6 481 19 500 4 4
7 1 7 475 25 500 8 2
8 1 8 469 31 500 6 0
9 1 9 471 29 500 2 4
10 1 10 469 31 500 6 4 sim time s.num i.num num si.flow is.flow
2491 5 491 318 182 500 20 17
2492 5 492 318 182 500 18 18
2493 5 493 319 181 500 20 21
2494 5 494 322 178 500 20 23
2495 5 495 319 181 500 22 19
2496 5 496 316 184 500 20 17
2497 5 497 327 173 500 12 23
2498 5 498 325 175 500 19 17
2499 5 499 327 173 500 17 19
2500 5 500 328 172 500 16 17Notice that the output above shows all compartment and flow sizes as integers, reinforcing this as an individual-level model.
The out argument may be changed to specify the output of means across the models (with out = "mean").
time s.num i.num num si.flow is.flow
1 1 490.0 10.0 500 NaN NaN
2 2 486.6 13.4 500 4.4 1.0
3 3 485.8 14.2 500 2.2 1.4
4 4 483.8 16.2 500 3.0 1.0
5 5 483.4 16.6 500 2.4 2.0
6 6 482.4 17.6 500 3.2 2.2
7 7 480.8 19.2 500 3.4 1.8
8 8 479.8 20.2 500 2.6 1.6
9 9 479.0 21.0 500 3.2 2.4
10 10 477.6 22.4 500 3.2 1.8 time s.num i.num num si.flow is.flow
491 491 318.6 181.4 500 20.4 19.2
492 492 319.8 180.2 500 16.8 18.0
493 493 321.8 178.2 500 18.6 20.6
494 494 318.4 181.6 500 20.2 16.8
495 495 319.0 181.0 500 18.8 19.4
496 496 320.8 179.2 500 17.8 19.6
497 497 323.6 176.4 500 18.2 21.0
498 498 323.2 176.8 500 19.8 19.4
499 499 326.6 173.4 500 17.6 21.0
500 500 326.8 173.2 500 19.4 19.6The networkDynamic objects are stored in the netsim object, and may be extracted with the get_network function. By default the dynamic networks are saved, and contain the full edge history for every node that has existed in the network, along with the disease status history of those nodes.
NetworkDynamic properties:
distinct change times: 502
maximal time range: 0 until Inf
Dynamic (TEA) attributes:
Vertex TEAs: testatus.active
Includes optional net.obs.period attribute:
Network observation period info:
Number of observation spells: 2
Maximal time range observed: 0 until 501
Temporal mode: discrete
Time unit: step
Suggested time increment: 1
Network attributes:
vertices = 500
directed = FALSE
hyper = FALSE
loops = FALSE
multiple = FALSE
bipartite = FALSE
net.obs.period: (not shown)
vertex.pid = tergm_pid
total edges= 1915
missing edges= 0
non-missing edges= 1915
Vertex attribute names:
active status tergm_pid testatus.active vertex.names
Edge attribute names not shown One thing you can do with that network dynamic object is to extract the timed edgelist of all ties that existed for that simulation.
onset terminus tail head onset.censored terminus.censored duration edge.id
1 0 18 2 10 FALSE FALSE 18 1
2 0 10 4 61 FALSE FALSE 10 2
3 0 140 7 280 FALSE FALSE 140 3
4 0 48 7 352 FALSE FALSE 48 4
5 0 44 7 450 FALSE FALSE 44 5
6 0 14 9 353 FALSE FALSE 14 6
7 0 41 12 22 FALSE FALSE 41 7
8 0 82 12 70 FALSE FALSE 82 8
9 0 150 14 268 FALSE FALSE 150 9
10 0 23 14 403 FALSE FALSE 23 10
11 0 2 15 167 FALSE FALSE 2 11
12 0 90 15 235 FALSE FALSE 90 12
13 0 129 17 265 FALSE FALSE 129 13
14 0 102 19 472 FALSE FALSE 102 14
15 0 3 19 475 FALSE FALSE 3 15
16 0 10 20 281 FALSE FALSE 10 16
17 0 17 23 157 FALSE FALSE 17 17
18 0 81 23 361 FALSE FALSE 81 18
19 0 23 24 380 FALSE FALSE 23 19
20 0 48 25 226 FALSE FALSE 48 20
21 0 1 26 289 FALSE FALSE 1 21
22 0 18 26 385 FALSE FALSE 18 22
23 0 20 27 248 FALSE FALSE 20 23
24 0 54 27 493 FALSE FALSE 54 24
25 0 36 28 193 FALSE FALSE 36 25One can also use the get_transmat function generate a record of some key details about each transmission event that occurred. Shown below are the first 10 transmission events for simulation number 1. The sus column shows the unique ID of the previously susceptible, newly infected node in the event. The inf column shows the ID of the transmitting node. The other columns show the duration of the transmitting node’s infection at the time of transmission, the per-act transmission probability, act rate during the transmission, and final per-partnership transmission probability at that time step.
# A tibble: 10 × 8
# Groups: at, sus [10]
at sus inf network infDur transProb actRate finalProb
<int> <int> <int> <int> <dbl> <dbl> <dbl> <dbl>
1 2 26 385 1 6 0.4 2 0.64
2 2 61 4 1 3 0.4 2 0.64
3 2 72 80 1 6 0.4 2 0.64
4 2 169 67 1 9 0.4 2 0.64
5 2 342 67 1 9 0.4 2 0.64
6 2 416 385 1 6 0.4 2 0.64
7 3 41 227 1 2 0.4 2 0.64
8 3 72 80 1 7 0.4 2 0.64
9 3 169 67 1 10 0.4 2 0.64
10 3 284 98 1 2 0.4 2 0.6421.2.3.5 Data Exporting and Plotting with ggplot
We built in plotting methods directly for netsim class objects so you can easily plot multiple types of summary statistics from the simulated model object. However, if you prefer an external plotting tool in R, such as ggplot, it is easy to extract the data in tidy format for analysis and plotting. Here is an example of how to do so for our model above. See the help for the ggplot if you are unfamiliar with this syntax.
Code
df <- as.data.frame(sim)
df.mean <- as.data.frame(sim, 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()