28 Lab

This lab on network model parameterization follows the introduction and description in the slides in this module.

Note

Download the R script to follow along with this tutorial here.

First, start by loading the EpiModel package.

Code

library(EpiModel)

28.1 Model Initialization

Start by initializing a network of 2000 persons in the population.

Code

mynet <- network_initialize(2000)

Add an attribute for attribute (1000 females, 1000 males). Note we call it group for the reasons mentioned in Module 5 that it is a “special” attribute in EpiModel.

Code

sex <- c(rep(1, 1000), rep(2, 1000))
mynet <- set_vertex_attribute(mynet, "group", sex)

Next, we add a community attribute (800 cmty = 1; 1200 cmty = 2; split evenly across sex).

Code

cmty <- c(rep(1, 400), rep(2, 600), rep(1, 400), rep(2, 600))
mynet <- set_vertex_attribute(mynet, "cmty", cmty)

Use the table function to confirm that the attribute combinations are as expected.

Code

table(get_vertex_attribute(mynet, "group"),
      get_vertex_attribute(mynet, "cmty"))

28.2 Fitting the Network Model

Following the terminology and syntax demonstrated previously, we would like to fit a TERGM where the formula includes terms for overall mean degree (edges), mixing by community (nodematch("cmty", diff = FALSE)), a restriction on the maximum total degree (degrange(from = 3)), having a degree exactly of one by sex (aka monogamy, degree(1, by = "group")), and a prohibition of mixing by group (nodematch("group", diff = FALSE)).

Code

formation <- ~edges + nodematch("cmty", diff = FALSE) + degrange(from = 3) +
        degree(1, by = "group") + nodematch("group", diff = FALSE)

The target statistics associated with these model terms are as follows (see the slides for more details).

Code

target.stats <- c(850, 600, 0, 550, 350, 0)

We then fit the TERGM with the above formula and data, plus a parameterization of a dissolution model with a mean relational duration of 60 time steps.

Code

myfit <- netest(mynet,
                formation = formation,
                target.stats = target.stats,
                coef.diss = dissolution_coefs(~offset(edges), 60))

Finally in this step, we run model diagnostics to ensure the fit is good.

Code

mydx <- netdx(myfit, nsims = 25, nsteps = 500, ncores = 5)

A printing of the mydx object.

Code

mydx

EpiModel Network Diagnostics
=======================
Diagnostic Method: Dynamic
Simulations: 25
Time Steps per Sim: 500

Formation Diagnostics
----------------------- 
                Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means)
edges              850  840.055   -1.170  1.123  -8.856         4.601
nodematch.cmty     600  595.077   -0.821  1.218  -4.043         6.087
deg3+                0    0.000      NaN    NaN     NaN         0.000
deg1.group1        550  545.042   -0.902  0.842  -5.887         4.024
deg1.group2        350  352.886    0.825  0.673   4.289         3.672
nodematch.group      0    0.000      NaN    NaN     NaN         0.000
                SD(Statistic)
edges                  19.308
nodematch.cmty         19.067
deg3+                   0.000
deg1.group1            16.519
deg1.group2            14.728
nodematch.group         0.000

Duration Diagnostics
----------------------- 
      Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges     60   60.103    0.172  0.156   0.663         0.672         2.061

Dissolution Diagnostics
----------------------- 
      Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges  0.017    0.017    0.221      0   0.948             0         0.004

There is a small difference between the simulations and data (target stats). We can also confirm with a plot.

Code

plot(mydx)

If it is of interest, you can also extract the raw data that goes into producing these numerical summaries and plots.

Code

head(get_nwstats(mydx), 10)

   time sim edges nodematch.cmty deg3+ deg1.group1 deg1.group2 nodematch.group
1     1   1   847            569     0         559         365               0
2     2   1   845            570     0         555         357               0
3     3   1   841            567     0         553         357               0
4     4   1   838            561     0         548         366               0
5     5   1   845            569     0         547         361               0
6     6   1   840            566     0         548         372               0
7     7   1   832            561     0         550         368               0
8     8   1   831            558     0         543         371               0
9     9   1   833            553     0         543         377               0
10   10   1   837            558     0         545         383               0

28.3 Epidemic Simulation

To get the epidemic simulation started, we first need to parameterize the model with epidemic model parameters, initial conditions, and control settings (details provided here).

Code

myparam <- param.net(inf.prob = 0.6, inf.prob.g2 = 0.6,
                     act.rate = 1.8,
                     rec.rate = 0.1, rec.rate.g2 = 0.1)
myinit <- init.net(i.num = 100, i.num.g2 = 100)
mycontrol <- control.net("SIS", nsteps = 50, nsims = 5, 
                         nwstats.formula = ~edges + nodematch("cmty") + 
                         degree(0:5, by = "group"), verbose = TRUE)

Run the epidemic simulation with netsim.

Code

mySIS <- netsim(myfit, param = myparam, init = myinit, control = mycontrol)

Plot infection prevalence by group. Note that the prevalence in group 1 (females) is higher than group 2 (males).

Code

plot(mySIS, y = c("i.num", "i.num.g2"), legend = TRUE)

Here are numerical summaries for the sex-specific disease prevalence at the final time step.

Code

mySIS <- mutate_epi(mySIS, prev.f = i.num / num,
                           prev.m = i.num.g2 / num.g2)
df <- as.data.frame(mySIS, out = "mean")
df$prev.f[df$time == 50]

[1] 0.204

Code

df$prev.m[df$time == 50]

[1] 0.1752

We can also calculate and plot network stats as a time series from the epidemic simulations (to ensure that the epidemic simulations did not have unintended consequences on the network structure).

Code

plot(mySIS, type = "formation", qnts = FALSE, sim.lines = TRUE)

--- output: html_document editor_options: chunk_output_type: console --- # Lab ```{r setup, include = FALSE} library(knitr) opts_chunk$set(comment = NA) par(mar = c(3, 3, 1, 1), mgp = c(2, 1, 0)) ``` This lab on network model parameterization follows the introduction and description in the [slides](mod6-Intro.html) in this module. ::: {.callout-note} Download the R script to follow along with this tutorial [here](mod6-Lab.R). ::: First, start by loading the `EpiModel` package. ```{r, message = FALSE, warning=FALSE} library(EpiModel) ``` ## Model Initialization Start by initializing a network of 2000 persons in the population. ```{r} mynet <- network_initialize(2000) ``` Add an attribute for attribute (1000 females, 1000 males). Note we call it `group` for the reasons mentioned in [Module 5](../5_nodal_attr/mod5-NodalAttributes.html) that it is a "special" attribute in EpiModel. ```{r} sex <- c(rep(1, 1000), rep(2, 1000)) mynet <- set_vertex_attribute(mynet, "group", sex) ``` Next, we add a community attribute (800 `cmty` = 1; 1200 `cmty` = 2; split evenly across sex). ```{r} cmty <- c(rep(1, 400), rep(2, 600), rep(1, 400), rep(2, 600)) mynet <- set_vertex_attribute(mynet, "cmty", cmty) ``` Use the `table` function to confirm that the attribute combinations are as expected. ```{r} table(get_vertex_attribute(mynet, "group"), get_vertex_attribute(mynet, "cmty")) ``` ## Fitting the Network Model Following the terminology and syntax demonstrated previously, we would like to fit a TERGM where the formula includes terms for overall mean degree (`edges`), mixing by community (`nodematch("cmty", diff = FALSE)`), a restriction on the maximum total degree (`degrange(from = 3)`), having a degree exactly of one by sex (aka monogamy, `degree(1, by = "group")`), and a prohibition of mixing by group (`nodematch("group", diff = FALSE)`). ```{r} formation <- ~edges + nodematch("cmty", diff = FALSE) + degrange(from = 3) + degree(1, by = "group") + nodematch("group", diff = FALSE) ``` The target statistics associated with these model terms are as follows (see the [slides](mod6-Intro.html) for more details). ```{r} target.stats <- c(850, 600, 0, 550, 350, 0) ``` We then fit the TERGM with the above formula and data, plus a parameterization of a dissolution model with a mean relational duration of 60 time steps. ```{r, results = "hide", message=FALSE, warning=FALSE} myfit <- netest(mynet, formation = formation, target.stats = target.stats, coef.diss = dissolution_coefs(~offset(edges), 60)) ``` Finally in this step, we run model diagnostics to ensure the fit is good. ```{r, results = "hide", message=FALSE, warning=FALSE} mydx <- netdx(myfit, nsims = 25, nsteps = 500, ncores = 5) ``` A printing of the `mydx` object. ```{r} mydx ``` There is a small difference between the simulations and data (target stats). We can also confirm with a plot. ```{r} plot(mydx) ``` If it is of interest, you can also extract the raw data that goes into producing these numerical summaries and plots. ```{r} head(get_nwstats(mydx), 10) ``` ## Epidemic Simulation To get the epidemic simulation started, we first need to parameterize the model with epidemic model parameters, initial conditions, and control settings (details provided [here](../4_basic_epimodel/mod4-BasicSIS.html)). ```{r} myparam <- param.net(inf.prob = 0.6, inf.prob.g2 = 0.6, act.rate = 1.8, rec.rate = 0.1, rec.rate.g2 = 0.1) myinit <- init.net(i.num = 100, i.num.g2 = 100) mycontrol <- control.net("SIS", nsteps = 50, nsims = 5, nwstats.formula = ~edges + nodematch("cmty") + degree(0:5, by = "group"), verbose = TRUE) ``` Run the epidemic simulation with `netsim`. ```{r, results = "hide", message=FALSE, warning=FALSE} mySIS <- netsim(myfit, param = myparam, init = myinit, control = mycontrol) ``` Plot infection prevalence by group. Note that the prevalence in group 1 (females) is higher than group 2 (males). ```{r} plot(mySIS, y = c("i.num", "i.num.g2"), legend = TRUE) ``` Here are numerical summaries for the sex-specific disease prevalence at the final time step. ```{r} mySIS <- mutate_epi(mySIS, prev.f = i.num / num, prev.m = i.num.g2 / num.g2) df <- as.data.frame(mySIS, out = "mean") df$prev.f[df$time == 50] df$prev.m[df$time == 50] ``` We can also calculate and plot network stats as a time series from the epidemic simulations (to ensure that the epidemic simulations did not have unintended consequences on the network structure). ```{r} plot(mySIS, type = "formation", qnts = FALSE, sim.lines = TRUE) ```