Excluding forbidden ties in starting networks in EpiModel
Steve Goodreau
1. Intro
The goal is simple; we need to:
- create a network of MSM in which men have the standard roles (I = exclusively insertive, R = exclusively receptive, V = versatile)
- estimate a simple stergm for that network with a given mean degree and where II and RR ties are forbidden
- simulate dynamic networks (eventually with transmission dynamics) on this model fit, ensuring that II and RR ties are absent from the very beginning of that simulation.
This seems like it should be a fairly basic (and important) task that can be done with relatively simple code. As we will see, it is in ergm and stergm (sans transmission, of course), but does not seem to be in EpiModel. Hopefully that is just because I am overlooking some option.
To demonstrate where the functionality seems to break down, we will first do this using just ergm, then again using stergm, then again using EpiModel. We load all needed packages:
library(EpiModel)We create an empty network and assign the role attribute:
msm_nw <- network.initialize(10000, directed = F)
role <- sample(c('I', 'R', 'V'), size = 10000, replace = T, prob = c(0.24, 0.27, 0.49))
msm_nw <- set.vertex.attribute(msm_nw, "role", role)And we specify our edges target statistic:
nEdges <- 0.7 * 50002. Starting without forbidden ties in ERGM
First, we demonstrate that our goal of generating a network without forbidden ties is easily doable in ergm. We fit the model, with -Inf offsets for the two types of forbidden ties (II and RR):
msm_mod_ergm <- ergm(formula = msm_nw ~ edges + offset(nodematch('role', diff = TRUE, keep = 1:2)),
target = nEdges,
offset.coef = c(-Inf, -Inf))
summary(msm_mod_ergm)
==========================
Summary of model fit
==========================
Formula: nw ~ edges + offset(nodematch("role", diff = TRUE, keep = 1:2))
<environment: 0x00000000210ab470>
Iterations: 5 out of 20
Monte Carlo MLE Results:
Estimate Std. Error MCMC % p-value
edges -9.43061 0.01653 0 <1e-04 ***
nodematch.role.I -Inf 0.00000 0 <1e-04 ***
nodematch.role.R -Inf 0.00000 0 <1e-04 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Null Deviance: 69307787 on 49995000 degrees of freedom
Residual Deviance: NaN on 49994997 degrees of freedom
AIC: NaN BIC: NaN (Smaller is better.)
The following terms are fixed by offset and are not estimated:
nodematch.role.I nodematch.role.R If we try to simulate using ergm with all default control parameters, we hit a snag; there are II and RR in the first dozen or more networks:
msm_sim_ergm_no <- simulate(msm_mod_ergm, nsim = 50)
summary(msm_sim_ergm_no~edges+nodemix("role"))[1:20,c(2,4)] mix.role.I.I mix.role.R.R
[1,] 18 30
[2,] 16 22
[3,] 11 21
[4,] 9 18
[5,] 9 16
[6,] 9 16
[7,] 9 13
[8,] 7 12
[9,] 4 9
[10,] 3 7
[11,] 2 4
[12,] 2 4
[13,] 1 3
[14,] 1 2
[15,] 1 2
[16,] 1 1
[17,] 0 1
[18,] 0 1
[19,] 0 1
[20,] 0 1But this is easily rectified by setting a longer simulation burnin:
msm_sim_ergm_yes <- simulate(msm_mod_ergm, nsim = 50,
control=control.simulate.ergm(MCMC.burnin = 1e5))
summary(msm_sim_ergm_yes~edges+nodemix("role"))[1:20,c(2,4)] mix.role.I.I mix.role.R.R
[1,] 0 0
[2,] 0 0
[3,] 0 0
[4,] 0 0
[5,] 0 0
[6,] 0 0
[7,] 0 0
[8,] 0 0
[9,] 0 0
[10,] 0 0
[11,] 0 0
[12,] 0 0
[13,] 0 0
[14,] 0 0
[15,] 0 0
[16,] 0 0
[17,] 0 0
[18,] 0 0
[19,] 0 0
[20,] 0 03. Starting without forbidden ties in STERGM
Now we see if we can do something similar for a dynamic network, adding a simple Bernoulli dissolution model with mean relational duration of 100. We first try using vanilla ergm and stergm functionality, without EpiModel. We are fine with using the dissolution approximation for estimation, so we don’t need to do a new estimation at all. But we do need to define the dissolution parameter and the starting network for the simulation:
coef.diss <- log(100-1)
coef.form <- msm_mod_ergm$coef[[1]] - coef.diss
msm_sim_stergm_startnw <- simulate(msm_mod_ergm, nsim = 1,
control=control.simulate.ergm(MCMC.burnin = 1e5))
summary(msm_sim_stergm_startnw~edges+nodemix("role")) edges mix.role.I.I mix.role.I.R mix.role.R.R mix.role.I.V
3564 0 534 0 925
mix.role.R.V mix.role.V.V
1097 1008 Looks good. Now the dynamic simulation:
msm_sim_stergm <- simulate(msm_sim_stergm_startnw,
formation = ~edges + offset(nodematch('role', diff = TRUE, keep = 1:2)),
dissolution = ~edges,
monitor = "all",
coef.form = c(coef.form, -Inf, -Inf),
coef.diss = coef.diss,
time.slices = 1000,
output = c("stats")
)
plot(msm_sim_stergm)
Not the most elegant code ever, especially the bit about having to manually combine the fitted and offset coefficients. Nonetheless, it does the trick. There may be a more elegant way to do this using only standard stergm syntax, but I haven’t found it.
4. Starting without forbidden ties in EpiModel
Since this is a reasonably common type of task to want to do, we anticpated that there would be also be relatively streamlined syntax for it in an EpiModel context. Herein are some of our attempts.
First, we need to re-estimate the model using netest so that our fit will be in the right format:
msm_mod_epimodel <- netest(nw = msm_nw,
formation = ~ edges + offset(nodematch('role',
diff = TRUE, keep = 1:2)),
target.stats = nEdges,
coef.form = c(-Inf, -Inf),
coef.diss = dissolution_coefs(~offset(edges),
duration = 30, d.rate = 0)
)We are wise, so before we get to the actual disease sim, we will check netdx.
msm_dx_epimodel_no <- netdx(msm_mod_epimodel, nsteps = 100, nwstats.formula = ~nodemix('role'))
Network Diagnostics
-----------------------
- Simulating 1 networks
- Calculating formation statistics
- Calculating duration statistics
- Calculating dissolution statistics
plot(msm_dx_epimodel_no) 
Ugh oh, forbidden ties! Let’s boost up the burnin like we did in ergm.
msm_dx_epimodel_yes <- netdx(msm_mod_epimodel, nsteps = 100, nwstats.formula = ~edges+nodemix('role'),
set.control.ergm = control.simulate.ergm(MCMC.burnin = 1e5))
Network Diagnostics
-----------------------
- Simulating 1 networks
- Calculating formation statistics
- Calculating duration statistics
- Calculating dissolution statistics
plot(msm_dx_epimodel_yes)
That works. So this seems like it will indeed be easy in EpiModel. Now let’s make it work for an actual epidemic simulation using netsim. Lets just start with basic EpiModel syntax without changing any default controls relating to burnin or starting network:
msm_init <- init.net(i.num=100, r.num=0)
msm_control <- control.net('SIR', nsteps=100)
msm_param <- param.net(inf.prob=0.1, act.rate=0.3, rec.rate=.01)
msm_sim_epimodel <- netsim(msm_mod_epimodel, msm_param, msm_init, msm_control)
plot(msm_sim_epimodel, type='formation', ylim=c(1,100))
Not surprisingly there are forbidden ties, since we didn’t set the burnin or starting network yet.
Let’s try setting the burnin. We’ll use the same argument we used in netdx:
msm_control_take2 <- control.net('SIR', nsteps=100,
set.control.ergm =
control.simulate.ergm(MCMC.burnin = 1e5))
msm_sim_epimodel_take2 <- netsim(msm_mod_epimodel, msm_param, msm_init, msm_control_take2)
plot(msm_sim_epimodel_take2, type='formation', ylim=c(1,100)) 
That didn’t work. But now we’re in stergm land; maybe we need that argument instead. Indeed, reading ?control.net tells us that there is a set.control.stergm argument but not a set.control.ergm one.
msm_control_take3 <- control.net('SIR', nsteps=100,
set.control.stergm =
control.simulate.network(MCMC.burnin.min = 1e5))
msm_sim_epimodel_take3 <- netsim(msm_mod_epimodel, msm_param, msm_init, msm_control_take3)
plot(msm_sim_epimodel_take3, type='formation', ylim=c(1,100))
This seems like a dead-end. Let’s instead try setting the starting network directly. We could, of course, go in and write our own initialize.FUN function that allows us to do this. But again, since this seems like a basic type of user request that should flow relatively automatically from having fit a model with offsets, we would like to avoid that. How about we simply replace the network that is saved in our netest object with the starting network we created for the simulation using stergm syntax? The issue is that there are three networks stored in our netsim object; we need to figure out which one is right:
msm_mod_epimodel_2 <- msm_mod_epimodel
msm_mod_epimodel_2$fit$newnetwork <- msm_sim_stergm_startnw
msm_sim_epimodel_take4 <- netsim(msm_mod_epimodel_2, msm_param, msm_init, msm_control)
plot(msm_sim_epimodel_take4, type='formation', ylim=c(1,100))
Nope.
msm_mod_epimodel_2$fit$network <- msm_sim_stergm_startnw
msm_sim_epimodel_take5 <- netsim(msm_mod_epimodel_2, msm_param, msm_init, msm_control)
plot(msm_sim_epimodel_take5, type='formation', ylim=c(1,100))
Nope. Third time’s going to be a charm?
msm_mod_epimodel_2$fit$newnetworks[[1]] <- msm_sim_stergm_startnw
msm_sim_epimodel_take6 <- netsim(msm_mod_epimodel_2, msm_param, msm_init, msm_control)
plot(msm_sim_epimodel_take6, type='formation', ylim=c(1,100))
Still no. Hmmm. Indeed, diving into the code for netsim and its children using debug commands indicates that the network is called using the environment in which the estimation occured. This is presumably a very wise thing, but it would seem to preclude the option of overriding the starting network using standard EpiModel calls.
So - is the only option to create a new initalize.FUN? If so, it seems like a lot of work for users who have a basic and reasonable need: making sure there are no forbidden ties in the dynamic network simulated from a fit with offsets. And given that netdx was able to do this just fine, it seems like it should not be that hard to add into netsim. Yes?
f