37 Adding a Latent Stage to Make an SEIR Model
EpiModel includes an integrated SIR model, but here we show how to model an SEIR disease like COVID-19. The E compartment in this disease is an exposed state in which a person has been infected but is not infectious to others. Some infectious diseases have this latent non-infectious stage, and in general it provides a general framework for transmission risk that is dependent on one’s stage of disease progression.
This tutorial has two R scripts you should download: a primary script containing the code below and a separate module script. Download both and put them in the same working directory.
37.1 Setup
First start by loading EpiModel and clearing your global environment.
37.2 EpiModel Model Extensions
37.2.1 Conceptualization
The first step to any EpiModel extension model is to conceptually identify what new functionality, above and beyond the built-in models, is desired and where that functionality should be added. There are often many “right” answers to these questions, and this aspect is only learned over time. But in general, it is helpful to map out a model extension on a state/flow diagram to pinpoint where the additions should be concentrated.
37.2.1.1 SIR Module Set
In this particular model, we will be adding a new disease state that occurs in between the susceptible disease state and an infectious disease state in an SIR model. The standard, built-in SIR model in EpiModel uses set modules (elements of the model), each with its own set of associated functions (the realization of those elements in code). For any built-in model, you can see what set of modules and functions have been used by running the simulation with netsim
and then printing the output.
Code
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)
param <- param.net(inf.prob = 0.3, rec.rate = 0.1)
init <- init.net(i.num = 10, r.num = 0)
control <- control.net(type = "SIR", nsteps = 25, nsims = 1, verbose = FALSE)
mod1 <- netsim(est1, param, init, control)
mod1
EpiModel Simulation
=======================
Model class: netsim
Simulation Summary
-----------------------
Model type: SIR
No. simulations: 1
No. time steps: 25
No. NW groups: 1
Fixed Parameters
---------------------------
inf.prob = 0.3
rec.rate = 0.1
act.rate = 1
groups = 1
Model Output
-----------------------
Variables: s.num i.num r.num num si.flow ir.flow
Networks: sim1
Transmissions: sim1
Formation Statistics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 50 46.96 -6.08 1.2 -2.534 NA 2.638
Duration Statistics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 20 20.242 1.21 0.752 0.322 NA 1.447
Dissolution Statistics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 0.05 0.061 21.407 0.006 1.699 NA 0.031
37.2.1.2 Module Classes
This includes a series of modules, which we classify into Standard and Flexible modules as shown on the table below. Any modules may be modified, but standard modules are those that are typically not modified because they generalize the core internal processes for simulation object generation, network resimulation, and epidemic bookkeeping. Flexible modules, in contrast, are those that are likely necessary to modify for an EpiModel extension.
In the example above, we used all of these modules except the arrivals and departures modules because we had a closed population. In addition to these set of built-in modules, any user can add more modules to this set, depending on what is needed. Again, this involves some conceptualization of how to organize the model processes, including whether those processes are similar to the built-in modules or something new.
37.2.1.3 Modules vs. Functions
All modules have associated functions, and these are passed into the epidemic model run in netsim
through control.net
. Printing out the arguments for this function, you will see that each of the standard modules have default associated functions as inputs, and flexible modules have a default of NULL
.
function (type, nsteps, start = 1, nsims = 1, ncores = 1, resimulate.network = FALSE,
tergmLite = FALSE, cumulative.edgelist = FALSE, truncate.el.cuml = 0,
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, tergmLite.track.duration = FALSE, set.control.ergm = control.simulate.formula(MCMC.burnin = 2e+05),
set.control.tergm = control.simulate.formula.tergm(), save.diss.stats = TRUE,
dat.updates = NULL, ...)
NULL
For built-in models, EpiModel selects which modules are needed based on the model parameters and initial conditions. For each time step, the modules run in the order in which they are specified in the output here; this also matches the order in which they are listed in control.net
.
37.2.1.4 Built-in Functions as Templates
Any of the built-in functions associated with flexible modules are intended to be templates for user inspection and extension for research-level models. So, for example, mod1
above shows that infection.net
was used as the infection module. That function has a help page briefly describing what it does. And you can also inspect the function contents with:
We’ll use an edited down version of that function, with some additional explanation, below. In addition, the disease progression state transition within an SIR model (and an SIS model too) is handled by the recovery module, with an associated recovery.net
function.
We will use an edited down version of that function as a template for a new, more general disease progression module.
37.2.2 The EpiModel API
All EpiModel module functions have a set of shared design requirements. This set of requirements defines the EpiModel Application Programming Interface (API) for extension. The best way to learn this is through a concrete example like this, but here are the general API rules:
- Each module has an associated R function, the general design of which is:
The function takes an object called dat
, which is the master data object passed around by netsim
, performs some processes (e.g., infection, recovery, aging, interventions), updates the dat
object, and then returns that object. The other input argument to each function must be at
, which is a time step counter.
Data are stored on the
dat
object in a particular way: in sublists that are organized by category. The main categories of data to interact with include model inputs (parameters, initial conditions, and controls) from those three associated input functions; nodal attributes (e.g., an individual disease status for each person); and summary statistics (e.g., the disease prevalence at at time step). There are accessor functions for reading (these are theget_
functions) and writing (these are theset_
functions) to thedat
object in the appropriate place.The typical function design involves three steps: a) reading the relevant inputs from the
dat
object; b) performing some micro-level process on the nodes that is usually a function of fixed parameters and time-varying nodal attributes; c) writing the updated objects back on to thedat
object.
Let’s see how this API works by extending our infection and recovery functions to transition an SIR model into an SEIR model.
37.2.3 Infection Module
The built-in infection module for an SIR model performs the following functions listed in the function below. The core process is determining which edges are eligible for a disease transmission to occur, and then randomly simulating that transmission process. Why is it necessary to update the infection function for an SEIR model? Because an SIR involves a transition between S and I disease statuses when a new infection occurs, but an SEIR model involves a transition between S and E disease statuses. It’s a small, but important change. Here’s the full modified function, with embedded comments. Note that we can use the browser
function to run this function in debug mode by uncommenting the third line (we will demonstrate this).
Code
infect <- function(dat, at) {
## Uncomment this to run environment interactively
# browser()
## Attributes ##
active <- get_attr(dat, "active")
status <- get_attr(dat, "status")
infTime <- get_attr(dat, "infTime")
## Parameters ##
inf.prob <- get_param(dat, "inf.prob")
act.rate <- get_param(dat, "act.rate")
## Find infected nodes ##
idsInf <- which(active == 1 & status == "i")
nActive <- sum(active == 1)
nElig <- length(idsInf)
## Initialize default incidence at 0 ##
nInf <- 0
## If any infected nodes, proceed with transmission ##
if (nElig > 0 && nElig < nActive) {
## Look up discordant edgelist ##
del <- discord_edgelist(dat, at)
## If any discordant pairs, proceed ##
if (!(is.null(del))) {
# Set parameters on discordant edgelist data frame
del$transProb <- inf.prob
del$actRate <- act.rate
del$finalProb <- 1 - (1 - del$transProb)^del$actRate
# Stochastic transmission process
transmit <- rbinom(nrow(del), 1, del$finalProb)
# Keep rows where transmission occurred
del <- del[which(transmit == 1), ]
# Look up new ids if any transmissions occurred
idsNewInf <- unique(del$sus)
nInf <- length(idsNewInf)
# Set new attributes and transmission matrix
if (nInf > 0) {
status[idsNewInf] <- "e"
infTime[idsNewInf] <- at
dat <- set_attr(dat, "status", status)
dat <- set_attr(dat, "infTime", infTime)
dat <- set_transmat(dat, del, at)
}
}
}
## Save summary statistic for S->E flow
dat <- set_epi(dat, "se.flow", at, nInf)
return(dat)
}
Each step is relatively self-explanatory with the comments, but we will run this interactively at the end of this tutorial to step through the updated data structures. Infection is one of the more complex processes because it involves a dyadic process (a contact between an S and I node). That involves a construction of a discordant edgelist; that is, a list of edges in which there is a disease discordant dyad.
The main point here is: we have made a change to the infection module function, and it consists of updating the disease status of a newly infected persons to "e"
instead of "i"
. Additionally, we are tracking a new summary statistic, se.flow
that tracks the size of the flow from S to E based on the number of new infected, nInf
at the time step.
37.2.4 Progression Module
Next up is the disease progression module. Here, we have generalized the built-in recovery module function to handle two disease progression transitions after infection: E to I (latent to infectious stages) and I to R (infectious to recovered stages). Like many individual-level transitions, this involves flipping a weighted coin with rbinom
: this performs a series of random Bernoulli draws based on the specified parameters. Here is the full model function.
Code
progress <- function(dat, at) {
## Uncomment this to function environment interactively
# browser()
## Attributes ##
active <- get_attr(dat, "active")
status <- get_attr(dat, "status")
## Parameters ##
ei.rate <- get_param(dat, "ei.rate")
ir.rate <- get_param(dat, "ir.rate")
## E to I progression process ##
nInf <- 0
idsEligInf <- which(active == 1 & status == "e")
nEligInf <- length(idsEligInf)
if (nEligInf > 0) {
vecInf <- which(rbinom(nEligInf, 1, ei.rate) == 1)
if (length(vecInf) > 0) {
idsInf <- idsEligInf[vecInf]
nInf <- length(idsInf)
status[idsInf] <- "i"
}
}
## I to R progression process ##
nRec <- 0
idsEligRec <- which(active == 1 & status == "i")
nEligRec <- length(idsEligRec)
if (nEligRec > 0) {
vecRec <- which(rbinom(nEligRec, 1, ir.rate) == 1)
if (length(vecRec) > 0) {
idsRec <- idsEligRec[vecRec]
nRec <- length(idsRec)
status[idsRec] <- "r"
}
}
## Write out updated status attribute ##
dat <- set_attr(dat, "status", status)
## Save summary statistics ##
dat <- set_epi(dat, "ei.flow", at, nInf)
dat <- set_epi(dat, "ir.flow", at, nRec)
dat <- set_epi(dat, "e.num", at,
sum(active == 1 & status == "e"))
dat <- set_epi(dat, "r.num", at,
sum(active == 1 & status == "r"))
return(dat)
}
This set of two progression processes involves querying who is eligible to transition, randomly transitioning some of those eligible, updating the status
attribute for those who have progressed, and then recording some new summary statistics.
37.3 Network Model
With the epidemic modules defined, we will now step back to parameterize, estimate, and diagnose the TERGM. Here we will use a relatively basic model with an edges
and degree(1)
. Here we are not using any nodal attributes in either the TERGM or the epidemic modules, but these could be added (we will get more practice with that tomorrow). Note that we are using a relatively high mean degree (2 per capita) compared to some of our prior models, with lower than expected isolates.
Code
# Initialize the network
nw <- network_initialize(500)
# Define the formation model: edges + degree terms
formation = ~edges + degree(1)
# Input the appropriate target statistics for each term
target.stats <- c(500, 180)
# Parameterize the dissolution model
coef.diss <- dissolution_coefs(dissolution = ~offset(edges), duration = 25)
coef.diss
Dissolution Coefficients
=======================
Dissolution Model: ~offset(edges)
Target Statistics: 25
Crude Coefficient: 3.178054
Mortality/Exit Rate: 0
Adjusted Coefficient: 3.178054
Next we fit the network model.
And then diagnose it. We are including a wide range of degree terms to monitor so we can see the full degree distribution . Although the degree(1)
term does not look great visually, it is only off by ~5 edges in absolute terms.
Code
Network Diagnostics
-----------------------
- Simulating 10 networks
- Calculating formation statistics
EpiModel Network Diagnostics
=======================
Diagnostic Method: Dynamic
Simulations: 10
Time Steps per Sim: 500
Formation Diagnostics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 500 502.029 0.406 1.930 1.051 7.299 22.487
degree0 NA 50.752 NA 0.336 NA 1.353 7.234
degree1 180 174.653 -2.971 0.691 -7.735 2.527 11.717
degree2 NA 114.996 NA 0.400 NA 2.326 9.579
degree3 NA 83.793 NA 0.349 NA 1.529 8.520
degree4 NA 45.501 NA 0.300 NA 1.806 7.106
degree5 NA 19.998 NA 0.203 NA 1.035 4.666
degree6 NA 7.218 NA 0.116 NA 0.541 2.789
degree7 NA 2.229 NA 0.056 NA 0.197 1.492
Duration Diagnostics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 25 25.093 0.371 0.086 1.078 0.232 1.032
Dissolution Diagnostics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 0.04 0.04 -0.218 0 -0.682 0 0.009
37.4 Epidemic Model Parameterization
The epidemic model parameterization for any extension model consists of using the same three input functions as with built-in models.
First we start with the model parameters in param.net
. Note that we have two new rates here: ei.rate
and ir.rate
. These control the transitions from E to I, and then I to R. We have selected these parameter names and input the values here, but note that the same parameters get pulled into the disease progression function above. So in general there must be consistency between the naming of the parameters as inputs and their references in the model functions.
Next we specify the initial conditions. Here we are specifying that there are 10 infectious individuals at the epidemic outset. If we wanted to initialize the model with persons only in the latent stage, we would need to set disease status as a nodal attribute on the network instead (following the same approach in the Serosorting Tutorial in Module 8.
Finally, we specify the control settings in control.net
. Extension models like this require some significant updates (compared to built-in models) here. First, type
is set to NULL
in any extension model because we are no longer using EpiModel to pre-select which modules to run; it is an entirely manual process. Second, the new input parameter for the infection.FUN
argument is infect
; in other words, our new function for our infection module is the one we built above. Note that it is called infect
. Third, we are defining a new module, progress.FUN
, with an argument that matches the name of our new progression function. We defined a new module, rather than replaced the function of the recovery module, just to show you how to define a new module here. We could have just as easily done it the latter way with the same results. Finally, we are explicitly setting resimulate.network
to FALSE
; this is the default, so this is not necessary, but it is just to remind you this is a model without network feedback.
In the R scripts that we include with this tutorial, you will see that we have two separate R script files (for the R Markdown build, they all go in this single file). One contains the module functions, and the other contains all the other code to parameterize and run the model. We do this because it allows for easier interaction with the functions in browser mode. I will demonstrate this live next. But in general, placing the functions in a separate file conceptually disentangles the model functionality from the model parameterization. It is critical, however, that you source in the file containing the functions before you run control.net
(otherwise, control.net
does not know what infect
and progress
are).
37.5 Epidemic Simulation
Finally we are ready to do the epidemic model simulation and analysis. This is done using the same approach as the built-in models. We will start with running one model in the browser mode interactively.
And then go back, comment out the browser
lines, re-source the functions, and run a full-scale model with 10 simulations.
Once the model simulation is complete, we can work with the model object just like a built-in model. Start by printing the output to see what is available.
EpiModel Simulation
=======================
Model class: netsim
Simulation Summary
-----------------------
Model type:
No. simulations: 10
No. time steps: 500
No. NW groups: 1
Fixed Parameters
---------------------------
inf.prob = 0.5
act.rate = 2
ei.rate = 0.01
ir.rate = 0.01
groups = 1
Model Functions
-----------------------
initialize.FUN
resim_nets.FUN
summary_nets.FUN
infection.FUN
nwupdate.FUN
prevalence.FUN
verbose.FUN
progress.FUN
Model Output
-----------------------
Variables: s.num i.num num ei.flow ir.flow e.num r.num
se.flow
Networks: sim1 ... sim10
Transmissions: sim1 ... sim10
Formation Statistics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 500 502.885 0.577 1.859 1.552 5.675 21.751
degree1 180 173.479 -3.623 0.752 -8.666 2.916 11.865
Duration Statistics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 25 25.018 0.073 0.081 0.227 0.327 1.017
Dissolution Statistics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 0.04 0.04 0.152 0 0.493 0 0.009
Here is the default plot with all the compartment sizes over time. This includes the new summary statistics we tracked in the disease progression function.
Here are the flow sizes, including the new se.flow
incidence tracker that we established in the new infection function.
Finally, here is the data frame output from the model, with rows limited to time step 100 across all 10 simulations.
sim time s.num i.num num ei.flow ir.flow e.num r.num se.flow
100 1 100 200 73 500 2 0 200 25 2
600 2 100 219 64 500 1 1 180 34 3
1100 3 100 250 51 500 4 2 168 26 5
1600 4 100 260 56 500 2 1 154 26 4
2100 5 100 318 40 500 0 0 124 16 2
2600 6 100 177 83 500 1 2 201 38 1
3100 7 100 184 72 500 2 1 210 32 2
3600 8 100 205 76 500 1 1 169 47 3
4100 9 100 280 45 500 0 1 157 16 2
4600 10 100 224 65 500 2 0 179 29 3
And here are the transmission matrices and their phylogram plots (one per infection seed).