Exponential Random Graph Models (ERGMs) using statnet

Prerequisites

This workshop assumes basic familiarity with R; experience with network concepts, terminology and data; and familiarity with the general framework for statistical modeling and inference. While previous experience with ERGMs is not required, some of the topics covered here may be difficult to understand without a strong background in linear and generalized linear models in statistics.

Software installation

Minimally, you will need to install the latest version of R (available here) and the statnet packages ergm and network to run the code presented here (ergm will automatically install network when it is loaded).
The workshops are conducted using the free version of Rstudio (available here).

The full set of installation intructions with details can be found on the statnet workshop wiki.

If you have not already downloaded the statnet packages for this workshop, the quickest way to install these (and the other most commonly used packages from the statnet suite), is to open an R session and type:

install.packages('ergm')

library(ergm)

You can check the version number with:

packageVersion("ergm")

[1] '4.5.0'

Throughout, we will set a random seed via set.seed() for commands in tutorial that require simulating random values—this is not necessary, but it ensures that you will get the same results as this tutorial (assuming that you are using the same ergm version or at least a version in which the algorithms you are using have not changed).

The general form for an ERGM

ERGMs are a class of models, superficially resembling linear regression or GLMs. The general form of the model specifies the probability of the entire network (the left hand side), as a function of terms that represent network features we hypothesize may occur more or less likely than expected by chance (the right hand side). The general form of the model is

where

• $Y$ is the random variable for the state of the network and $y$ is a particular realization $Y$ could take,

• $g(y)$ is a vector of model statistics for network $y$,

• $h(y)$ is the reference measure (which defines the baseline behavior of the model when $\theta=0$)

• $\theta$ is the vector of coefficients for those statistics, and

• $k(\theta)$ is the summation of the numerator’s value over the set of all possible networks $y$, typically taken to be all networks with the same node set as the observed network.

In particular, the model implies that the probability attached to a network $y$ only depends on the network via the vector of statistics $g(y)$. Among other things, this means that maximum likelihood estimation may be carried out even if we don’t observe the network itself, as long as we know the observed value of $g(y)$.

If you’re not familiar with the compact notation above, the numerator represents a formula that is linear in the log form:

$h(y)$ is important for specifying model behavior for valued ERGMs, for controlling for network size in extrapolative or multi-network settings, and other applications. When modeling single, unvalued networks, however, it is common to take $h(y)\propto 1$ (the counting measure). We will adopt that convention here. Under the counting measure, the baseline behavior of an ERGM approaches a uniform random graph when $\theta \to 0$, and the ERGM terms serve to modify this baseline distribution. We will see many examples below.

The model statistics $g(y)$: ERGM terms

The statistics $g(y)$ can be thought of as the “covariates” in the model. In the network modeling context, these represent network features like density, homophily, triads, etc. In one sense, they are like covariates you might use in other statistical models. But they are different in one important respect: these $g(y)$ statistics are functions of the network itself — each is defined by the frequency of a specific configuration of dyads observed in the network — so they are not measured by a question you include in a survey (e.g., the income of a node), but instead need to be computed on the specific network you have, after you have collected the data.

As a result, every term in an ERGM must have an associated algorithm for computing its value for your network. The ergm package in statnet includes about 150 term-computing algorithms. We will explore some of these terms in this tutorial, and links to more information are provided in section 3.

You can get an up-to-date list of all available terms, and the syntax for using them, by typing ?ergmTerm. When using RStudio, it is possible to press the tab key after starting a line with ?ergm to view the wide range of possible help options beginning with the letters ergm.

Available keywords and their meanings can be obtained by typing ?ergmKeyword. You can also search for terms with keywords, as in

search.ergmTerms(keyword='curved')

Found  8  matching ergm terms:
altkstar(lambda, fixed=FALSE) (binary)
    Alternating k-star

gwb1degree(decay, fixed=FALSE, attr=NULL, cutoff=30, levels=NULL) (binary)
    Geometrically weighted degree distribution for the first mode in a bipartite network

gwb1dsp(decay=0, fixed=FALSE, cutoff=30) (binary)
    Geometrically weighted dyadwise shared partner distribution for dyads in the first bipartition

gwb2degree(decay, fixed=FALSE, attr=NULL, cutoff=30, levels=NULL) (binary)
    Geometrically weighted degree distribution for the second mode in a bipartite network

gwb2dsp(decay=0, fixed=FALSE, cutoff=30) (binary)
    Geometrically weighted dyadwise shared partner distribution for dyads in the second bipartition

gwdegree(decay, fixed=FALSE, attr=NULL, cutoff=30, levels=NULL) (binary)
    Geometrically weighted degree distribution

gwidegree(decay, fixed=FALSE, attr=NULL, cutoff=30, levels=NULL) (binary)
    Geometrically weighted in-degree distribution

gwodegree(decay, fixed=FALSE, attr=NULL, cutoff=30, levels=NULL) (binary)
    Geometrically weighted out-degree distribution

To obtain help for a specific term, use either help("[name]-ergmTerm") or the shorthand version ergmTerm?[name], where [name] is the name of the term.

For more guidance on ergm terms, there is a vignette in the ergm package entitled ergm-term-crossRef.Rmd that can be compiled as an RMarkdown document.

One key categorization of model terms is worth keeping in mind: terms are either dyad independent or dyad dependent. Dyad independent terms (like nodal homophily terms) imply no dependence between dyads—the presence or absence of a tie may depend on nodal attributes, but not on the state of other ties. Dyad dependent terms (like degree terms, or triad terms), by contrast, imply dependence between dyads. Dyad dependent terms have very different effects, and much of what is different about network models comes from these terms. They introduce complex cascading effects that can often lead to counter-intuitive and highly non-linear outcomes. In addition, a model with at least one dyad dependent term requires a different estimation algorithm, so when we use these terms below you will see some different components in the output.

ERGM probabilities: at the tie level

The ERGM expression for the probability of the entire graph shown above can be re-expressed in terms of the conditional log-odds (that is, the logit of the conditional probability) of a single tie between two actors:

• $Y_{ij}$ is the random variable for the state of the actor pair $i,j$ (with realization $y_{ij}$), and

• $y^{c}_{ij}$ signifies the complement of $y_{ij}$, i.e. the entire network $y$ except for $y_{ij}$.

• $\delta_{ij}(y)$ is a vector of the “change statistics” for each model term. The change statistic records how the $g(y)$ term changes if the $y_{ij}$ tie is toggled from off to on while fixing the rest of the network. So

  $$\delta_{ij}(y) = g(y^{+}_{ij})-g(y^{-}_{ij}),$$ where

• $y^{+}_{ij}$ is defined as $y^{c}_{ij}$ along with $y_{ij}$ set to 1, and

• $y^{-}_{ij}$ is defined as $y^{c}_{ij}$ along with $y_{ij}$ set to 0.

So $\delta_{ij}(y)$ equals the value of $g(y)$ when $y_{ij}=1$ minus the value of $g(y)$ when $y_{ij}=0$, but all other dyads are as in $y$. When this vector of change statistics is multiplied by the vector of coefficients $\theta$, the equation above shows that this dot product is the log-odds of the tie between $i$ and $j$, conditional on all other dyads remaining the same.

In other words, for an individual statistic, its change value for $Y_{ij}$ times its corresponding coefficient can be interpreted as that term’s contribution to the log-odds of that tie, conditional on all other dyads remaining the same.

We will see exactly how this works in the sections that follow.

Loading network data

Network data can come in many different forms — ties can be stored as edgelists or sociomatrices in .csv files, or as exported data from other programs like Pajek. Attributes for the nodes, ties, and dyads can also come in various forms. All can be read into R using either standard R tools (e.g., for .csv files), or methods from the network package. For more information, refer to the following:

?read.paj
?read.paj.simplify
?loading.attributes

However you read them in, the data will need to be transformed into a network object, the format that Statnet packages use to store and work with network data. For information on how to do this, refer to:

?network

The ergm package also contains several network data sets, and we will use those here for demonstration purposes.

data(package='ergm') # tells us the datasets in our packages

We’ll start with Padgett’s data on Renaissance Florentine families for our first example. As with all data analysis, it is good practice to start by summarizing our data using graphical and numerical descriptives.

set.seed(123) # The plot.network function uses random values
data(florentine) # loads flomarriage and flobusiness data
flomarriage # Equivalent to print.network(flomarriage): Examine properties

 Network attributes:
  vertices = 16 
  directed = FALSE 
  hyper = FALSE 
  loops = FALSE 
  multiple = FALSE 
  bipartite = FALSE 
  total edges= 20 
    missing edges= 0 
    non-missing edges= 20 

 Vertex attribute names: 
    priorates totalties vertex.names wealth 

No edge attributes

par(mfrow=c(1,2)) # Set up a 2-column (and 1-row) plot area
plot(flomarriage, 
     main="Florentine Marriage", 
     cex.main=0.8, 
     label = network.vertex.names(flomarriage)) # Equivalent to plot.network(...)
wealth <- flomarriage %v% 'wealth' # %v% references vertex attributes
wealth

 [1]  10  36  55  44  20  32   8  42 103  48  49   3  27  10 146  48

plot(flomarriage, 
     vertex.cex=wealth/25, # Make vertex size proportional to wealth attribute
     main="Florentine marriage by wealth", cex.main=0.8) 

The summary and ergm functions, and supporting functions

We’ll start by running some simple models to demonstrate the most commonly used functions for ERG modeling.

The syntax for specifying a model in the ergm package follows R’s formula convention:

This syntax is used for both the summary and ergm functions. The summary function simply returns the numerical values of the network statistics in the model. The ergm function estimates the model with those statistics.

It is good practice to run a summmary command on any model before fitting it with ergm. This is the ERGM equivalent of performing some descriptive analysis on your covariates. This can help you make sure you understand what the term represents, and it can help to flag potential problems that will lead to poor modeling results. We will now demonstrate the summary and ergm commands using a simple model.

summary(flomarriage ~ edges) # Calculate the edges statistic for this network

edges 
   20 

flomodel.01 <- ergm(flomarriage ~ edges) # Estimate the model 

Starting maximum pseudolikelihood estimation (MPLE):

Obtaining the responsible dyads.

Evaluating the predictor and response matrix.

Maximizing the pseudolikelihood.

Finished MPLE.

Evaluating log-likelihood at the estimate. 

summary(flomodel.01) # Look at the fitted model object

Call:
ergm(formula = flomarriage ~ edges)

Maximum Likelihood Results:

      Estimate Std. Error MCMC % z value Pr(>|z|)    
edges  -1.6094     0.2449      0  -6.571   <1e-04 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

     Null Deviance: 166.4  on 120  degrees of freedom
 Residual Deviance: 108.1  on 119  degrees of freedom
 
AIC: 110.1  BIC: 112.9  (Smaller is better. MC Std. Err. = 0)

set.seed(321)
summary(flomarriage~edges+triangle) # Look at the g(y) statistics for this model

   edges triangle 
      20        3 

flomodel.02 <- ergm(flomarriage~edges+triangle) # Estimate the theta coefficients
summary(flomodel.02)

Call:
ergm(formula = flomarriage ~ edges + triangle)

Monte Carlo Maximum Likelihood Results:

         Estimate Std. Error MCMC % z value Pr(>|z|)    
edges     -1.6900     0.3620      0  -4.668   <1e-04 ***
triangle   0.1901     0.5982      0   0.318    0.751    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

     Null Deviance: 166.4  on 120  degrees of freedom
 Residual Deviance: 108.1  on 118  degrees of freedom
 
AIC: 112.1  BIC: 117.6  (Smaller is better. MC Std. Err. = 0.01102)

plogis(coef(flomodel.02)[[1]] + (0:2) * coef(flomodel.02)[[2]])

[1] 0.1557799 0.1824455 0.2125265

class(flomodel.02) # this has the class ergm

[1] "ergm"

names(flomodel.02) # the ERGM object contains lots of components.

 [1] "coefficients"    "sample"          "iterations"      "MCMCtheta"      
 [5] "loglikelihood"   "gradient"        "hessian"         "covar"          
 [9] "failure"         "newnetwork"      "coef.init"       "est.cov"        
[13] "coef.hist"       "stats.hist"      "steplen.hist"    "control"        
[17] "etamap"          "MCMCflag"        "nw.stats"        "call"           
[21] "network"         "ergm_version"    "info"            "MPLE_is_MLE"    
[25] "drop"            "offset"          "estimable"       "formula"        
[29] "reference"       "constraints"     "obs.constraints" "estimate"       
[33] "estimate.desc"   "null.lik"        "mle.lik"        

coef(flomodel.02) # you can extract/inspect individual components

    edges  triangle 
-1.689969  0.190103 

summary(wealth) # summarize the distribution of wealth

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   3.00   17.50   39.00   42.56   48.25  146.00 

# plot(flomarriage, 
#      vertex.cex=wealth/25, 
#      main="Florentine marriage by wealth", 
#      cex.main=0.8) # network plot with vertex size proportional to wealth
summary(flomarriage~edges+nodecov('wealth')) # observed statistics for the model

         edges nodecov.wealth 
            20           2168 

flomodel.03 <- ergm(flomarriage~edges+nodecov('wealth'))

Starting maximum pseudolikelihood estimation (MPLE):

Obtaining the responsible dyads.

Evaluating the predictor and response matrix.

Maximizing the pseudolikelihood.

Finished MPLE.

Evaluating log-likelihood at the estimate. 

summary(flomodel.03)

Call:
ergm(formula = flomarriage ~ edges + nodecov("wealth"))

Maximum Likelihood Results:

                Estimate Std. Error MCMC % z value Pr(>|z|)    
edges          -2.594929   0.536056      0  -4.841   <1e-04 ***
nodecov.wealth  0.010546   0.004674      0   2.256   0.0241 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

     Null Deviance: 166.4  on 120  degrees of freedom
 Residual Deviance: 103.1  on 118  degrees of freedom
 
AIC: 107.1  BIC: 112.7  (Smaller is better. MC Std. Err. = 0)

data(faux.mesa.high) 
mesa <- faux.mesa.high

set.seed(1)
mesa

 Network attributes:
  vertices = 205 
  directed = FALSE 
  hyper = FALSE 
  loops = FALSE 
  multiple = FALSE 
  bipartite = FALSE 
  total edges= 203 
    missing edges= 0 
    non-missing edges= 203 

 Vertex attribute names: 
    Grade Race Sex 

No edge attributes

par(mfrow=c(1,1)) # Back to 1-panel plots
plot(mesa, vertex.col='Grade')
legend('bottomleft',fill=7:12,
       legend=paste('Grade',7:12),cex=0.75)

fauxmodel.01 <- ergm(mesa ~edges + 
        nodefactor('Grade') + nodematch('Grade',diff=T) +
        nodefactor('Race') + nodematch('Race',diff=T))

Observed statistic(s) nodematch.Race.Black and nodematch.Race.Other are at their smallest attainable values. Their coefficients will be fixed at -Inf.

Starting maximum pseudolikelihood estimation (MPLE):

Obtaining the responsible dyads.

Evaluating the predictor and response matrix.

Maximizing the pseudolikelihood.

Finished MPLE.

Evaluating log-likelihood at the estimate. 

summary(fauxmodel.01)

Call:
ergm(formula = mesa ~ edges + nodefactor("Grade") + nodematch("Grade", 
    diff = T) + nodefactor("Race") + nodematch("Race", diff = T))

Maximum Likelihood Results:

                      Estimate Std. Error MCMC % z value Pr(>|z|)    
edges                  -8.0538     1.2561      0  -6.412  < 1e-04 ***
nodefactor.Grade.8      1.5201     0.6858      0   2.216 0.026663 *  
nodefactor.Grade.9      2.5284     0.6493      0   3.894  < 1e-04 ***
nodefactor.Grade.10     2.8652     0.6512      0   4.400  < 1e-04 ***
nodefactor.Grade.11     2.6291     0.6563      0   4.006  < 1e-04 ***
nodefactor.Grade.12     3.4629     0.6566      0   5.274  < 1e-04 ***
nodematch.Grade.7       7.4662     1.1730      0   6.365  < 1e-04 ***
nodematch.Grade.8       4.2882     0.7150      0   5.997  < 1e-04 ***
nodematch.Grade.9       2.0371     0.5538      0   3.678 0.000235 ***
nodematch.Grade.10      1.2489     0.6233      0   2.004 0.045111 *  
nodematch.Grade.11      2.4521     0.6124      0   4.004  < 1e-04 ***
nodematch.Grade.12      1.2987     0.6981      0   1.860 0.062824 .  
nodefactor.Race.Hisp   -1.6659     0.2963      0  -5.622  < 1e-04 ***
nodefactor.Race.NatAm  -1.4725     0.2869      0  -5.132  < 1e-04 ***
nodefactor.Race.Other  -2.9618     1.0372      0  -2.856 0.004296 ** 
nodefactor.Race.White  -0.8488     0.2958      0  -2.869 0.004112 ** 
nodematch.Race.Black      -Inf     0.0000      0    -Inf  < 1e-04 ***
nodematch.Race.Hisp     0.6912     0.3451      0   2.003 0.045153 *  
nodematch.Race.NatAm    1.2482     0.3550      0   3.517 0.000437 ***
nodematch.Race.Other      -Inf     0.0000      0    -Inf  < 1e-04 ***
nodematch.Race.White    0.3140     0.6405      0   0.490 0.623947    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

     Null Deviance: 28958  on 20889  degrees of freedom
 Residual Deviance:  1798  on 20868  degrees of freedom
 
AIC: 1836  BIC: 1987  (Smaller is better. MC Std. Err. = 0)

 Warning: The following terms have infinite coefficient estimates:
  nodematch.Race.Black nodematch.Race.Other 

table(mesa %v% 'Race') # Frequencies of race

Black  Hisp NatAm Other White 
    6   109    68     4    18 

mixingmatrix(mesa, "Race")

      Black Hisp NatAm Other White
Black     0    8    13     0     5
Hisp      8   53    41     1    22
NatAm    13   41    46     0    10
Other     0    1     0     0     0
White     5   22    10     0     4

Note:  Marginal totals can be misleading for undirected mixing matrices.

summary(mesa ~edges  + 
          nodefactor('Grade') + nodematch('Grade',diff=T) +
          nodefactor('Race') + nodematch('Race',diff=T))

                edges    nodefactor.Grade.8    nodefactor.Grade.9 
                  203                    75                    65 
  nodefactor.Grade.10   nodefactor.Grade.11   nodefactor.Grade.12 
                   36                    49                    28 
    nodematch.Grade.7     nodematch.Grade.8     nodematch.Grade.9 
                   75                    33                    23 
   nodematch.Grade.10    nodematch.Grade.11    nodematch.Grade.12 
                    9                    17                     6 
 nodefactor.Race.Hisp nodefactor.Race.NatAm nodefactor.Race.Other 
                  178                   156                     1 
nodefactor.Race.White  nodematch.Race.Black   nodematch.Race.Hisp 
                   45                     0                    53 
 nodematch.Race.NatAm  nodematch.Race.Other  nodematch.Race.White 
                   46                     0                     4 

set.seed(2)
data(samplk) # directed data: Sampson's Monks
ls() 

 [1] "faux.magnolia.high" "faux.mesa.high"     "fauxmodel.01"      
 [4] "fit"                "flobusiness"        "flomarriage"       
 [7] "flomodel.01"        "flomodel.02"        "flomodel.03"       
[10] "flomodel.03.gof"    "flomodel.03.sim"    "magnolia"          
[13] "mesa"               "mesamodel.02"       "mesamodel.02.gof"  
[16] "missnet"            "missnet_bad"        "missnetmat"        
[19] "samplk1"            "samplk2"            "samplk3"           
[22] "sampmodel.01"       "tempnet"            "wealth"            

samplk3

 Network attributes:
  vertices = 18 
  directed = TRUE 
  hyper = FALSE 
  loops = FALSE 
  multiple = FALSE 
  bipartite = FALSE 
  total edges= 56 
    missing edges= 0 
    non-missing edges= 56 

 Vertex attribute names: 
    cloisterville group vertex.names 

No edge attributes

plot(samplk3)

summary(samplk3~edges+mutual)

 edges mutual 
    56     15 

set.seed(3)
sampmodel.01 <- ergm(samplk3~edges+mutual)
summary(sampmodel.01)

Call:
ergm(formula = samplk3 ~ edges + mutual)

Monte Carlo Maximum Likelihood Results:

       Estimate Std. Error MCMC % z value Pr(>|z|)    
edges   -2.1639     0.2211      0  -9.789   <1e-04 ***
mutual   2.3118     0.4860      0   4.757   <1e-04 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

     Null Deviance: 424.2  on 306  degrees of freedom
 Residual Deviance: 268.1  on 304  degrees of freedom
 
AIC: 272.1  BIC: 279.6  (Smaller is better. MC Std. Err. = 0.3199)

Terms provided with ergm

The ergm package provides myriad terms, and it can be difficult to absorb the full array of available model terms in any one place. This is particularly true with the release of ergm version 4.0, which expands the user’s ability to create terms even further, for example through the use of term operators. As mentioned above in Section 1, it is possible to search for specific topics using search.ergmTerms; to obtain help on a particular term called [name] using ergmTerm?[name], where [name] is the name of the term; or to see the full list of available terms using ?ergmTerm.

The list of all terms is quite lengthy, so it may be helpful to start with a more concise list such as the one found here. A more detailed discussion can be found in volume 24, issue 4 of the Journal of Statistical Software.

To appreciate the expanded capabilities of the ergm package as of the release of version 4.0, we recommend Krivitsky et al (2021). In this article, Section 3 describes the enhanced flexibility to create specialized model terms involving functions of nodal covariates, and Section 4 explains how operators further extend the types of terms at the user’s disposal.

Coding new ergm-terms

There is a statnet package called ergm.userterms that provides the utilities needed to write new ergm-terms. The package is available via GitHub at https://github.com/statnet/ergm.userterms, and installing it will include the tutorial, called ergmuserterms.pdf. A tutorial can also be found in the Journal of Statistical Software 52(2), and some introductory slides and installation instructions from the workshop we teach on coding ergm-terms can be found on GitHub.

Writing up new ergm terms requires some knowledge of C and the ability to build R from source.

What it looks like when a model converges properly

We will first consider a simple dyadic dependent model where the algorithm works using the program defaults, with a degree(1) term that captures whether there are more (or less) degree 1 nodes than we would expect, given the density.

set.seed(314159)
summary(flobusiness~edges+degree(1))

  edges degree1 
     15       3 

fit <- ergm(flobusiness~edges+degree(1))
summary(fit)

Call:
ergm(formula = flobusiness ~ edges + degree(1))

Monte Carlo Maximum Likelihood Results:

        Estimate Std. Error MCMC % z value Pr(>|z|)    
edges    -2.1177     0.3032      0  -6.984   <1e-04 ***
degree1  -0.6272     0.6010      0  -1.044    0.297    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

     Null Deviance: 166.36  on 120  degrees of freedom
 Residual Deviance:  89.39  on 118  degrees of freedom
 
AIC: 93.39  BIC: 98.96  (Smaller is better. MC Std. Err. = 0.03364)

mcmc.diagnostics(fit)

Sample statistics summary:

Iterations = 7168:131072
Thinning interval = 512 
Number of chains = 1 
Sample size per chain = 243 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

           Mean    SD Naive SE Time-series SE
edges    0.2140 3.601   0.2310         0.2310
degree1 -0.1276 1.817   0.1166         0.1166

2. Quantiles for each variable:

        2.5% 25% 50% 75% 97.5%
edges     -8  -2   1   2     7
degree1   -3  -1   0   1     4


Are sample statistics significantly different from observed?
               edges    degree1    (Omni)
diff.      0.2139918 -0.1275720        NA
test stat. 0.9262347 -1.0944112 1.4827084
P-val.     0.3543240  0.2737747 0.4790184

Sample statistics cross-correlations:
             edges    degree1
edges    1.0000000 -0.3929828
degree1 -0.3929828  1.0000000

Sample statistics auto-correlation:
Chain 1 
                edges     degree1
Lag 0     1.000000000  1.00000000
Lag 512  -0.018373256  0.06527360
Lag 1024 -0.008853872 -0.00419178
Lag 1536 -0.006593784 -0.05395258
Lag 2048  0.033260731  0.02580333
Lag 2560 -0.059894956  0.02109630

Sample statistics burn-in diagnostic (Geweke):
Chain 1 

Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5 

    edges   degree1 
1.3768812 0.1355053 

Individual P-values (lower = worse):
    edges   degree1 
0.1685490 0.8922124 
Joint P-value (lower = worse):  0.1257365 

Note: MCMC diagnostics shown here are from the last round of
  simulation, prior to computation of final parameter estimates.
  Because the final estimates are refinements of those used for this
  simulation run, these diagnostics may understate model performance.
  To directly assess the performance of the final model on in-model
  statistics, please use the GOF command: gof(ergmFitObject,
  GOF=~model).

What this shows is a summary of the statistics generated by the MCMC process, with each row summarizing a different statistic and with each statistic measured in terms of its value relative to the corresponding value for the original observed network. On the left is a “traceplot” in which the values are plotted as a function of iteration number; while on the right is a histogram-like plot of the whole sample of statistics without regard to their order in the MCMC process.

This example exhibits the sort of behavior that we want to see in the MCMC diagnostics: The MCMC sample statistics are varying randomly around the observed values at each step; we might say that the chain is “mixing well”. The sampled values show little serial correlation, indicating that they are independent draws, and they have a roughly bell-shaped distribution, centered at zero. The sawtooth pattern visible in the degree term deviation plot is due to the combination of discrete values and small range in the statistics: the observed number of degree 1 nodes is 3, and only a few discrete values are produced by the simulations. So the sawtooth pattern is is an inherent property of the statistic, not a problem with the model.

There are many control parameters for the MCMC algorithm (help(snctrl)), and we’ll play with some of these below. To see what the algorithm is doing at each step, we can drop the sampling interval down to 1:

set.seed(271828)
fit <- ergm(flobusiness~edges+degree(1),
            control=snctrl(MCMC.interval=1))

This runs an MCMC algorithm where every network’s statistics are returned, which might be useful if we are trying to debug a bad model fit.

In the last section we’ll look at some models that don’t converge properly, and how to use MCMC diagnostics to identify and address this.

What it looks like when a model fails

For this purpose, we’ll use a larger network, faux.magnolia.high, and look at a simple model for triad closure that includes only edges and triangle terms.

set.seed(10)
data('faux.magnolia.high')
magnolia <- faux.magnolia.high
magnolia

 Network attributes:
  vertices = 1461 
  directed = FALSE 
  hyper = FALSE 
  loops = FALSE 
  multiple = FALSE 
  bipartite = FALSE 
  total edges= 974 
    missing edges= 0 
    non-missing edges= 974 

 Vertex attribute names: 
    Grade Race Sex vertex.names 

 Edge attribute names not shown 

plot(magnolia, vertex.cex=.5)

summary(magnolia~edges+triangle) # Simple model for triad closure

   edges triangle 
     974      169 

We now try to fit this “simple” model:

set.seed(100)
fit <- ergm(magnolia~edges+triangle,
            control=snctrl(MCMLE.effectiveSize=NULL))

Starting maximum pseudolikelihood estimation (MPLE):
Evaluating the predictor and response matrix.
Maximizing the pseudolikelihood.
Finished MPLE.
Starting Monte Carlo maximum likelihood estimation (MCMLE):
...
Iteration 4 of at most 60:
Optimizing with step length 0.3963.
The log-likelihood improved by 1.1568.
Estimating equations are not within tolerance region.
Iteration 5 of at most 60:
Error in ergm.MCMLE(init, nw, model, initialfit = (initialfit <- NULL),  : 
  Number of edges in a simulated network exceeds that in the observed by a factor of more than 20. This is a strong indicator of model degeneracy or a very poor starting parameter configuration. If you are reasonably certain that neither of these is the case, increase the MCMLE.density.guard control.ergm() parameter.

Very interesting. Instead of converging, the algorithm heads off into networks that are much much more dense than the observed network. This is such a clear indicator of a degenerate model specification that the algorithm stops after 3 iterations, to avoid storage problems. To peek a bit more under the hood, we can stop the algorithm earlier, by setting MCMLE.maxit=2, to catch where it’s heading:

set.seed(1000)
fit <- ergm(magnolia~edges+triangle, 
            control=snctrl(MCMLE.maxit=2,MCMLE.effectiveSize=NULL))

Starting maximum pseudolikelihood estimation (MPLE):
Evaluating the predictor and response matrix.
Maximizing the pseudolikelihood.
Finished MPLE.
Starting Monte Carlo maximum likelihood estimation (MCMLE):
Iteration 1 of at most 2:
Optimizing with step length 0.2805.
The log-likelihood improved by 3.0798.
Estimating equations are not within tolerance region.
Iteration 2 of at most 2:
Optimizing with step length 0.0420.
The log-likelihood improved by 4.6627.
Estimating equations are not within tolerance region.
MCMLE estimation did not converge after 2 iterations. The estimated coefficients may not be accurate. Estimation may be resumed by passing the coefficients as initial values; see 'init' under ?control.ergm for details.
Finished MCMLE.
Evaluating log-likelihood at the estimate. Fitting the dyad-independent submodel...
Bridging between the dyad-independent submodel and the full model...
Setting up bridge sampling...
Using 16 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 .
Bridging finished.
This model was fit using MCMC.  To examine model diagnostics and check for degeneracy, use the mcmc.diagnostics() function.

Let’s use the MCMC diagnostics from Section 4 to get a sense of what happened:

mcmc.diagnostics(fit)

For the diagnostic plots, the simulated network statistics are subtracted from their observed values so that the observed values equal zero. Clearly, this Markov chain is heading somewhere very bad!

The edges + triangle model class turns out to be one of the classic degenerate model specifications, and we now understand much more about why it does not produce reasonable levels of triadic closure.

We also now have a more robust way of modeling triangles: the geometrically-weighed edgewise shared partner term (GWESP). For a technical introduction to GWESP, see Hunter and Handcock, 2006; for a more intuitive description and empirical application, see Goodreau, Kitts & Morris, 2009 )

Let’s see what using gwesp instead of triangle can do. We can also control the number of Metropolis-Hastings (MCMC) proposals between sampled statistics in our Markov chain, one of the many control parameters that may be passed to functions in the ergm package using the control=snctrl() syntax. (To see the many control parameters that may be set by the user in the ergm package, type ?snctrl.)

set.seed(10101)
fit <- ergm(magnolia~edges+gwesp(0.25, fixed=T), 
            control=snctrl(MCMC.interval = 10000),
            verbose=T)

Evaluating network in model.
Initializing unconstrained Metropolis-Hastings proposal: ‘ergm:MH_TNT’.
Initializing model...
Model initialized.
Using initial method 'MPLE'.
Fitting initial model.
Starting maximum pseudolikelihood estimation (MPLE):
Evaluating the predictor and response matrix.
Maximizing the pseudolikelihood.
Finished MPLE.
Starting Monte Carlo maximum likelihood estimation (MCMLE):

 ... (output snipped)

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

mcmc.diagnostics(fit)

Sample statistics summary:

Iterations = 2800000:55400000
Thinning interval = 40000 
Number of chains = 1 
Sample size per chain = 1316 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                  Mean    SD Naive SE Time-series SE
edges            7.866 39.98   1.1020          4.141
gwesp.fixed.0.25 7.206 31.99   0.8819          3.360

2. Quantiles for each variable:

                   2.5%    25%   50%   75% 97.5%
edges            -66.12 -18.25 6.000 34.00 93.12
gwesp.fixed.0.25 -57.14 -13.50 8.166 27.07 71.32


Are sample statistics significantly different from observed?
                edges gwesp.fixed.0.25     (Omni)
diff.      7.86626140       7.20579674         NA
test stat. 1.89941461       2.14437967 6.19661638
P-val.     0.05750998       0.03200248 0.04675916

Sample statistics cross-correlations:
                     edges gwesp.fixed.0.25
edges            1.0000000        0.7833691
gwesp.fixed.0.25 0.7833691        1.0000000

Sample statistics auto-correlation:
Chain 1 
               edges gwesp.fixed.0.25
Lag 0      1.0000000        1.0000000
Lag 40000  0.5460541        0.8587880
Lag 80000  0.4618254        0.7501801
Lag 120000 0.4129546        0.6642087
Lag 160000 0.3832516        0.5940199
Lag 2e+05  0.3082655        0.5300815

Sample statistics burn-in diagnostic (Geweke):
Chain 1 

Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5 

           edges gwesp.fixed.0.25 
       0.1445195       -0.1056854 

Individual P-values (lower = worse):
           edges gwesp.fixed.0.25 
       0.8850903        0.9158320 
Joint P-value (lower = worse):  0.3733633 

Note: MCMC diagnostics shown here are from the last round of
  simulation, prior to computation of final parameter estimates.
  Because the final estimates are refinements of those used for this
  simulation run, these diagnostics may understate model performance.
  To directly assess the performance of the final model on in-model
  statistics, please use the GOF command: gof(ergmFitObject,
  GOF=~model).

MORAL: Unwanted degeneracy is an indicator of a poorly specified model. It is not a property of all ERGMs (nor is it unique to ERGMs!), but it is associated with some dyadic-dependent terms, in particular, the reduced homogeneous Markov specifications (e.g., 2-stars and triangle terms). For a good technical discussion of unstable terms, see Schweinberger 2012. For a discussion of alternative terms that exhibit more stable behavior, see Snijders et al. 2006.. For the gwesp term and the curved exponential family terms in general, see Hunter and Handcock 2006.. Note that there are cases in which degenerate behavior may be realistic and desired (e.g., in modeling small groups that do genuinely collapse into cliques or other regular structures, or in some physical applications), but these are not frequently encountered in modeling medium-to-large social networks. If your system is not degenerate but your model is, this suggests that you may need to rethink the theoretical motivation behind your choice of model terms.

Current statnet packages

Packages developed by the Statnet team that are not covered in this tutorial:

• sna — classical social network analysis utilities
• tsna — descriptive statistics for temporal network data
• tergm — temporal ergms for dynamic networks
• ergm.ego— estimation/simulation of ergms from egocentrically sampled data
• ergm.count — models for tie count network data
• ergm.rank — models for tie rank network data
• ergmgp — continuous time network dynamics with ERGM equilibria
• relevent — relational event models for networks
• latentnet — latent space and latent cluster analysis
• degreenet — MLE estimation for degree distributions (negative binomial, Poisson, scale-free, etc.)
• networksis — simulation of bipartite networks with given degree distributions
• ndtv package — network movie maker
• EpiModel — network modeling of infectious disease and social diffusion processes
• ergm.multi — ERGMs for multiple or multilayer networks
• ergm.userterms — template for users who want to implement their own new ERGM terms. (available on GitHub only)

Many of these packages have associated training workshops. Our tutorials can be found online, on the GitHub statnet Workshops wiki.

Additional functionality in base ergm

The ergm package has considerable functionality beyond what has been discussed here. This includes support for:

• ERGMs for valued ties
• Constraints, and multi-network estimation
• Offsets and size correction
• Fine-tuning of simulation and estimation
• Changescore calculation and other utilities

Some of these topics are covered in our advanced ergm workshops.

Extensions by other developers

There are now a number of excellent packages developed by others that extend the functionality of statnet. The easiest way to find these is to look at the “reverse depends” of the ergm package on CRAN. Examples include:

• Bergm — Bayesian Exponential Random Graph Models
• btergm — Temporal Exponential Random Graph Models by Bootstrapped Pseudolikelihood
• hergm — hierarchical ERGMs for multi-level network data
• xergm — extensions to ERGM modeling

Statnet Commons: The development group

Mark S. Handcock <handcock@stat.ucla.edu>

David R. Hunter <dhunter@stat.psu.edu>

Carter T. Butts <buttsc@uci.edu>

Steven M. Goodreau <goodreau@u.washington.edu>

Skye Bender-deMoll <skyebend@skyeome.net>

Martina Morris <morrism@u.washington.edu>

Pavel N. Krivitsky <p.krivitsky@unsw.edu.au>

Samuel M. Jenness <samuel.m.jenness@emory.edu>

Chad Klumb <cklumb@gmail.com>

Michał Bojanowski <mbojanowski@kozminski.edu.pl>