32 Epidemic Models with Demography
This tutorial demonstrates how to model a Susceptible-Infected (SI) epidemic in an open population. An example of an SI epidemic is HIV, with infected persons never leaving the infected stage before mortality. Modeling HIV is quite complicated, and this tutorial is not meant to be a full-scale HIV model: it does not include many components of HIV models, such as disease stage, drug therapy, and so on. This only represents SI dynamics generally with the possibility of vital dynamics. Although we use HIV as the running example, the vital-dynamics machinery shown here (arrivals, departures, and the coefficient adjustments that accompany them) is general to any open-population disease: it applies equally to, say, an endemic infection followed over many years, or a childhood infection in a population with ongoing births and deaths.
The population is open because we are now including arrivals and departures in the epidemic process. Simulating an epidemic in an open population requires adjustments that accommodate changes to the network structure over time, as detailed in the Module slides.
For simplicity and to keep the overall population size stable, our initial example will have the departure processes not depend on disease status (no disease-dependent mortality). We will then show the implications of the “death correction” to the estimation to demonstrate why we implement that model correction, and finally explore models in which the population size does vary.
Unlike in the tutorial in Basic Epidemics Module tutorial, we are using an attribute here but not the special group attribute. So we will allow for heterogeneity in network structure but assume homogeneous epidemic parameters. It is straightforward to extend EpiModel (we will cover this in NME-II) to allow for any arbitrarily specified attributes to impact both the network structure and epidemic processes.
Download the R script to follow along with this tutorial here.
Start with loading EpiModel:
32.1 Network Structure
The network is initialized as before, but now we set a vertex attribute of risk group onto the network. Risk group will be a binary variable, with the risk groups evenly sized in the population.
32.2 Model 1: Stable Population Size
We will proceed by fitting a network model in which Risk Group 1 has a higher mean degree and there is strong assortative (within-group) mixing (aka, homophily).
32.2.1 Network Parameterization
We will add heterogeneity in mean degree by risk group, and homophily in risk group mixing, using nodefactor and nodematch terms and their associated target statistics.
The formula for the formation component of the model will be as follows. The nodefactor term represents variation in the overall mean degree by the named attribute, risk. The nodematch represents mixing by that named attribute.
The overall mean degree will be 0.5, but will vary by risk group. We will specify that the mean degree for risk group 1 is 0.75. Because the overall mean degree is 0.5, and the proportions of the risk groups are equal, that implies that the mean degree for risk group 0 will be 0.25.
The target statistic for the edges term is the overall mean degree times the population size, divided by two because each edge involves two nodes: 0.5 * 500 / 2 = 125. The target statistic for the nodefactor term is the product of the mean degree of risk group 1 and the size of risk group 1: 0.75 * 250 = 187.5. Finally, we will specify that 90% of the partnerships occur between persons of the same risk group: 0.9 * 125 = 112.5 edges.
The dissolution model components will be as follows. We specify a homogeneous dissolution model with an average relational duration of 40 weeks. The dissolution coefficients must be adjusted to accommodate the exogenous probability of edge dissolution due to mortality (or other forms of nodal departure) – this is the “death correction.” As shown in the Module slides, this adjustment has the impact of increasing the coefficient. Recall that this is because the dissolution model is actually implemented as a persistence model: a larger coefficient = longer partnerships.
Code
Dissolution Coefficients
=======================
Dissolution Model: ~offset(edges)
Target Statistics: 40
Crude Coefficient: 3.663562
Mortality/Exit Rate: 0.001
Adjusted Coefficient: 3.7469The d.rate parameter here is the average mortality (or other departure) rate per time step. This will correspond to the epidemic parameters specified below in param.net. In the case that there is disease-induced mortality, this parameter should represent a weighted average of the different departure rates in the full system (we will explore this below).
In this tutorial example a d.rate parameter of 0.001 corresponds to an average “lifespan” (time spent in the population) of 1/0.001 = 1000 time steps. This is probably unrealistic for most disease models (especially since the average relations last 40 time steps), but we are using a small value of lifespan for demonstration purposes. In fact, you may encounter the following message if the departure rate is too high relative to the average relational duration.
In this test case below, the average “lifespan” is half as long as the average relationship (which probably doesn’t make any sense). We have tried to catch this problem for you with an informative error message.
Ok, so let’s use our smaller d.rate parameter of 0.001.
32.2.2 Estimation and Diagnostics
The network model is estimated using the netest function with all the inputs we specified above.
We can use the summary function to display the coefficient estimates for the nodefactor and nodematch terms, which are strong (> 0 on the log odds scale) and highly significant. The specific values of those coefficients do not matter much here (what matters is how they generate simulations of the complete network, below), but the positive sign of both terms shows these are greater than expected under the null model, precisely as we have parameterized them.
Call:
ergm(formula = formation, constraints = constraints, offset.coef = coef.form,
target.stats = target.stats, eval.loglik = FALSE, control = set.control.ergm,
verbose = verbose, basis = nw)
Monte Carlo Maximum Likelihood Results:
Estimate Std. Error MCMC % z value Pr(>|z|)
edges -9.1616 0.2974 0 -30.810 <1e-04 ***
nodefactor.risk.1 0.6389 0.1186 0 5.389 <1e-04 ***
nodematch.risk 2.0261 0.2868 0 7.064 <1e-04 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Log-likelihood was not estimated for this fit. To get deviances, AIC, and/or BIC, use ‘*fit* <-logLik(*fit*, add=TRUE)’ to add it to the object or rerun this function with eval.loglik=TRUE.
Dissolution Coefficients
=======================
Dissolution Model: ~offset(edges)
Target Statistics: 40
Crude Coefficient: 3.663562
Mortality/Exit Rate: 0.001
Adjusted Coefficient: 3.7469For the diagnostics for this model, we set the nwstats.formula to expand the formation formula to specify that we want to monitor the mean degree and mixing by the risk attribute. We do this with nodefactor and nodematch terms, respectively. By adding them to the nwstats.formula, but not in the formation formula for the TERGM itself, we can get summary statistics for simulations to see what the expected values of these terms are for this model. Setting levels = NULL in the nodefactor term will output statistics for both risk groups in the model.
We can also monitor the nodematch term with a diff = TRUE argument, which allows for differential homophily by group. In this model, there would never be a reason to specify differential homophily by group in the ERGM formula because it would be over-specified: because risk is a binary attribute, the two within-group (diagonal) cells of the mixing matrix are already determined by the edges, nodefactor, and overall nodematch terms together. Differential homophily adds information only for models without a nodefactor term, or for models whose mixing attribute has more than two levels. But again we show it here for comparison, and getting you familiar with the different syntax options.
Printing the diagnostics shows a good fit between the targets and the simulated outputs.
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)
edges 125.0 124.554 -0.356 0.887 -0.502 2.882
nodefactor.risk.0 NA 61.498 NA 0.904 NA 2.204
nodefactor.risk.1 187.5 187.611 0.059 1.528 0.073 5.174
nodematch.risk.0 NA 24.698 NA 0.394 NA 1.055
nodematch.risk.1 NA 87.755 NA 0.767 NA 2.748
SD(Statistic)
edges 11.097
nodefactor.risk.0 10.845
nodefactor.risk.1 19.403
nodematch.risk.0 4.977
nodematch.risk.1 9.639
Duration Diagnostics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 40 39.514 -1.214 0.259 -1.872 0.686 3.366
Dissolution Diagnostics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges 0.025 0.025 0.905 0 1.626 0 0.014Plotting the diagnostics shows that the simulations are right on target for the terms in the formation model. We will skip plotting the dissolution model diagnostics for this example.
Recall that this TERGM estimation and diagnostics has not yet been confronted with the complexities of an open population in which there is arrival and departure. So we should not expect any challenges with corrections to any changing population or mortality, because the population size and composition remains fixed during this step. Only when we move on to the next step of the epidemic simulation do these exogenous processes introduce these factors.
32.2.3 Epidemic Simulation
This epidemic simulation will investigate how variation in mean degree by risk group and highly assortative risk group mixing jointly impact epidemic prevalence overall and by risk group. We will use one set of epidemic parameters, initial conditions, and control settings shared across the two counterfactual models.
The parameters are entered into the param.net function should be familiar by now. Since this is an SI epidemic, there will be no recovery rate parameter. The three new parameters will control the arrival rate, departure rate among the susceptibles, and departure rate in the infected group. Note that this is the same departure rate as what we specified in dissolution_coefs. For this step, we will assume that all three rates are equal, which implies that there is no disease-induced mortality and a stable population size. These departure rates and the d.rate in dissolution_coefs describe the same demographic process from two sides: the parameters here govern who actually leaves the population during the epidemic, while d.rate tells the network model how much extra edge persistence to build in (the death correction) to offset the ties those departures would otherwise break. The two must therefore be kept consistent.
Similar to the previous tutorial, we only need to specify the number infected at the outset. To generate stable epidemic conditions quickly, we set the initial prevalence to 10% (50 of the 500 nodes), distributed randomly among the population.
The inputs for the control settings share some similarities with the control settings from the Basic EpiModel tutorials, but we will be adding some new parameters here:
resimulate.network: whether the network is re-drawn at every time step (TRUE) rather than simulated once as a fixed dynamic timeline (FALSE, the default). You need it whenever the network itself must change during the run: here because arrivals and departures add and remove nodes, and in the next tutorial because a formation term depends on an attribute (disease status) that changes over time. A closed model with a static formula can leave itFALSE; an open population cannot, so we set it toTRUE.epi.by: a categorical attribute by which to break out the prevalence and flow outputs (here"risk"), so we can compare epidemic outcomes across groups. It is the general-attribute counterpart to the automatic by-group output that the specialgroupattribute provides. By default it isNULL(no stratification).tergmLite: a streamlined, memory-light network representation that speeds simulation 20-50 fold. Use it for the large or long runs typical of epidemic models, when you do not need the full network history. The tradeoff is exactly that: because only a sparse representation is retained, you cannot reconstruct or plot the network at past time points. It defaults toFALSE(full network retained); otherwise the two are equivalent.
For the main epidemic model comparison, we are going to run 20 simulations on 5 cores with tergmLite = TRUE. Note that we are using the nwstats.formula argument to add a couple of network terms to monitor for the post-simulation diagnostics.
To simulate the epidemic model, we again use the netsim function.
32.2.4 Epidemic Model Analysis
Our model analysis will consist of model diagnostics, and then the examination of the epidemiological outcomes overall and by risk group.
32.2.4.1 Post-Simulation Diagnostics
It is important to first examine the network model diagnostics after the epidemic simulation, since the vital dynamics within these simulations may have changed the structure of the network in unexpected ways.
Printing the netsim output object shows its contents. The output compartments now include the post-simulation diagnostics, similar to the output from a netdx object (pre-simulation diagnostics).
EpiModel Simulation
=======================
Model class: netsim
Simulation Summary
-----------------------
Model type: SI
No. simulations: 20
No. time steps: 300
No. NW groups: 1
Fixed Parameters
---------------------------
inf.prob = 0.1
act.rate = 5
a.rate = 0.001
ds.rate = 0.001
di.rate = 0.001
groups = 1
Model Output
-----------------------
Variables: s.num s.num.risk0 s.num.risk1 i.num i.num.risk0
i.num.risk1 num num.risk0 num.risk1 si.flow ds.flow di.flow
a.flow
Networks: sim1 ... sim20
Transmissions: sim1 ... sim20
Formation Statistics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means)
edges 125.0 128.250 2.600 1.179 2.757 5.321
nodefactor.risk.0 NA 64.230 NA 0.894 NA 3.830
nodefactor.risk.1 187.5 192.270 2.544 1.800 2.651 9.096
nodematch.risk.0 NA 25.678 NA 0.421 NA 1.838
nodematch.risk.1 NA 89.698 NA 0.893 NA 4.157
meandeg NA 0.510 NA 0.005 NA 0.021
SD(Statistic)
edges 12.471
nodefactor.risk.0 10.498
nodefactor.risk.1 20.849
nodematch.risk.0 4.948
nodematch.risk.1 10.137
meandeg 0.049
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)Note that even with the introduction of vital dynamics (births and deaths), the simulation means are very close to the target statistics. While there may be a small difference here, it is quite small on the absolute scale (3-5 edges per time step), which is due to the smaller sample size (smaller network and low number of simulations) rather than a model bias. Note that dissolution statistics are not currently available in this model, but this is something we plan to add in the future.
Here we show how to plot model diagnostics recorded within the netsim object. Although we have saved all the network statistics in the formation formula, here we plot the number of edges over time. We are happy with the diagnostics because the interquartile range of the simulations overlaps with the target statistic lines.
We can also pull out specific formula terms, like edges and meandeg (for mean degree), to plot separately here. We also show the unique values of each of the 20 simulations.
Code
Overall, these diagnostics suggest that any of the exogenous vital dynamic processes that we added to the model did not have unintended consequences on the network structure. Good news!
32.2.4.2 Demographic Outcomes
Next, let’s explore the demographic outcomes that are now part of the epidemic model processes. First, we can look at the overall population size, num, and the population size by risk group attribute. Because we have balanced birth and death rates, we expect the total population size to remain stable at 500 (our initial network size); it does! We also expect the 50/50 distribution of the two groups to remain stable; it does!
Code
Note however that we are evaluating stability based on the means across simulations. Some individual outlier simulations have the population changing +/- 30 nodes, and the same goes for the risk group-specific numbers. That’s ok.
Next, we can look at the flows for the three vital dynamics processes. a.flow is the count of new arrivals into the model, ds.flow is the number of departures from the susceptible group, and di.flow is the number of departures from the infected group.
Code
The overall count for the arrivals is stable, and right around the expected value of 0.5 per time step (because the per capita arrival rate is 0.001 and the population size is 500, then 500 * 0.001 = 0.5). The two death flow counts are changing. Why? Because (as we will see below), the number of susceptibles is shrinking as the number of infecteds is growing.
The raw departure counts are hard to read while the compartments themselves are changing size: a falling ds.flow could mean a lower per-capita death rate, or simply that fewer susceptibles remain to die. To separate the per-capita risk from the size of the pool, we compute standardized death rates, dividing each departure flow by the current size of its compartment, and add them to the model object using the mutate_epi function. If the underlying per-capita departure rates are truly constant, these standardized rates should be flat and equal.
With that added, let’s plot it. Now, the numbers converge to stable and equal values (however there is some fluctuation at the beginning of the model due to low numbers of infecteds).
32.2.4.3 Epidemic Outcomes
Finally, let’s take a look at the epidemic outcomes for disease prevalence. Overall, total population prevalence is contained in i.num, and we can plot that here.
The next plot shows the prevalence by risk group. Our model demonstrates that Risk Group 1, with the higher mean degree and strong assortative mixing, has the higher disease prevalence.
32.3 Model 2: Impact of Death Correction
To clearly demonstrate the importance of the “death correction” for the dissolution component of the model, we will re-estimate the network model from above, but incorrectly setting the d.rate parameter in dissolution_coefs to 0 (whereas it was properly specified before as the weighted average of the death rates in the epidemic parameters). Notice that the crude and adjusted coefficients are the same.
Code
Dissolution Coefficients
=======================
Dissolution Model: ~offset(edges)
Target Statistics: 40
Crude Coefficient: 3.663562
Mortality/Exit Rate: 0
Adjusted Coefficient: 3.663562We next reestimate the network model (new diagnostics with netdx are not needed because the death correction does not impact the TERGM estimation), and then resimulate the epidemic model with netsim (all the other inputs are recycled from Model 1).
Let’s first review the demographics of the model as we did above, looking at the total population size and the group-specific population size. Both are stable, so any changes to the network stats should not be influenced by population size variation.
Code
Next, we plot the edges in the post-simulation network diagnostics. Now there is a clear downward trend in the count of edges.
Code
If the population size declined, we might expect the total edge count to decline, but the mean degree (meandeg) should remain constant (because of that network-size correction). Alas, the mean degree declined too.
Code
The numerical summary of the statistics in the print of the netsim object reveals that the simulations for edges are a bit lower than the target, and the meandeg is lower than the expected 0.5.
EpiModel Simulation
=======================
Model class: netsim
Simulation Summary
-----------------------
Model type: SI
No. simulations: 20
No. time steps: 300
No. NW groups: 1
Fixed Parameters
---------------------------
inf.prob = 0.1
act.rate = 5
a.rate = 0.001
ds.rate = 0.001
di.rate = 0.001
groups = 1
Model Output
-----------------------
Variables: s.num s.num.risk0 s.num.risk1 i.num i.num.risk0
i.num.risk1 num num.risk0 num.risk1 si.flow ds.flow di.flow
a.flow
Networks: sim1 ... sim20
Transmissions: sim1 ... sim20
Formation Statistics
-----------------------
Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means)
edges 125.0 117.504 -5.997 1.008 -7.434 5.825
nodefactor.risk.0 NA 59.836 NA 0.952 NA 5.722
nodefactor.risk.1 187.5 175.171 -6.576 1.323 -9.322 7.709
nodematch.risk.0 NA 24.253 NA 0.466 NA 2.990
nodematch.risk.1 NA 81.920 NA 0.644 NA 3.998
meandeg NA 0.472 NA 0.004 NA 0.021
SD(Statistic)
edges 11.727
nodefactor.risk.0 11.408
nodefactor.risk.1 17.595
nodematch.risk.0 5.474
nodematch.risk.1 8.667
meandeg 0.045
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)This example highlights why we add that death rate correction in dissolution_coefs: to increase the persistence of edges from natural dissolution to account for the competing risk of edge dissolution from departure (death).
32.4 Model 3: Varying Population Size
Finally, let’s build a model in which there is a higher rate of disease-specific mortality. Specifically, we will assume here that the mortality rate from disease is twice that of the natural mortality rate (due to all other causes). Therefore, we will parameterize a di.rate = 0.002 and a ds.rate = 0.001. We will not change any other epidemic parameters.
Going back to the network model parameterization, no changes are needed for the formation model target statistics. However, we do need to update the d.rate parameter in dissolution_coefs to account for the higher average mortality in this new model parameterization.
What value of d.rate should we choose, since there are now two death rates in param.net? One natural starting point would be 0.0015, the average of the two rates above. However, we also need to consider the weighted average of person-time spent in each disease state. If there were equal time spent in S and I, then the simple average would do. But in our model, more total person-time is spent in the infected state.
In practice, coming up with a d.rate that works here involves some trial and error, guided by one specific post-simulation diagnostic: the meandeg (mean degree) plot. If the death correction matches the realized flow of departures, the mean degree holds flat at its target; if d.rate is too low, edges are under-corrected and the mean degree drifts downward.
Let us show one round. We start with the simple average, d.rate = 0.0015, then estimate and simulate.
Plotting the mean degree from this first attempt shows that it still drifts below the target of 0.5: the simple average under-corrects, because more person-time is spent in the higher-mortality I state than an equal split would imply.
Code
Raising d.rate and re-checking this same diagnostic, we settle on d.rate = 0.0018, which holds the mean degree flat (as we confirm at the end of this section).
Code
Dissolution Coefficients
=======================
Dissolution Model: ~offset(edges)
Target Statistics: 40
Crude Coefficient: 3.663562
Mortality/Exit Rate: 0.0018
Adjusted Coefficient: 3.818895Next, we estimate the network model and simulate the epidemic model again.
First, we review the total population size. Now there is a pretty steep decline in the population size as prevalence increases.
Code
What is the impact on the edges count in the network simulations? It also declines as the population size declines. That is OK! The edges are a function of the mean degree and the population size (N), so as N decreases the edges should decrease.
Code
This is where the second coefficient adjustment enters. As the population size N changes during the simulation, netsim automatically applies a network-size correction: it shifts the edges coefficient by a log(N)-based offset at each time step so that the model holds the mean degree (edges per capita) fixed, rather than the raw edge count. This is the counterpart to the death correction introduced at the start of the module, and it is why the edge count above is free to fall with the population while the mean degree is not. If both the network-size correction and the death correction are working, the mean degree (meandeg) should remain constant. And voila! This is a good model diagnostic.
