Last updated: June 19, 2015 [//]: (Last updated 06-19-2015)

This tutorial is a joint product of the Statnet Development Team:

Mark S. Handcock (University of California, Los Angeles)
Carter T. Butts (University of California, Irvine)
David R. Hunter (Penn State University)
Steven M. Goodreau (University of Washington)
Skye Bender de-Moll (Oakland)
Pavel N. Krivitsky (University of Wollongong)
Martina Morris (University of Washington) For general questions and comments, please refer to statnet users group and mailing list
http://statnet.csde.washington.edu/statnet_users_group.shtml

1. Getting Started

Open an R session, and set your working directory to the location where you would like to save this work.

To install all of the packages in the statnet suite:

install.packages('statnet')
library(statnet)

Or, to only install the specific statnet packages needed for this tutorial:

install.packages('tergm') # will install the network package
install.packages('sna')

After the first time, to update the packages one can either repeat the commands above, or use:

update.packages('name.of.package')

For this tutorial, we will need one additional package (coda), which is recommended (but not required) by ergm:

install.packages('coda')

Make sure the packages are attached:

library(statnet)

or

library(tergm)
library(sna)
library(coda)

Check package version

# latest versions:  tergm 3.3.0, ergm 3.4.0, network 1.12.0, networkDynamic 0.7.1 (as of 6/22/2015)
sessionInfo()

Set seed for simulations – this is not necessary, but it ensures that we all get the same results (if we execute the same commands in the same order).

set.seed(0)

2. A quick review of static ERGMs

Exponential-family random graph models (ERGMs) represent a general class of models based in exponential-family theory for specifying the probability distribution underlying a set of random graphs or networks. Within this framework, one can—among other tasks—obtain maximum-likehood estimates for the parameters of a specified model for a given data set; simulate additional networks with the underlying probability distribution implied by that model; test individual models for goodness-of-fit, and perform various types of model comparison.

The basic expression for the ERGM class can be written as:

\(P(Y=y)=\frac{\exp(\theta'g(y))}{k(\theta)}\)


where:

This can be re-expressed in terms of the conditional log-odds of a single actor pair:

\(\operatorname{logit}{P(Y_{ij}=1 \mid y^{c}_{ij}) = \ln{\frac{P(Y_{ij}=1 \mid y^{c}_{ij})}{P(Y_{ij}=0 \mid y^{c}_{ij})}} = \theta'\delta(g(y))_{ij}}\)


where:

Fitting an ERGM usually begins with obtaining data:

data("florentine")
ls()
plot(flobusiness)

To refresh our memories on ERGM syntax, let us fit a cross-sectional example. Just by looking at the plot of flobusiness, we might guess that there are more triangles than expected by chance for a network of this size and density, and thus that there is some sort of explicit triangle closure effect going on. One useful way to model this effect in ERGMs that has been explored in the literature is with a gwesp statistic. When the alpha parameter of this statistic equals 0, then the statistic counts the number of edges in the network that are in at least one triangle. You can see this by summarizing this network and then counting them for yourself in the visualization:

summary(flobusiness~edges+gwesp(0, fixed=T))
        edges gwesp.fixed.0 
           15            12 
fit1 <- ergm(flobusiness~edges+gwesp(0,fixed=T))
Starting maximum likelihood estimation via MCMLE:
Iteration 1 of at most 20: 
The log-likelihood improved by 0.2804 
Step length converged once. Increasing MCMC sample size.
Iteration 2 of at most 20: 
The log-likelihood improved by 0.003476 
Step length converged twice. Stopping.

This model was fit using MCMC.  To examine model diagnostics and check for degeneracy, use the mcmc.diagnostics() function.
summary(fit1)

==========================
Summary of model fit
==========================

Formula:   flobusiness ~ edges + gwesp(0, fixed = T)

Iterations:  2 out of 20 

Monte Carlo MLE Results:
              Estimate Std. Error MCMC % p-value    
edges          -3.3792     0.6010      0 < 1e-04 ***
gwesp.fixed.0   1.5797     0.5674      0 0.00625 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

     Null Deviance: 166.4  on 120  degrees of freedom
 Residual Deviance:  78.1  on 118  degrees of freedom
 
AIC: 82.1    BIC: 87.67    (Smaller is better.) 

How to interpret the coefficients? Conditional on the rest of the network, a given edge’s log-odds of existing depends on how much its presence changes the number of edges in at least one triangle in the network.

With the estimation in place, we can simulate a new network from the given model:

sim1 <- simulate(fit1,nsim=1,
          control=control.simulate.ergm(MCMC.burnin=1000))
plot(sim1)

3. An Introduction to STERGMs (non-technical)

\(prevalence = incidence x duration\)


\(\ln{\frac{P(Y_{ij,t+1}=1 \mid y^{c}_{ij}, Y_{ij,t}=0)}{P(Y_{ij,t+1}=0\mid y^{c}_{ij}, Y_{ij,t}=0)}} = \theta^+ \delta(g^+(y))_{ij}\)


\(\ln{\frac{P(Y_{ij,t+1}=1 \mid y^{c}_{ij}, Y_{ij,t}=1)}{P(Y_{ij,t+1}=0\mid y^{c}_{ij}, Y_{ij,t}=1)}} = \theta^- \delta(g^-(y))_{ij}\)


4. An Introduction to STERGMs (more formal)

To understand STERGMs formally, we first review the ERGM framework for cross-sectional or static networks again, this time with more complete notation. Let \(\mathbb{Y}\subseteq \{1,\dotsc,n\}^2\) be the set of potential relations (dyads) among \(n\) nodes, ordered for directed networks and unordered for undirected. We can represent a network \(\mathbf{y}\) as a set of ties, with the set of possible sets of ties, \(\mathcal{Y}\subseteq 2^\mathbb{Y}\), being the sample space: \(\mathbf{y} \in \mathcal{Y}\). Let \(\mathbf{y}_{ij}\) be 1 if \((i,j)\in\mathbf{y}\) — a relation of interest exists from \(i\) to \(j\) — and 0 otherwise.

The network also has an associated covariate array \(\mathbf{X}\) containing attributes of the nodes, the dyads, or both. An exponential-family random graph model (ERGM) represents the pmf of \(\mathbf{Y}\) as a function of a \(p\)-vector of network statistics \(g(\mathbf{Y},\mathbf{X})\), with parameters \(\theta \in \mathbb{R}^p\), as follows:

\(\mbox{Pr}\left(\mathbf{Y} = \mathbf{y} \mid \mathbf{X}\right) = \frac{\exp\left\{\theta\cdot g(\mathbf{y}, \mathbf{X})\right\} } {k(\theta, \mathbf{X}, \mathcal{Y})},\)


where the normalizing constant \(k(\theta, \mathbf{X}, \mathcal{Y}) = \sum\limits_{\mathbf{y}' \in \mathcal{Y}}\exp\left\{\theta\cdot g(\mathbf{y}', \mathbf{X})\right\}\) is a summation over the space of possible networks on \(n\) nodes, \(\mathcal{Y}\). Where \(\mathcal{Y}\) and \(\mathbf{X}\) are held constant, as in a typical cross-sectional model, they may be suppressed in the notation. Here, on the other hand, the dependence on \(\mathcal{Y}\) and \(\mathbf{X}\) is made explicit.

In modeling the transition from a network \(\mathbf{Y}^{t}\) at time \(t\) to a network \(\mathbf{Y}^{t+1}\) at time \(t+1\), the separable temporal ERGM assumes that the formation and dissolution of ties occur independently from each other within each time step, with each half of the process modeled as an ERGM. For two networks (sets of ties) \(\mathbf{y},\mathbf{y}'\in \mathcal{Y}\), let \(\mathbf{y} \supseteq \mathbf{y}'\) if any tie present in \(\mathbf{y}'\) is also present in \(\mathbf{y}\). Define \(\mathcal{Y}^+(\mathbf{y})=\{\mathbf{y}'\in\mathcal{Y}:\mathbf{y}'\supseteq\mathbf{y}\}\), the networks that can be constructed by forming ties in \(\mathbf{y}\); and \(\mathcal{Y}^-(\mathbf{y})=\{\mathbf{y}'\in\mathcal{Y}:\mathbf{y}'\subseteq\mathbf{y}\}\), the networks that can be constructed dissolving ties in \(\mathbf{y}\).

Given \(\mathbf{y}^{t}\), a formation network \(\mathbf{Y}^+\) is generated from an ERGM controlled by a \(p\)-vector of formation parameters \(\theta^+\) and formation statistics \(g^+(\mathbf{y}^+, \mathbf{X})\), conditional on only adding ties:

\(\Pr\left(\mathbf{Y}^+ = \mathbf{y}^+\mid \mathbf{Y}^{t};\theta^+\right) = \frac{\exp\left\{\theta^+\cdot g^+(\mathbf{y}^+, \mathbf{X})\right\} } {k\left(\theta^+, \mathbf{X}, \mathcal{Y}^+(\mathbf{Y}^{t})\right)}, \quad \mathbf{y}^+ \in \mathcal{Y}^+(\mathbf{y}^{t})\).


A dissolution network \(\mathbf{Y}^-\) is simultaneously generated from an ERGM controlled by a (possibly different) \(q\)-vector of dissolution parameters \(\theta^-\) and corresponding statistics \(g^-(\mathbf{y}^-, \mathbf{X})\), conditional on only dissolving ties from \(\mathbf{y}^{t}\), as follows:

\(\Pr\left(\mathbf{Y}^- = \mathbf{y}^-\mid \mathbf{Y}^{t};\theta^-\right) = \frac{\exp\left\{\theta^-\cdot g^-(\mathbf{y}^-, \mathbf{X})\right\} } {k\left(\theta^-, \mathbf{X}, \mathcal{Y}^-(\mathbf{Y}^{t})\right)}, \quad \mathbf{y}^- \in \mathcal{Y}^-(\mathbf{y}^{t}).\)


The cross-sectional network at time \(t+1\) is then constructed by applying the changes in \(\mathbf{Y}^+\) and \(\mathbf{Y}^-\) to \(\mathbf{y}^{t}\), as follows: \(\mathbf{Y}^{t+1} = \mathbf{Y}^{t}\cup (\mathbf{Y}^+ - \mathbf{Y}^{t}) \, - \, ( \mathbf{Y}^{t} - \mathbf{Y}^-).\) which simplifies to either: \(\mathbf{Y}^{t+1} = \mathbf{Y}^+ - (\mathbf{Y}^{t} - \mathbf{Y}^-)\) \(\mathbf{Y}^{t+1} = \mathbf{Y}^-\cup (\mathbf{Y}^+ - \mathbf{Y}^{t})\)

Visually, we can sum this up as:


5. Notes on model specification and syntax

Within statnet, an ERGM involves one network and one set of network statistics, so these are specified together using R’s formula notation:

my.network \(\sim\) my.vector.of.model.terms


For a call to stergm, there is still one network (or a list of networks), but two formulas. These are now passed as three separate arguments: the network(s) (argument nw), the formation formula (argument formation), and the dissolution formula (argument dissolution). The latter two both take the form of a one-sided formula. E.g.:

        stergm(my.network,
            formation= ~edges+kstar(2),
            dissolution= ~edges+triangle
        )

There are other features of a call to either ergm or stergm that will be important; we list these features here, and illustrate them in one or more examples during the workshop.

dissolution = ~offset(edges)


            and the corresponding argument for passing the parameter value would be:

offset.coef.diss = 4.2


control=control.stergm(MCMC.burnin=10000)

            For a list of options, type ?control.stergm

6. Example 1: Multiple network cross-sections

One common form of data for modeling dynamic network processes consists of observations of network structure at two or more points in time on the same node set. Many classic network studies were of this type, and data of this form continue to be collected and analyzed.

Let us consider the famous Sampson monastery data:

data(samplk)
ls()
[1] "fit1"        "flobusiness" "flomarriage" "samplk1"     "samplk2"    
[6] "samplk3"     "sim1"       

To pass them into stergm, we need to combine them into a list:

samp <- list()
samp[[1]] <- samplk1
samp[[2]] <- samplk2
samp[[3]] <- samplk3

Now we must decide on a model to fit to them. Plotting one network:

plot(samplk1)