1. Preparation

First, we install the required packages.

• specify the settings that work for your coding environment)
• these packages are already installed if you are using the binder server

# install.packages('network', .libPaths(), repos='http://cran.us.r-project.org')
# install.packages('ergm', .libPaths(), repos='http://cran.us.r-project.org')
# install.packages('devtools', .libPaths(), repos='http://cran.us.r-project.org')

We have a workshop package that contains the data we will be working with and some utility functions.

devtools::install_github('tedhchen/ergmWorkshopTools')

Next, we load the libraries.

library(ergm)
library(ergmWorkshopTools)

As you can see from the output, loading the ergm package also loaded the network package.

There are also these references to the statnet project. They are the main developers of ERGM functionalities in R.

2. Building our networks

First, load the data.

The data sets we will be using are part of the workshop package. We start with the militarized interstate dispute data from Zeev Maoz.

data(mid_mat)
?mid_mat

checkmat() is a utility function that plots the supplied adjacency matrix.

Note about plotting: You can ignore the plotting options if you are not using a jupyter notebook.

options(repr.plot.width=10, repr.plot.height=10)
checkmat(mid_mat)

Next, we use the network function to create a network object.

?network

nw <- network(mid_mat, directed = F)
nw

options(repr.plot.width=10, repr.plot.height=10)
par(mfrow = c(1, 1))
set.seed(210615); plot(nw)

There are different ways network data is usually stored.

The adjacency matrix format is generally the easiest to work with, but it can take a lot of space.

Sometimes, we will have edgelists instead. In fact, unless you have pre-prepared network data, edgelists are usually what you'll have.

We need to be careful when working with them.

data(mid_edgelist)
nw.b <- network(mid_edgelist, directed = F)
nw.b

What's the difference?

options(repr.plot.width=20, repr.plot.height=10)
par(mfrow = c(1, 2))
set.seed(210615); plot(nw); plot(nw.b)

Where are the isolates?

The edgelist does not contain information about them.

When working with this kind of data, we have to make sure we have a list of the full actors in the network.

data(mid_node_attr)
head(mid_node_attr)
?mid_node_attr

Let's manually make the adjacency matrix from the edgelist.

• Start by making an empty matrix.

mat <- matrix(0, ncol = 189, nrow = 189, dimnames = list(row.names(mid_node_attr), row.names(mid_node_attr)))
mat[1:10, 1:10]

• Next, loop through the edgelist and fill the correct cells.

for(i in 1:nrow(mid_edgelist)){
    mat[mid_edgelist[i, 1], mid_edgelist[i, 2]] <- 1
    mat[mid_edgelist[i, 2], mid_edgelist[i, 1]] <- 1
}

options(repr.plot.width=14, repr.plot.height=7)
par(mfrow = c(1, 2))
checkmat(mat);checkmat(mid_mat)

Let's get back to working with our network.

nw

Let's work with vertex attributes.

The %v% operator is how to access vertex/node variables from the network object.

nw%v%'vertex.names'

We can look at how democracies and nondemocracies differ in conflict behavior.

nw%v%'dem' <- mid_node_attr$'dem'
nw

nw%v%'dem'

Red nodes are democracies and black are nondemocracies.

options(repr.plot.width=10, repr.plot.height=10)
set.seed(210615); plot(nw, vertex.col = ifelse(nw%v%'dem', 2, 1)); legend('bottomright', legend = c('Democracy', 'Nondemocracy'), fill = c(2, 1), bty = 'n', cex = 1.1)

Edge attributes can also be specified in a similar way but it's easier to just use a separate matrix.

These are the geographical contiguity and joint-democracy networks.

data(contig); data(joint_dem)

options(repr.plot.width=16, repr.plot.height=8)
par(mfrow = c(1, 2))
set.seed(210615); plot(network(contig, directed = F), sub = 'Contiguity', cex.sub = 2); plot(network(joint_dem, directed = F), sub = 'Joint-Democracy', cex.sub = 2)

3. Specifying the ERGM

Let's take a look at the basic function for ERGM fitting: ergm.

ergm(network ~ ergm-terms)

What are the ERGM terms that we can use?

3.1 ERGM Terms

• use ?`ergm-terms' to look at all the existing terms in the ergm package
• the search.ergmTerms function also works well

?`ergm-terms`

search.ergmTerms(keyword = 'transitive')

3.2 Starting to fit an ERGM

We start with some node and dyad variables.

• node-level democractic regime nodefactor('dem')
• joint-democracies edgecov(joint_dem)

m0 <- ergm(nw ~ edges + nodefactor('dem') + edgecov(joint_dem))
summary(m0)

Let's compare it with the logistic regression.

dyad_df <- ergmMPLE(nw ~ edges + nodefactor('dem') + edgecov(joint_dem))
dyad_df$predictor

summary(glm(dyad_df$response ~ dyad_df$predictor - 1, weights = dyad_df$weights, family = 'binomial'))

They are exactly the same.

• when there are no network effects, the ERGM and a logistic regression is equivalent

Let's move to some network effects.

• start with the "pile-on" effect -kstar(2) (which is the 2-star term) is the basic term for this

m1 <- ergm(nw ~ edges + nodefactor('dem') + edgecov(joint_dem) + kstar(2), control = control.ergm(seed = 210616))

What happened here?

• This is one of the difficulties of ERGMs: they run into degeneracy issues.
• Degeneracy is a problem that arises during the simulation-based model estimation procedure.
• For simple model terms (e.g. 2-star), a single edge contribute to too many of the corresponding local configurations.
• We only end up simulating very dense or very sparse networks (i.e. the model is not a good representation of our network).

So instead, we can use alternative network statistics designed to address these problems.

• For the preferential attachment (or pile-on) effect, the commonly-used one is the geometrically-weighted degree.
• These kind of "geometrically weighted" terms are used to down-weigh the contribution of each edge to network configurations.
• In the ergm package, it is called gwdegree.

References:

Hunter, DR et al. 2008. "ergm: A package to fit, simulate and diagnose exponentiall-family models for networks." Journal of Statistical Software.

Snijders, T et al. 2006. "New specifications for exponential random graph models." Sociological Methodology.

m1 <- ergm(nw ~ edges + nodefactor('dem') + edgecov(joint_dem) + gwdegree(decay = 0, fixed = T), control = control.ergm(seed = 210616))

Looks like it converged.

summary(m1)

• the gwdeg.fixed.0 term has coefficient of -1.16
• for the geometrically weighted degree term, negative means popular
  • This is an easy mistake: a review in 2016 shows that only 4/16 papers got it right.
• compare the model fit statistics: AIC and BIC are both much better than without fitting the pile-on effect

Let's also add the triadic closure effect before looking at the model in more detail.

m2 <- ergm(nw ~ edges + nodefactor('dem') + edgecov(joint_dem) + gwdegree(decay = 0, fixed = T) + triangle, control = control.ergm(seed = 210616))

Again, seems like we have degeneracy.

• we can use the geometrically weighted edgewise shared partners statistic
• in the ergm package, it is called gwesp

m2 <- ergm(nw ~ edges + nodefactor('dem') + edgecov(joint_dem) + gwdegree(decay = 0, fixed = T) + gwesp(decay = 0, fixed = T), control = control.ergm(seed = 210616))

summary(m2)

Looks like it converged.

• gwesp.fixed.0 is positive here
• this is against what we expected, but it's hard to draw conclusions since the model is quite rudimentary
• model fit is improved again; this is generally what happens with the gwdegree and gwesp terms

3.2.1 Let's consider model fit, degeneracy, and decay parameters

How to understand the decay parameters?

• the lower the decay value, the less likely you are going to run into degeneracy problems
• at the same time, a lower decay value also makes the model fit worse because you are removing information (generally speaking)
• for the geometrically weighted degree, lower decay values makes additional clustered ties matter less (the same interpretation works for additional edgewise shared partners in the gwesp)
• the model is less able to discriminate between the network effect of small and large clusters
• The most useful tool when working with geometrically weighted degrees: https://github.com/michaellevy/gwdegree

Reference:

Levy, M. 2016. "gwdegree: Improving interpretation of geometrically-weighted degree estimates in exponential random graph models." Journal of Open Source Software.

This is easiest to demonstrate using the gwesp.

• let's make a fully saturated direct network with 10 nodes
• this network will have $10\times 9 = 90$ edges
• every edge will have 8 edgewise shared partners

tenclique <- network(matrix(1, ncol = 10, nrow = 10), directed = T)
options(repr.plot.width=7, repr.plot.height=7)
set.seed(210615); plot(tenclique)

The summary function is very useful when provided with an ERGM model.

summary(network ~ ergm-terms)

summary(tenclique ~ edges + gwesp(0, fixed = T) + gwesp(0.5, fixed = T) + gwesp(1, fixed = T) + gwesp(2, fixed = T) + gwesp(5, fixed = T) + gwesp(10, fixed = T) + ttriple)

What happened?

• gwesp with decay set to 0 means every edge has "one or more" shared partner - this makes it the same value as the edge count term
• as we increase the decay value, it starts to converge toward the ttriple term, which is a count of all edgewise shared partners of all edges (i.e. $10\times 9\times 8 = 720$)
• on the one end, you will have an unidentified model (i.e. edges and gwesp are the same) and on the other end, you will run into degeneracy issues

Let's look at a real example using our MIDs network.

• let's increase the decay value to 0.2

m2.b <- ergm(nw ~ edges + nodefactor('dem') + edgecov(joint_dem) + gwdegree(decay = 0.2, fixed = T) + gwesp(decay = 0, fixed = T), control = control.ergm(seed = 210616))

summary(m2.b)

Slightly increasing the decay value improved fit and didn't run into degeneracy issues.

3.3 Assessing ERGM fit

Back to our MIDs model; it seems like our model converged, but we should still check the fit using simulations.

• the ERGM is a generative model, so it can be used to simulate networks that have the same generative features (parameters and configurations)
• simulate a bunch of networks, which lets us check how well our model captures the observed network
• the gof function is what we use

m2fit <- gof(m2, control = control.gof.ergm(seed = 210616))
m2fit

What are we looking at?

• we simulated a bunch of networks, and computed some important network statistics
• the output above are summaries of the simulations compared to the observed value
• it's easier to see when we plot them

options(repr.plot.width=32, repr.plot.height=8)
par(mfrow = c(1, 4))
plot(m2fit)

• black line shows the observed values
• the box and whisker plots show the distributions of the simulated networks
• ideally, we want the observed values to fall in the gray boxes
• we also don't want to "overfit"
• there isn't an exact answer to this task

But it helps to compare.

m1fit <- gof(m1, control = control.gof.ergm(seed = 210616))
options(repr.plot.width=32, repr.plot.height=8)
par(mfrow = c(1, 4))
plot(m1fit)

3.4 Interpreting the ERGM output

Let's bring up the model output again.

summary(m2)

• Coefficients are log odds.

At the network level:

• positive coefficient means positive contribution to the network generative process
• for example, positive gwesp means triadic closure is a feature of the network
• $exp(\theta_r)$ is the odds of observing a network that is one-unit greater in configuration $r$ (while keeping all other configuration counts constant)

At the dyad level:

• how many of each configurations are each edge a part of? (this is the change statistic)
• nodefactor.dem will be 0 if two nondemocracies, 1 if mixed regime, and 2 if joint democracies
• gwesp.fixed.0 will be 1 if the dyad is part of at least one triangle

$P(y_{ij} = 1 | Y, \mathbf{\theta}) = logit ^{-1}(\sum^k_{r=1}\theta_r \delta_r^{(ij)}Y)$

• conditional log odds of a tie given a one unit change in the statistic
• for the less complicated terms, it becomes the conditional log odds of a tie given that it is part of the local configuration one more time

Let's illustrate using a simple example and the same coefficients.

Just a simple three node network with two ties to show tie formation on the empty dyad

tri <- matrix(c(0, 1, 1,
                1, 0, 0,
                1, 0, 0), byrow = T, ncol = 3)
nw.tri <- network(tri, directed = F)
nw.tri%v%'dem' <- c(0, 1, 1)
jd.tri <- matrix(c(0, 0, 0,
                   0, 0, 1,
                   0, 1, 0), byrow = T, ncol = 3)
options(repr.plot.width=10, repr.plot.height=10)
set.seed(210615); plot(nw.tri, vertex.col = ifelse(nw.tri%v%'dem', 2, 1), displaylabels = T)

What's the probability of a tie forming on the empty dyad (2,3) if this network has the same generative features of our conflict network?

• we need to see how the addition of the tie will change the network statistics
• the ergmMPLE function can help with this, especially with more complicated terms like the geometrically weighted ones

ergmMPLE(nw.tri ~ edges + nodefactor('dem') + edgecov(jd.tri) + gwdegree(decay = 0, fixed = T) + gwesp(decay = 0, fixed = T))$predictor

There are two types of dyads here:

• (1,2) and (1,3)
• (2,3)

Adding an edge on the (2,3) dyad will add...

• 1 edge
• 2 democracies
• 1 joint democracy dyad
• 0 to the gwdegree.0 term
• 3 to the gwesp.0 term

Why does gwdegree not change?

Let's look at the degree distribution:

nodes with degree:	0	1	2
w/o (2,3) edge	0	2	1
w/ (2,3) edge	0	0	3

• When the decay value is 0, all node degrees are treated the same
• Adding the (2, 3) edge removes two degree-1 nodes, and adds two degree-2 nodes so gwdegree.0 does not change

round(coefficients(m2), 2)

Calculate the log-odds.

$(-4.70 \times 1) + (-0.15 \times 2) + (-0.21 \times 1) + (-0.64 \times 0) + (1.15 \times 3) = -1.76$

Calculate the corresponding probability.

plogis(-1.76)

• generally speaking, it will get confusing to calculate the probability of a tie forming as the networks get more complex and the model terms more nested

At the block level:

• Subnetwork-based sampling that can be used to simulate networks from different ERGMs
• Allows for comparison on outcomes of interest when varying certain model terms
• Beyond the scope of what we can do in this tutorial

Reference:

Desmarais, BA & Cranmer, SJ. 2012. "Micro-level interpretations of exponential random graph models with applications to estuary networks." Policy Studies Journal.

3.5 ERGM Settings

Let's conclude this section by looking at some settings and useful things to do.

• turn off eval.loglik if you are exploring large and complex models
• set a seed to ensure consistency
• increase MCMC.burnin, MCMC.samplesize, and MCMC.interval for higher quality estimates
• use parallel to speed up estimation when working with larger models

See all the settings with ?control.ergm

m2.c <- ergm(nw ~ edges + nodefactor('dem') + edgecov(joint_dem) + gwdegree(decay = 0, fixed = T) + gwesp(decay = 0, fixed = T),
             eval.loglik = F,
             control = control.ergm(seed = 210615,
                                    MCMC.burnin = 50000,
                                    MCMC.samplesize = 2500,
                                    MCMC.interval = 2500,
                                    parallel = 0))

summary(m2.c)

We don't get any likelihood-based statistics. We can follow the instructions to add them.

m2.c <- logLik(m2.c, add = T)

summary(m2.c)

Increasing the MCMC samples can improve our estimates.

• we didn't talk about estimation, but the primary way of doing estimation for ERGMs is using MCMC-based methods
• these settings increase the MCMC samples and interval

options(repr.plot.width=20, repr.plot.height=10)
mcmc.diagnostics(m2)

options(repr.plot.width=20, repr.plot.height=10)
mcmc.diagnostics(m2.c)

4. Resources and References

Packages

statnet suite of packages have lots of extensions of the ERGM and other network methods

xergm package: lots of extensions to the ERGM

gwdegree package: helps with understanding gwdegree terms (https://github.com/michaellevy/gwdegree)

References

Desmarais, BA & Cranmer, SJ. 2012. "Micro-level interpretations of exponential random graph models with applications to estuary networks." Policy Studies Journal.

Hunter, DR et al. 2008. "ergm: A package to fit, simulate and diagnose exponentiall-family models for networks." Journal of Statistical Software.

Levy, M. 2016. "gwdegree: Improving interpretation of geometrically-weighted degree estimates in exponential random graph models." Journal of Open Source Software.

Snijders, T et al. 2006. "New specifications for exponential random graph models." Sociological Methodology.

5. Questions or comments?

Please visit the Github repository for this notebook.

Or email me at ted.hsuanyun.chen@gmail.com.