library(EpiModel)
rm(list = ls())41 Adding An Asymptomatic Pathway and Screening-Based Interventions to the COVID Model
In this tutorial, we will build on our COVID-19 model by adding a more complex disease progression structure, with clinical versus sub-clinical (asymptomatic) pathways, as well as disease screening and diagnosis-based interventions.
41.1 Setup
First start by loading EpiModel and clearing your global environment.
We will use the network model parameterization and extension functions from the previous tutorial as the starting point for this tutorial.
41.2 Network Initialization
The network will have the same age and aging structure as the previous model. We have to reinitialize the mortality rates as in the last model.
ages <- 0:85
departure_rate <- c(588.45, 24.8, 11.7, 14.55, 47.85, 88.2, 105.65, 127.2,
154.3, 206.5, 309.3, 495.1, 736.85, 1051.15, 1483.45,
2294.15, 3642.95, 6139.4, 13938.3)
dr_pp_pd <- departure_rate / 1e5 / 365
age_spans <- c(1, 4, rep(5, 16), 1)
dr_vec <- rep(dr_pp_pd, times = age_spans)The network size and age attribute are next added.
n <- 1000
nw <- network_initialize(n)
ageVec <- sample(ages, n, replace = TRUE)
nw <- set_vertex_attribute(nw, "age", ageVec)We will also initialize the status vector as in the past model, plus add three new nodal attributes that are needed for the new extension epidemic models. statusTime will track when a node changes disease stage; it is initialized as NA for everyone and then set to time step 1 for those who are infected. clinical will be a binary attribute to record whether nodes are in a clinical or subclinical pathway; this will be initialized as NA for everyone because it is updated at the move of transition out of the latent (E) state. Finally, dxStatus will keep track of the diagnosed nodal status; it will have three values (0 for never screened, 1 for screened negative, and 2 for screened positive), but everyone will be initialized as 0. As with the others, all nodal attributes are set on the network object with set_vertex_attribute.
statusVec <- rep("s", n)
init.latent <- sample(1:n, 50)
statusVec[init.latent] <- "e"
statusTime <- rep(NA, n)
statusTime[which(statusVec == "e")] <- 1
clinical <- rep(NA, n)
dxStatus <- rep(0, n)
nw <- set_vertex_attribute(nw, "status", statusVec)
nw <- set_vertex_attribute(nw, "statusTime", statusTime)
nw <- set_vertex_attribute(nw, "clinical", clinical)
nw <- set_vertex_attribute(nw, "dxStatus", dxStatus)
nw Network attributes:
vertices = 1000
directed = FALSE
hyper = FALSE
loops = FALSE
multiple = FALSE
bipartite = FALSE
total edges= 0
missing edges= 0
non-missing edges= 0
Vertex attribute names:
age clinical dxStatus status statusTime vertex.names
No edge attributes41.3 New Modules
The new epidemic modules will involve some minor and major updates to our extension modules from the previous tutorial, and the addition of a new module to handle COVID screening and diagnosis.
41.3.1 Progression
The more sophisticated progression module that we develop below corresponds to this flow diagram, which is an extension of the SEIR diagram. Here, persons who are infected enter one of two pathways: a clinical (symptomatic) pathway and a sub-clinical (asymptomatic) pathway. This follows emerging research suggesting that a substantial fraction of COVID cases never experience any symptoms yet still transmit (although likely at a lower rate). In the clinical pathway, persons are further subdivided into a pre-symptomatic and symptomatic phase to reflect that infectiousness may occur just prior to symptoms.
In the updated progression function below, we build out in code what is shown in the figure above. Each transition from one state to another involves a random Bernoulli process with the probability equal to the transition rate (which is the reciprocal of the average time spent in the current state before transition).
To determine how people enter the subclinical versus clinical pathway, we will use the prop.clinical parameter, which is actually a vector of probabilities corresponding to decade of age. This will allow the level of asymptomatic infection to decrease over age, as evidence suggests. The age.group calculation involves rounding down the continous age in years to a decade, which is then used to index the prop.clinical vector. The clinical attribute is then assigned and stays with each person throughout their infection.
Note that we also include a feature here of tracking the individual statusTime, which is the time step of transition from one state to another. This allows us to prevent immediate transitions across multiple states within a single time step, since the condition for transition requires statusTime < at.
progress2 <- function(dat, at) {
## Attributes
active <- get_attr(dat, "active")
status <- get_attr(dat, "status")
age <- get_attr(dat, "age")
statusTime <- get_attr(dat, "statusTime")
clinical <- get_attr(dat, "clinical")
## Parameters
prop.clinical <- get_param(dat, "prop.clinical")
ea.rate <- get_param(dat, "ea.rate")
ar.rate <- get_param(dat, "ar.rate")
eip.rate <- get_param(dat, "eip.rate")
ipic.rate <- get_param(dat, "ipic.rate")
icr.rate <- get_param(dat, "icr.rate")
## Determine Subclinical (E to A) or Clinical (E to Ip to Ic) pathway
ids.newInf <- which(active == 1 & status == "e" & statusTime <= at & is.na(clinical))
num.newInf <- length(ids.newInf)
if (num.newInf > 0) {
age.group <- pmin((floor(age[ids.newInf] / 10)) + 1, 8)
prop.clin.vec <- prop.clinical[age.group]
if (any(is.na(prop.clin.vec))) stop("error in prop.clin.vec")
vec.new.clinical <- rbinom(num.newInf, 1, prop.clin.vec)
clinical[ids.newInf] <- vec.new.clinical
}
## Subclinical Pathway
# E to A: latent move to asymptomatic infectious
num.new.EtoA <- 0
ids.Es <- which(active == 1 & status == "e" & statusTime < at & clinical == 0)
num.Es <- length(ids.Es)
if (num.Es > 0) {
vec.new.A <- which(rbinom(num.Es, 1, ea.rate) == 1)
if (length(vec.new.A) > 0) {
ids.new.A <- ids.Es[vec.new.A]
num.new.EtoA <- length(ids.new.A)
status[ids.new.A] <- "a"
statusTime[ids.new.A] <- at
}
}
# A to R: asymptomatic infectious move to recovered
num.new.AtoR <- 0
ids.A <- which(active == 1 & status == "a" & statusTime < at & clinical == 0)
num.A <- length(ids.A)
if (num.A > 0) {
vec.new.R <- which(rbinom(num.A, 1, ar.rate) == 1)
if (length(vec.new.R) > 0) {
ids.new.R <- ids.A[vec.new.R]
num.new.AtoR <- length(ids.new.R)
status[ids.new.R] <- "r"
statusTime[ids.new.R] <- at
}
}
## Clinical Pathway
# E to Ip: latent move to preclinical infectious
num.new.EtoIp <- 0
ids.Ec <- which(active == 1 & status == "e" & statusTime < at & clinical == 1)
num.Ec <- length(ids.Ec)
if (num.Ec > 0) {
vec.new.Ip <- which(rbinom(num.Ec, 1, eip.rate) == 1)
if (length(vec.new.Ip) > 0) {
ids.new.Ip <- ids.Ec[vec.new.Ip]
num.new.EtoIp <- length(ids.new.Ip)
status[ids.new.Ip] <- "ip"
statusTime[ids.new.Ip] <- at
}
}
# Ip to Ic: preclinical infectious move to clinical infectious
num.new.IptoIc <- 0
ids.Ip <- which(active == 1 & status == "ip" & statusTime < at & clinical == 1)
num.Ip <- length(ids.Ip)
if (num.Ip > 0) {
vec.new.Ic <- which(rbinom(num.Ip, 1, ipic.rate) == 1)
if (length(vec.new.Ic) > 0) {
ids.new.Ic <- ids.Ip[vec.new.Ic]
num.new.IptoIc <- length(ids.new.Ic)
status[ids.new.Ic] <- "ic"
statusTime[ids.new.Ic] <- at
}
}
# Ic to R: clinical infectious move to recovered (if not mortality first)
num.new.IctoR <- 0
ids.Ic <- which(active == 1 & status == "ic" & statusTime < at & clinical == 1)
num.Ic <- length(ids.Ic)
if (num.Ic > 0) {
vec.new.R <- which(rbinom(num.Ic, 1, icr.rate) == 1)
if (length(vec.new.R) > 0) {
ids.new.R <- ids.Ic[vec.new.R]
num.new.IctoR <- length(ids.new.R)
status[ids.new.R] <- "r"
statusTime[ids.new.R] <- at
}
}
## Save updated status attribute
dat <- set_attr(dat, "status", status)
dat <- set_attr(dat, "statusTime", statusTime)
dat <- set_attr(dat, "clinical", clinical)
## Save summary statistics
dat <- set_epi(dat, "ea.flow", at, num.new.EtoA)
dat <- set_epi(dat, "ar.flow", at, num.new.AtoR)
dat <- set_epi(dat, "eip.flow", at, num.new.EtoIp)
dat <- set_epi(dat, "ipic.flow", at, num.new.IptoIc)
dat <- set_epi(dat, "icr.flow", at, num.new.IctoR)
dat <- set_epi(dat, "e.num", at, sum(status == "e"))
dat <- set_epi(dat, "a.num", at, sum(status == "a"))
dat <- set_epi(dat, "ip.num", at, sum(status == "ip"))
dat <- set_epi(dat, "ic.num", at, sum(status == "ic"))
dat <- set_epi(dat, "r.num", at, sum(status == "r"))
return(dat)
}At the end of the function, we reset the relevant attributes that have changed on the dat object and keep track of all the flow sizes and state sizes. There are several more flows and states to track now, compared to the earlier SEIR model!
41.3.2 Diagnosis
The diagnosis module will handle the process for screening of cases, which here is controlled by two screening rate parameters, dx.rate.sympt and dx.rate.other. The former parameter controls the rate of screening for persons currently with symptomatic infection (that is, in the ic disease state), while the latter parameter controls the rate for all other persons. This reflects the higher rates of symptoms-based diagnosis of active cases.
Additionally, we have a logical parameter, allow.rescreen, that controls whether persons who have previously had a negative COVID test can subsequently retest (this is why we wanted to track dxStatus as a three-level variables of never-tested, tested-negative, and tested-positive). Finally, because COVID diagnostics are imperfect, we incorporate PCR sensitive parameter, pcr.sens, to simulate the process of false-negative test results.
dx_covid <- function(dat, at) {
## Pull attributes
active <- get_attr(dat, "active")
status <- get_attr(dat, "status")
dxStatus <- get_attr(dat, "dxStatus")
## Pull parameters
dx.rate.sympt <- get_param(dat, "dx.rate.sympt")
dx.rate.other <- get_param(dat, "dx.rate.other")
allow.rescreen <- get_param(dat, "allow.rescreen")
pcr.sens <- get_param(dat, "pcr.sens")
## Initialize trackers
idsDx.sympt <- idsDx.other <- NULL
idsDx.sympt.pos <- idsDx.other.pos.true <- NULL
idsDx.sympt.neg <- idsDx.other.pos.false <- NULL
## Determine screening eligibility
idsElig.sympt <- which(active == 1 & dxStatus %in% 0:1 & status == "ic")
if (allow.rescreen == TRUE) {
idsElig.other <- which(active == 1 & dxStatus %in% 0:1 &
status %in% c("s", "e", "a", "ip", "r"))
} else {
idsElig.other <- which(active == 1 & dxStatus == 0 &
status %in% c("s", "e", "a", "ip", "r"))
}
## Symptomatic testing
nElig.sympt <- length(idsElig.sympt)
if (nElig.sympt > 0) {
vecDx.sympt <- which(rbinom(nElig.sympt, 1, dx.rate.sympt) == 1)
idsDx.sympt <- idsElig.sympt[vecDx.sympt]
nDx.sympt <- length(idsDx.sympt)
if (nDx.sympt > 0) {
vecDx.sympt.pos <- rbinom(nDx.sympt, 1, pcr.sens)
idsDx.sympt.pos <- idsDx.sympt[which(vecDx.sympt.pos == 1)]
idsDx.sympt.neg <- idsDx.sympt[which(vecDx.sympt.pos == 0)]
dxStatus[idsDx.sympt.pos] <- 2
dxStatus[idsDx.sympt.neg] <- 1
}
}
## Asymptomatic screening
nElig.other <- length(idsElig.other)
if (nElig.other > 0) {
vecDx.other <- which(rbinom(nElig.other, 1, dx.rate.other) == 1)
idsDx.other <- idsElig.other[vecDx.other]
nDx.other <- length(idsDx.other)
if (nDx.other > 0) {
idsDx.other.neg <- intersect(idsDx.other, which(status == "s"))
idsDx.other.pos.all <- intersect(idsDx.other,
which(status %in% c("e", "a", "ip", "r")))
vecDx.other.pos <- rbinom(length(idsDx.other.pos.all), 1, pcr.sens)
idsDx.other.pos.true <- idsDx.other.pos.all[which(vecDx.other.pos == 1)]
idsDx.other.pos.false <- idsDx.other.pos.all[which(vecDx.other.pos == 0)]
dxStatus[idsDx.other.neg] <- 1
dxStatus[idsDx.other.pos.false] <- 1
dxStatus[idsDx.other.pos.true] <- 2
}
}
## Set attr
dat <- set_attr(dat, "dxStatus", dxStatus)
## Summary statistics
dat <- set_epi(dat, "nDx", at, length(idsDx.sympt) + length(idsDx.other))
dat <- set_epi(dat, "nDx.pos", at, length(idsDx.sympt.pos) +
length(idsDx.other.pos.true))
dat <- set_epi(dat, "nDx.pos.sympt", at, length(idsDx.sympt.pos))
dat <- set_epi(dat, "nDx.pos.fn", at, length(idsDx.sympt.neg) +
length(idsDx.other.pos.false))
return(dat)
}At the end of the function, we updated the modified dxStatus attribute, and calculate some summary statistics for new cases.
41.3.3 Infection
It is also necessary to update the infection module function in a couple of ways. The first will reflect that we now have 3 infectious disease states, of varying infectiousness, compared to the earlier SEIR model’s one infectious state. This requires modifying the code querying the definition of infectious nodes, and the construction of the discordant edgelist.
Second, onto the discordant edgelist data frame, we add the disease state of the infectious node and the diagnostic status of that node. Those new data are then used in two ways. First, being in the asymptomatic disease state, a, is associated with a lower probability of transmission compared to the other two (clinical) infectious disease states. That is accomplished by modifying the base transmission probability by a relative risk parameter, inf.prob.a.rr.
Second, being infectious and diagnosed positive (which a dxState of 2), may trigger behavioral interventions reflect case isolation. This type of behavioral change may be accomplished in several different ways. Here it involves a modification of the act.rate parameter that controls the number of individual exposure events between active dyads in the current time step. This type of intervention may start at a particular time step, act.rate.dx.inter.time, and result in a relative reduction in the current act rate of act.rate.dx.inter.rr.
infect2 <- 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")
dxStatus <- get_attr(dat, "dxStatus")
statusTime <- get_attr(dat, "statusTime")
## Parameters ##
inf.prob <- get_param(dat, "inf.prob")
act.rate <- get_param(dat, "act.rate")
inf.prob.a.rr <- get_param(dat, "inf.prob.a.rr")
act.rate.dx.inter.time <- get_param(dat, "act.rate.dx.inter.time")
act.rate.dx.inter.rr <- get_param(dat, "act.rate.dx.inter.rr")
## Find infected nodes ##
infstat <- c("a", "ic", "ip")
idsInf <- which(active == 1 & status %in% infstat)
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, infstat = infstat)
## If any discordant pairs, proceed ##
if (!(is.null(del))) {
del$status <- status[del$inf]
del$dxStatus <- dxStatus[del$inf]
# Set parameters on discordant edgelist data frame
del$transProb <- inf.prob
del$transProb[del$status == "a"] <- del$transProb[del$status == "a"] *
inf.prob.a.rr
del$actRate <- act.rate
if (at >= act.rate.dx.inter.time) {
del$actRate[del$dxStatus == 2] <- del$actRate[del$dxStatus == 2] *
act.rate.dx.inter.rr
}
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
statusTime[idsNewInf] <- at
dat <- set_transmat(dat, del, at)
}
}
}
dat <- set_attr(dat, "status", status)
dat <- set_attr(dat, "infTime", infTime)
dat <- set_attr(dat, "statusTime", statusTime)
## Save summary statistics
dat <- set_epi(dat, "se.flow", at, nInf)
return(dat)
}There are no other modifications of the infection module other than to track the statusTime upon infection and then resetting that attribute on the dat object.
41.3.4 Births
The birth module requires a very minor change to update the three new nodal attributes on the dat object for any incoming nodes.
afunc2 <- function(dat, at) {
## Parameters ##
n <- get_epi(dat, "num", at - 1)
a.rate <- get_param(dat, "arrival.rate")
## Process ##
nArrivalsExp <- n * a.rate
nArrivals <- rpois(1, nArrivalsExp)
# Update attributes
if (nArrivals > 0) {
dat <- append_core_attr(dat, at = at, n.new = nArrivals)
dat <- append_attr(dat, "status", "s", nArrivals)
dat <- append_attr(dat, "infTime", NA, nArrivals)
dat <- append_attr(dat, "age", 0, nArrivals)
dat <- append_attr(dat, "statusTime", NA, nArrivals)
dat <- append_attr(dat, "clinical", NA, nArrivals)
dat <- append_attr(dat, "dxStatus", 0, nArrivals)
}
## Summary statistics ##
dat <- set_epi(dat, "a.flow", at, nArrivals)
return(dat)
}41.3.5 Deaths
The death module requires an even more minor modification from the last tutorial, which involves simulating COVID-related mortality and tracking the number of covid.deaths based on deaths that occurs within the ic (clinical symptomatic) state (previously these were in the i state only).
dfunc2 <- function(dat, at) {
## Attributes
active <- get_attr(dat, "active")
exitTime <- get_attr(dat, "exitTime")
age <- get_attr(dat, "age")
status <- get_attr(dat, "status")
## Parameters
dep.rates <- get_param(dat, "departure.rates")
dep.dis.mult <- get_param(dat, "departure.disease.mult")
## Query alive
idsElig <- which(active == 1)
nElig <- length(idsElig)
## Initialize trackers
nDepts <- 0
idsDepts <- NULL
if (nElig > 0) {
## Calculate age-specific departure rates for each eligible node ##
## Everyone older than 85 gets the final mortality rate
whole_ages_of_elig <- pmin(ceiling(age[idsElig]), 86)
departure_rates_of_elig <- dep.rates[whole_ages_of_elig]
## Multiply departure rates for diseased persons
idsElig.inf <- which(status[idsElig] == "ic")
departure_rates_of_elig[idsElig.inf] <- departure_rates_of_elig[idsElig.inf] * dep.dis.mult
## Simulate departure process
vecDepts <- which(rbinom(nElig, 1, departure_rates_of_elig) == 1)
idsDepts <- idsElig[vecDepts]
nDepts <- length(idsDepts)
## Update nodal attributes
if (nDepts > 0) {
active[idsDepts] <- 0
exitTime[idsDepts] <- at
}
}
## Reset attributes
dat <- set_attr(dat, "active", active)
dat <- set_attr(dat, "exitTime", exitTime)
## Summary statistics ##
dat <- set_epi(dat, "total.deaths", at, nDepts)
# covid deaths
covid.deaths <- length(intersect(idsDepts, which(status == "ic")))
dat <- set_epi(dat, "covid.deaths", at, covid.deaths)
return(dat)
}41.4 Network Model Estimation
With epidemic modules designed, we parameterize the exact same network model as the previous tutorial.
formation <- ~edges + degree(0) + absdiff("age")
mean_degree <- 2
edges <- mean_degree * (n/2)
avg.abs.age.diff <- 2
isolates <- n * 0.08
absdiff <- edges * avg.abs.age.diff
target.stats <- c(edges, isolates, absdiff)
target.stats[1] 1000 80 2000coef.diss <- dissolution_coefs(~offset(edges), 20, mean(dr_vec))
coef.dissDissolution Coefficients
=======================
Dissolution Model: ~offset(edges)
Target Statistics: 20
Crude Coefficient: 2.944439
Mortality/Exit Rate: 3.159205e-05
Adjusted Coefficient: 2.945703We estimate the model with netest. Here we demonstrate how to increase the maximum number of MCMLE iterations (the default is 20), which was sometimes necessary to get this model to converge.
est <- netest(nw, formation, target.stats, coef.diss,
set.control.ergm = control.ergm(MCMLE.maxit = 100))Model diagnostics look similar to last time.
dx <- netdx(est, nsims = 10, ncores = 5, nsteps = 500,
nwstats.formula = ~edges + absdiff("age") + degree(0:6),
set.control.tergm =
control.simulate.formula.tergm(MCMC.burnin.min = 20000))
Network Diagnostics
-----------------------
- Simulating 10 networks
- Calculating formation statisticsprint(dx)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 1000 982.169 -1.783 1.771 -10.068 5.586 25.833
absdiff.age 2000 1963.437 -1.828 6.519 -5.609 13.908 84.133
degree0 80 90.395 12.993 0.443 23.478 1.275 9.996
degree1 NA 323.522 NA 0.692 NA 1.851 15.502
degree2 NA 291.009 NA 0.534 NA 2.168 14.251
degree3 NA 175.786 NA 0.468 NA 1.103 12.182
degree4 NA 79.719 NA 0.379 NA 1.365 9.091
degree5 NA 28.365 NA 0.231 NA 1.090 5.625
degree6 NA 8.416 NA 0.116 NA 0.707 3.043
Duration Diagnostics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 20 19.938 -0.308 0.052 -1.182 0.125 0.633
Dissolution Diagnostics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 0.05 0.05 0.183 0 0.939 0 0.007plot(dx)
41.5 Epidemic Model Simulation
We are now using more realistic COVID parameters for disease progression, based on the current literature and other models. Note that the current parameters allow for 50% lower transmissibility in the asymptomatic stage, age-varying clinical pathways, a PCR sensitivity of 80%, a diagnosis rate that is 20-fold higher for those with symptomatic infection, and no rescreening. We have also prevented any case isolation by setting the act.rate.dx.inter.time to Inf.
param <- param.net(inf.prob = 0.1,
act.rate = 3,
departure.rates = dr_vec,
departure.disease.mult = 1000,
arrival.rate = 1/(365*85),
inf.prob.a.rr = 0.5,
act.rate.dx.inter.time = Inf,
act.rate.dx.inter.rr = 0.05,
# proportion in clinical pathway by age decade
prop.clinical = c(0.40, 0.25, 0.37, 0.42, 0.51, 0.59, 0.72, 0.76),
ea.rate = 1/4.0,
ar.rate = 1/5.0,
eip.rate = 1/4.0,
ipic.rate = 1/1.5,
icr.rate = 1/3.5,
pcr.sens = 0.8,
dx.rate.sympt = 0.2,
dx.rate.other = 0.01,
allow.rescreen = FALSE)
init <- init.net()For the control settings, it is necessary to define all the relevant modules for our system, and input the associated functions. We will simulate the model only over 100 days, with 10 simulations. Here we use tergmLite, but this can be set to FALSE to retain the full network data (these simulations will take a bit longer).
source("d5-s6-COVIDScreen-fx.R")
control <- control.net(type = NULL,
nsims = 10,
ncores = 5,
nsteps = 100,
infection.FUN = infect2,
progress.FUN = progress2,
dx.FUN = dx_covid,
aging.FUN = aging,
departures.FUN = dfunc2,
arrivals.FUN = afunc2,
resimulate.network = TRUE,
tergmLite = TRUE,
set.control.tergm =
control.simulate.formula.tergm(MCMC.burnin.min = 10000))The model is then simulated with netsim.
sim <- netsim(est, param, init, control)41.6 Model Analysis
Let’s print out the netsim object to review the available data variables.
simEpiModel Simulation
=======================
Model class: netsim
Simulation Summary
-----------------------
Model type:
No. simulations: 10
No. time steps: 100
No. NW groups: 1
Fixed Parameters
---------------------------
inf.prob = 0.1
act.rate = 3
departure.rates = 1.612192e-05 6.794521e-07 6.794521e-07 6.794521e-07
6.794521e-07 3.205479e-07 3.205479e-07 3.205479e-07 3.205479e-07 3.205479e-07
...
departure.disease.mult = 1000
arrival.rate = 3.223207e-05
inf.prob.a.rr = 0.5
act.rate.dx.inter.time = Inf
act.rate.dx.inter.rr = 0.05
prop.clinical = 0.4 0.25 0.37 0.42 0.51 0.59 0.72 0.76
ea.rate = 0.25
ar.rate = 0.2
eip.rate = 0.25
ipic.rate = 0.6666667
icr.rate = 0.2857143
pcr.sens = 0.8
dx.rate.sympt = 0.2
dx.rate.other = 0.01
allow.rescreen = FALSE
groups = 1
Model Functions
-----------------------
initialize.FUN
resim_nets.FUN
infection.FUN
departures.FUN
arrivals.FUN
nwupdate.FUN
prevalence.FUN
verbose.FUN
progress.FUN
dx.FUN
aging.FUN
Model Output
-----------------------
Variables: s.num i.num num ea.flow ar.flow eip.flow
ipic.flow icr.flow e.num a.num ip.num ic.num r.num nDx
nDx.pos nDx.pos.sympt nDx.pos.fn meanAge se.flow
total.deaths covid.deaths a.flow
Transmissions: sim1 ... sim10
Formation Statistics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 1000 949.798 -5.020 4.590 -10.938 12.438 30.568
degree0 80 91.969 14.961 0.954 12.548 3.206 9.755
absdiff.age 2000 1907.666 -4.617 14.302 -6.456 44.835 96.952
Duration and Dissolution Statistics
-----------------------
Not available when:
- `control$tergmLite == TRUE`
- `control$save.network == FALSE`
- `control$save.diss.stats == FALSE`
- dissolution formula is not `~ offset(edges)`
- `keep.diss.stats == FALSE` (if merging)Here is a plot of disease state sizes over time. The default plot makes it difficult to see the prevalence because the susceptible and recovered state sizes are much larger overall, so we plot just the three infectious states alone.
par(mfrow = c(1,2))
plot(sim)
plot(sim, y = c("a.num", "ip.num", "ic.num"), legend = TRUE)
Here are the first three transitions after leaving the suceptible state.
par(mfrow = c(1, 1))
plot(sim, y = c("se.flow", "ea.flow", "eip.flow"), legend = TRUE)
Finally, we can export the mean data, averaged across the 10 simulations.
df <- as.data.frame(sim, out = "mean")
head(df) time s.num i.num num ea.flow ar.flow eip.flow ipic.flow icr.flow e.num
1 1 950.0 0 1000.0 NaN NaN NaN NaN NaN NaN
2 2 945.1 0 1000.0 5.4 0.0 6.2 0.0 0.0 38.4
3 3 937.6 0 999.8 4.6 0.8 5.0 4.3 0.0 33.7
4 4 929.1 0 999.5 4.1 1.5 5.1 3.3 1.4 31.9
5 5 920.7 0 999.5 6.0 2.8 5.3 5.3 1.8 29.1
6 6 913.4 0 998.9 3.1 3.1 5.0 6.1 2.5 29.6
a.num ip.num ic.num r.num nDx nDx.pos nDx.pos.sympt nDx.pos.fn meanAge
1 NaN NaN NaN NaN NaN NaN NaN NaN NaN
2 5.4 6.2 0.0 0.0 9.1 0.5 0.0 0.0 43.25074
3 9.2 6.9 4.3 0.8 9.4 1.8 1.3 0.2 43.24598
4 11.8 8.7 6.1 3.7 9.8 1.5 1.0 0.4 43.24125
5 15.0 8.7 9.3 8.3 11.5 1.3 0.8 0.3 43.23986
6 15.0 7.6 12.7 13.9 11.3 1.9 1.2 0.6 43.22710
se.flow total.deaths covid.deaths a.flow
1 NaN NaN NaN NaN
2 4.9 0.1 0.0 0.1
3 7.4 0.2 0.1 0.0
4 8.5 0.3 0.3 0.0
5 8.6 0.2 0.2 0.2
6 7.2 0.6 0.4 0.0And calculate the average cumulative incidence.
sum(df$se.flow, na.rm = TRUE)[1] 334.9