Day 11: Advanced ERGM Specification and Temporal ERGMs

ICPSR Network Analysis: Advanced Topics · Day 11

Author

Shahryar Minhas

Published

July 30, 2026

How to use this document

This is a walkthrough, not slides. Every number, table, and figure is produced by a chunk in this file and is already rendered in the HTML version. During class, use the rendered walkthrough and teaching deck. The .qmd remains available so you can reproduce the analysis later without relying on screenshots or copied console output.

The rendered checkpoints show the short lines of output you should obtain if you rerun the analysis later.
Code wrapped in fit_or_load("name", { ... }) either returns the shipped cached result or fits the model when a cached result is unavailable. Read the code inside the braces as the reproducible specification, not as something that must be estimated during class.
Folded depth boxes are optional. Folded stuck? boxes are hints. Open whichever one you need and ignore the other. I wrote both; nobody is expected to read both.

Before You Reproduce the Walkthrough

Setup check for running the code yourself

If you plan to rerun the walkthrough, begin with the block below. SETUP OK confirms that the required packages and example data are available. The class presentation uses the rendered results, so this check is for later reproduction rather than estimation during class.

require_pkgs(c("ergm", "network", "sna", "ergm.count", "igraph", "netify"))
library(ergm)
data(faux.mesa.high); data(sampson); data(samplk)
# Installed-but-not-loaded check for the temporal section (C4); packageVersion()
# reads metadata only, so neither package's namespace loads here. Deliberate --
# see the load-order note in C4.
stopifnot(network::network.size(faux.mesa.high) == 205,
          network::network.size(samplike) == 18,
          nzchar(system.file(package = "btergm")),
          nzchar(system.file(package = "RSiena")))
cat("SETUP OK ·  R", as.character(getRversion()),
    "· ergm", as.character(packageVersion("ergm")),
    "· ergm.count", as.character(packageVersion("ergm.count")),
    "· btergm", as.character(packageVersion("btergm")),
    "· RSiena", as.character(packageVersion("RSiena")),
    "\n")
#> SETUP OK ·  R 4.3.3 · ergm 4.12.0 · ergm.count 4.1.3 · btergm 1.11.1 · RSiena 1.5.0

That is the whole dependency list. No data files, no install_github, nothing off CRAN. The school and monastery examples ship inside ergm, and the panel data for the final section ship inside RSiena. Keeping the materials self-contained makes the later reproduction less fragile.

Already comfortable with some of this? Where to jump

Already fluent in the exponential-family form and change statistics? Start at block B.
Fit ERGMs before, and here for degeneracy, diagnostics and goodness of fit? Start at block C.
Have panel data and here for the temporal extension? C4 needs C1’s fixed = TRUE warning and nothing else.
Want the mathematics the core path deliberately walks past? Appendix A.
Here mainly to read ERGM papers rather than fit them? The checklist in the wrap is the thing to take away, and appendix B has the referee exercise.

A note on the reading list

The Krivitsky (2012) reading on valued ERGMs belongs with the valued-network material in section C3. Nothing on the SRM half of Day 9 requires it.

Also: Hunter, Handcock, Butts, Goodreau & Morris (2008) remains one of the best descriptions of the goodness-of-fit logic in print. But it describes ergm version 2. We are on 4.12.0. Several things it shows you will not run today, and a couple will run but mean something different. I list exactly what changed in Appendix A. Read it alongside, not instead.

Start Here: What We Are Refreshing and Building On

Olga’s ERGM sessions already did the foundational work. You have seen why a random-graph baseline matters, why an ordinary dyadic logit is not a model of the joint graph, how an ERGM assigns probability to networks, what a change statistic is, why the normalizing constant requires simulation, and how edges, nodematch, mutual, triangle, and gwesp enter a specification. You have also fitted models, inspected MCMC output, used gof(), and worked through the Mesa and rebel-cooperation applications.

We will not repeat Olga’s full introduction, but we will refresh the pieces needed to make today’s argument. That matters in a compressed course. Recognizing gwesp in a formula is not the same as remembering what question the model answers, why its estimation requires simulation, or how it differs from the latent models we just fitted. The detailed derivations remain below as a reference, but they are not part of the central lecture route.

Already covered	What we do with it today
Change statistics and conditional log-odds	Translate fitted terms into comparisons a substantive reader can evaluate
MCMLE, MCMC diagnostics, and GOF	Separate estimation failure, degeneracy, and substantive misspecification
Triangles and GWESP	Ask whether apparent closure survives degree heterogeneity and rival representations
ERGM specification and applications	Audit claims that change across nested specifications
Cross-sectional ERGMs	Move to conditional models for network transitions

A Fast Refresher: Question, Contrast, Estimation, and Interpretation

Here are the four things to recover before we move forward.

1. What Is an ERGM Meant to Answer?

An ERGM asks which kinds of networks are more or less plausible under a declared set of actor, dyadic, and structural statistics. It does not merely ask whether the observed graph contains triangles or reciprocal ties. Nearly every observed network does. It asks whether networks with more of a feature receive greater probability after the other statistics in the model are held fixed.

The whole-network probability model is

\[ \Pr_{\theta}(Y=y)= \frac{\exp\{\theta^{\mathsf T}g(y)\}} {\sum_{y'\in\mathcal Y}\exp\{\theta^{\mathsf T}g(y')\}}. \]

For interpretation, we usually turn that joint distribution into a local comparison: how do the conditional log-odds change when we toggle one eligible tie and leave the rest of the network fixed? If adding $i \to j$ creates a mutual dyad, closes two-paths, or adds a tie between similar actors, the corresponding change statistics tell us how the model scores that toggle.

So the answer is not simply “closure matters.” In the Mesa friendship example below, the useful question is whether a friendship that would connect two students through shared friends receives greater conditional probability after grade, sex, density, and uneven numbers of friends are represented. That sentence identifies the actors, the action, the rival explanations, and the comparison the coefficient actually makes.

2. How Is This Different From What We Just Did?

The recent models all take dependence seriously, but they place it in different parts of the model and therefore answer different questions.

Model	What it primarily asks	Where network dependence goes
Dyadic regression	How are measured covariates associated with a tie?	Usually nowhere unless dependence is modeled separately
SRM or AME	How are covariates associated with ties after accounting for sender, receiver, reciprocity, and latent relational structure?	Additive and multiplicative latent effects in the conditional mean, plus the model’s dyadic error structure
Latent distance model	Can unobserved social positions explain which actors are likely to connect?	Distance between estimated latent positions
ERGM	Which actor, dyadic, and graph configurations make the whole observed network more plausible?	Explicit statistics of the observed graph distribution

The contrast with Day 10 is especially important. AME can represent clustering through latent factors, re-estimate covariate associations conditional on that structure, and sometimes improve prediction without naming the clustering process. An ERGM can put a coefficient directly on a shared-partner configuration. That makes closure available as an estimand, but it does not prove that actors intentionally seek closure. Unobserved homophily can still produce the same pattern.

3. How Is an ERGM Estimated?

For a realistic graph, the denominator above sums over far too many possible networks to calculate directly. MCMC maximum likelihood estimation replaces that impossible enumeration with a simulation loop:

Start with a provisional coefficient vector $\theta$.
Use MCMC tie toggles to simulate whole networks under that provisional model.
Calculate $g(Y)$ for every simulated network.
Compare the simulated average statistics with $g(y_{obs})$ from the observed network.
Update $\theta$ so the simulation moves toward the observed statistics, then repeat.

At the MLE, the fitted distribution should make the observed sufficient statistics look typical, which is captured by the moment condition $E_{\hat\theta}[g(Y)] \approx g(y_{obs})$. This explains both the power and the difficulty of ERGM estimation. We are fitting the whole graph distribution by repeatedly generating graphs, so poor MCMC mixing or a specification that runs toward nearly empty and nearly complete networks cannot be separated from the estimation problem.

Convergence still does not establish substantive adequacy. It tells us the algorithm found a stable solution for the specified statistics. We then simulate again and ask whether the model reproduces features that matter for the argument but were not guaranteed by the fit.

4. Why Does This Matter for Applied Research?

Many social outcomes are relational and self-organizing. Alliances may be reciprocated, coalitions may expand through common partners, organizations may accumulate ties because they are already prominent, and segregation may emerge from many local choices. In those settings, dependence is not only a nuisance that makes standard errors inconvenient. It may be part of the social process the theory is trying to describe.

ERGMs let us express those ideas as competing features of a probability model rather than as informal descriptions of a graph. They can show that a covariate association changes once network structure is represented, or that an apparent structural pattern disappears once a rival process enters the model. They also force us to state the comparison we mean: closure relative to what density, degree structure, homophily, and risk set?

The limit is equally important. An ERGM coefficient tells us about conditional compatibility between a statistic and the observed network under the specification. It does not, on its own, identify a causal behavioral mechanism. That is why interpretation, simulation, and research design remain part of the analysis after the model converges.

The Applied Model Compass

Model	Question we want to understand	How it helps	What makes it distinct	What we should be able to say after fitting
Cross-sectional ERGM	Which actor, dyadic, and structural configurations make this graph more plausible?	Puts explicit graph statistics into a joint probability model	Treats network dependence as part of the estimand rather than only latent residual structure	“Tie toggles that close shared partners receive higher conditional probability after density, homophily, and degree structure are represented.”
Valued ERGM	What explains whether a relationship exists and how large it is?	Models counts or bounded tie values using an explicit reference measure	The support and baseline distribution are part of the model, not an afterthought	“Most dyads never exchange, but conditional on exchange occurring, observed volume is higher than the reference process predicts.”
TERGM	How are current ties and configurations associated with the preceding network?	Adds lagged dyad and network statistics to a transition distribution	Represents ordered network change rather than one graph snapshot	“Receiving a tie last period is associated with a higher probability of returning it this period, conditional on the rest of the transition model.”
STERGM	Do the processes that create ties differ from those that keep existing ties alive?	Gives formation and persistence separate risk sets and equations	Does not force absent and present dyads into one transition mechanism	“Common threats are associated with alliance formation, while institutionalization is associated with alliance persistence.”

The Advanced Route for Today

The lecture is organized around six questions:

When does a structural term represent the process we care about, and when is it absorbing a rival process?
How do we communicate an ERGM result without giving readers an uninterpretable log-odds table?
Which diagnostic problem are we seeing: data support, MCMC mixing, degeneracy, or poor substantive fit?
What changes when the outcome is a transition from one network to the next rather than one cross-section?
How do roster decisions, missing dyads, lag construction, and the number of transitions define the TERGM estimand?
When should formation and persistence receive separate equations, and when is a SAOM asking the more defensible question?

The cross-sectional material we will emphasize is interpretation on the friendship-probability scale, advanced failure diagnosis, and the ERGM versus latent-model comparison. The main new model begins at Temporal ERGMs.

ERGM Reference Bridge

The next derivation is included for review. The main path begins with the substantive specification problem in A1, then moves directly to B2.

A1 · Where We Left Off, and What Today Is For

On Friday (Day 10) you fitted latent-position models: the distance model as a useful comparison, then the latent factor and AME framework. In the longitudinal ICEWS application, the dynamic rank-two AME model represented changing country-level activity and pair-specific compatibility, and the bootstrap intervals made clear which time paths were stable enough to discuss. That model could absorb clustering into additive and multiplicative effects, but it did not hand us a coefficient called “closure.” Even good prediction would not turn its latent surface into a directly parameterized shared-partner mechanism.

Today is the opposite premise. Instead of absorbing the dependence through latent terms in the conditional mean, we make it the substance. We write down the probability of the whole graph as a function of the structural features we have a theory about: edges, triangles, and shared partners. We get a parameterized, testable coefficient on the closure term AME absorbs, though one cross-section cannot determine whether that term represents transitivity or homophily on an unobserved attribute (we return to this in the wrap). On the course map Day 8 drew, this is the “put it in the model” cell: the dependence is not a nuisance to soak up, it is the quantity of interest.

So let me be precise about what this day is not. It is not another explanation of why logit fails on network data. Earlier material established that dyads are dependent, and Day 9 made that problem formal through the SRM. The question now is what you do about dependence when it is the thing you want to measure.

The Substantive Version of the Question

Suppose we observe an alliance, legislative cosponsorship, or rebel-cooperation network and find substantial closure. One possible explanation is that partners introduce partners or that collaboration becomes easier through a trusted intermediary. The same triangles can also arise because actors share ideology, geography, institutional membership, or simply have many ties.

The ERGM question is therefore conditional: after accounting for density, actor characteristics, dyadic covariates, and other declared network features, is the observed amount of closure compatible with a positive closure parameter? That is sharper than saying “the network has many triangles,” but it is still an associational model of one observed graph unless the design supplies a stronger interpretation.

The applied workflow is:

State the hypothesized relational process, such as reciprocity, homophily, preferential attachment, or closure.
Translate that process into a graph statistic.
Identify rival processes that can produce the same visible pattern.
Put those rivals into the same probability model where the data support them.
Interpret the focal coefficient conditionally through a change statistic.
Simulate from the fitted model and check whether the model reproduces features the argument needs but the fit was not guaranteed to match.

The question to hold onto is: how clustered is faux.mesa.high relative to a clearly defined comparison? The answer depends entirely on that comparison. Below I take the network’s observed global transitivity and divide it by what we would expect under two null models: an Erdős–Rényi null (same number of nodes, same number of edges, ties placed at random) and a degree-preserving null (same nodes, same degree sequence, so each student keeps their observed number of friends, but partners are otherwise rewired at random).

nulls <- fit_or_load("a1_nulls", {
  # igraph is used with :: prefixes only; attaching it would mask network/sna.
  A   <- as.matrix.network(faux.mesa.high)
  g   <- igraph::graph_from_adjacency_matrix(A, mode = "undirected")
  obs <- igraph::transitivity(g, type = "global")
  n   <- igraph::gorder(g); m <- igraph::gsize(g)
  deg <- igraph::degree(g)
  set.seed(6886)
  er <- replicate(500, igraph::transitivity(
          igraph::sample_gnm(n, m), type = "global"))
  # Rewire the observed simple graph so every draw preserves the degree
  # sequence exactly. Sampling a multigraph and then simplifying it would not.
  dp <- replicate(500, igraph::transitivity(
          igraph::rewire(g, with = igraph::keeping_degseq(niter = 20 * m)),
          type = "global"))
  list(obs = obs, er_ratio = obs / mean(er), dp_ratio = obs / mean(dp))
})
data.frame(quantity = c("observed transitivity",
                        "ratio vs Erdos-Renyi null",
                        "ratio vs degree-preserving null"),
           value = round(c(nulls$obs, nulls$er_ratio, nulls$dp_ratio), 2))
#>                          quantity value
#> 1           observed transitivity  0.28
#> 2       ratio vs Erdos-Renyi null 27.37
#> 3 ratio vs degree-preserving null 11.71

You do not need to run this; the ratios are the point, not the igraph null-model code.

Against an Erdős–Rényi null the network is about 27× more transitive than chance. Against the degree-preserving null, which already knows how many friends each student has, it is about 12×. Same network, same statistic, two nulls, and the number more than halves: from 27× against the Erdős–Rényi null to 12× against the degree-preserving one. The Erdős–Rényi comparison was crediting triad closure with work that the degree sequence alone already does.

The choice of null model does a large part of the work in any “more clustered than chance” claim, and the first thing an ERGM buys you is control over exactly that null.

But notice what did not happen: it did not flip to “unremarkable.” Even against the demanding degree-preserving null, Mesa is still about an order of magnitude more clustered than the rewired draws. That residual pattern is large relative to this particular null, and modelling it is most of what today is about. ERGMs are the machine for building the null you actually want. You write down a probability distribution over entire graphs, you put in the features you think generated the network, such as degree, homophily, and closure, and then you ask whether the features you care about survive once the others are accounted for.

Two Nulls, Two Different Stories

We are going to spend most of today on one network, and it is a school.

data(faux.mesa.high)
mesa <- faux.mesa.high
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

This is faux.mesa.high, a simulated-but-realistic adolescent friendship network built to match an Add Health school in Mesa, Arizona. 205 students, 203 friendships, undirected. Every student has a Grade (7–12), a Sex, and a Race.

A school example keeps the actors and relationship easy to identify: students and friendships. That lets us spend our time on the model rather than first learning an unfamiliar institutional setting. The example also gives us recognizable actor attributes, including grade, sex, and race, so we can ask concrete questions about homophily, unequal levels of social activity, and clustering while remaining careful about what this simulated network can establish.

par(mar = c(0,0,1,0))
set.seed(6886)
plot(mesa, vertex.col = "Grade", vertex.cex = 0.9,
     main = "faux.mesa.high: vertex color = Grade")

Three plain questions, in the order you would actually ask them:

Do students befriend others in their own grade?
Do they befriend others of the same sex?
Are friendships transitive? Is my friend’s friend more likely to be my friend?

The first two are questions about node attributes. You could, if you squinted, answer them with a logistic regression on dyads and a lot of hand-waving about standard errors. The third one you cannot, because it is a statement about triples, and there is no way to write it as a property of dyad $ij$ alone. That asymmetry is why this day exists.

checkpoint(n = network.size(mesa),
           edges = network.edgecount(mesa),
           density = network.density(mesa),
           label = "CHECKPOINT A1")
#> ------------------------------------------------------------------
#> CHECKPOINT A1: n = 205   |   edges = 203   |   density = 0.01
#> ------------------------------------------------------------------

Note that density: 0.0097. Under one percent of possible friendships exist. Hold that number; it comes back in B2 and it changes how you report everything.

A2 · The Joint Model, and the Constant That Ruins Your Life

Here is the model. It is one line.

\[\Pr(Y = y \mid \theta) \;=\; \frac{\exp\{\theta^{\top} g(y)\}}{\kappa(\theta, \mathcal{Y})}\]

$\theta^{\top} g(y)$ is just each coefficient times its statistic, summed into a single weighted score for the whole graph. Read it slowly, because every practical property of the software falls out of it:

$Y$ is the whole network, treated as a single random variable. Not a collection of dyads. One draw from this distribution is one entire graph.
$y$ is the graph you actually observed.
$g(y)$ is a vector of network statistics, or counts of configurations. How many edges, how many same-grade edges, how many triangles.
$\theta$ is the vector of parameters you want.
$\kappa(\theta, \mathcal{Y})$ is the normalizing constant: the sum of $\exp\{\theta^\top g(y)\}$ over every graph in the sample space $\mathcal{Y}$.

That is an exponential family. If $g(y)$ were a row of covariates for an independent observation, this would be logistic regression and you would be done in a second. It is not, and you are not.

$g(y)$ Is Just Counting, and You Can Do It Right Now

The most useful thing to know about summary() in ergm is that when you hand it a formula, it counts, it does not fit. No estimation happens. This is how you inspect $g(y)$ before you ever run a model.

summary(mesa ~ edges + nodematch("Grade") +
          nodematch("Sex") + triangle)
#>           edges nodematch.Grade   nodematch.Sex        triangle 
#>             203             163             132              62

So: 203 edges; of those, 163 join students in the same grade and 132 join students of the same sex; and there are 62 triangles. That vector, 203, 163, 132, 62, is $g(y)$ for that model.

Stuck?: summary() on a formula vs summary() on a fit

These are two completely different operations that share a name, and mixing them up is the single most common early confusion:

summary(mesa ~ edges + triangle): mesa ~ ... is a formula. You get counts. Instant. Nothing is estimated. This is $g(y)$.
summary(m1) where m1 <- ergm(...): m1 is a fitted model. You get a coefficient table with standard errors. This is $\hat\theta$.

If you ever see a coefficient table where you expected counts, you accidentally fed it a fit. If you see counts where you expected coefficients, you forgot to run ergm().

And Now the Constant

$\kappa$ is a sum over the sample space. How big is the sample space? For an undirected network on $n$ nodes there are $n(n-1)/2$ possible ties, each present or absent, so there are $2^{n(n-1)/2}$ possible graphs.

n     <- network.size(mesa)
ndyad <- n * (n - 1) / 2
cat("Possible ties in Mesa:",
    format(ndyad, big.mark = ","), "\n")
#> Possible ties in Mesa: 20,910
cat("Graphs in the sample space: 2^",
    format(ndyad, big.mark = ","), "\n", sep = "")
#> Graphs in the sample space: 2^20,910
cat("For scale, atoms in the observable universe ~ 2^266\n")
#> For scale, atoms in the observable universe ~ 2^266

$2^{20{,}910}$. You are not summing that. Nobody is summing that, ever, on any hardware.

ergm does four things, and this is the working summary you need to use the software competently:

We cannot compute $\kappa$, so we cannot compute the likelihood directly.
But we can simulate networks from the model at a given $\theta$, using MCMC.
So: guess $\theta$, simulate a pile of networks, compare their statistics to the observed $g(y)$, move $\theta$ in the direction that closes the gap, repeat.
That is Monte Carlo maximum likelihood estimation (MCMC-MLE).

And now every annoying property of ergm is explained. Fits are slow because each iteration simulates networks. Fits are stochastic: run it twice without a seed and you get slightly different numbers, which is why every model today sets one. And fits can fail, because if the simulated networks do not look anything like your observed network, the comparison in step 3 carries no information and the optimizer has nothing to walk toward. That last failure has a name, degeneracy, and we will produce one on purpose in C1.

Depth: what MCMC-MLE is actually doing, and why I am not deriving it

Optional; nothing later depends on it.

The log-likelihood ratio between a candidate $\theta$ and a reference $\theta_0$ is

\[\ell(\theta) - \ell(\theta_0) = (\theta - \theta_0)^{\top} g(y) - \log \mathrm{E}_{\theta_0}\!\left[\exp\{(\theta-\theta_0)^{\top} g(Y)\}\right]\]

And the expectation on the right is approximated by an average over networks simulated at $\theta_0$ (Geyer & Thompson 1992). The approximation is only trustworthy when $\theta$ is close to $\theta_0$, which is why estimation iterates.

I am not deriving this in class because modern ergm no longer follows the simplest textbook version. As of 4.0, MCMLE.termination = "confidence" uses a multivariate equivalence test rather than a fixed iteration count, and MCMC.samplesize defaults to NULL, meaning the sample size is chosen adaptively to hit a target effective sample size (Vats, Flegal, and Jones 2019). A short whiteboard derivation would describe an outdated implementation. The current defaults are shown below:

f <- formals(control.ergm)
str(list(MCMLE.termination   = eval(f$MCMLE.termination),
         MCMLE.maxit         = eval(f$MCMLE.maxit),
         MCMC.samplesize     = eval(f$MCMC.samplesize),
         MCMLE.effectiveSize = eval(f$MCMLE.effectiveSize)))
#> List of 4
#>  $ MCMLE.termination  : chr [1:5] "confidence" "Hummel" "Hotelling" "precision" ...
#>  $ MCMLE.maxit        : num 60
#>  $ MCMC.samplesize    : NULL
#>  $ MCMLE.effectiveSize: num 64

Practical consequence, and it is the opposite of the folklore: do not hand-tune these. I have measured the same model on this dataset at 53 seconds using defaults and at over 24 minutes with MCMC.samplesize = 10000. The expensive one was worse because setting it by hand overrides the adaptive targeting that exists precisely to pick it for you. If you inherit a script full of inflated MCMC controls, the first thing to try is deleting them.

gy <- summary(mesa ~ edges + nodematch("Grade") + triangle)
checkpoint(gy = gy, label = "CHECKPOINT A2")
#> ------------------------------------------------------------------
#> CHECKPOINT A2: gy = 203, 163,  62
#> ------------------------------------------------------------------

A3 · Change Statistics: The Idea the Whole Day Rests On

If you remember one thing from today, this is the one. It is also the single most misunderstood point in applied ERGM work, and you will catch published papers getting it wrong once you have it.

$\theta$ lives in a model over whole graphs. But we want to say something about a tie. The bridge between them is the change statistic.

For dyad $ij$, define

\[\delta(y_{ij}) \;=\; g(y^{+}_{ij}) \;-\; g(y^{-}_{ij})\]

In words: take your network, force tie $ij$ on, compute $g$. Take the same network, force tie $ij$ off, compute $g$. Subtract. Everything else in the graph is held exactly as observed.

Then, and this is the result, it follows from the model in A2 by pure algebra, no approximation:

\[\log \frac{\Pr(Y_{ij} = 1 \mid Y^{c}_{ij} = y^{c}_{ij})}{\Pr(Y_{ij} = 0 \mid Y^{c}_{ij} = y^{c}_{ij})} \;=\; \theta^{\top} \delta(y_{ij})\]

An ERGM coefficient is the contribution to the conditional log-odds of a tie, given the rest of the graph. The phrase “given the rest of the graph” is not decoration and it is not a technicality. $Y^c_{ij}$ means every other dyad in the network, held at its observed value.

And that is the naive dyadic logit, the “just run a logistic regression on the dyads” that Day 9 opened by warning you against. Same object. The naive approach writes it down and estimates it directly, pretending the rest of the graph is fixed and independent; today it is instead a derived consequence of a single model over whole graphs.

Let Us Actually Compute One

No hand-waving. Pick two students and do the arithmetic.

grade <- get.vertex.attribute(mesa, "Grade")
sex   <- get.vertex.attribute(mesa, "Sex")
i <- 1; j <- 2
cat("student", i, ": grade", grade[i], "sex", sex[i], "\n")
#> student 1 : grade 7 sex F
cat("student", j, ": grade", grade[j], "sex", sex[j], "\n")
#> student 2 : grade 7 sex F
cat("tie currently present?", mesa[i, j] == 1, "\n")
#> tie currently present? FALSE

Now build the two networks, one with the tie forced on and one with it forced off, and count.

delta_ij <- function(net, i, j, form) {
  count <- function(g) {              # Count g(.) on one net
    ff <- statnet.common::nonsimp_update.formula(form, g ~ .)  # Plumbing: ignore on first read
    environment(ff) <- environment()    # Plumbing: ignore on first read
    summary(ff)
  }
  on  <- net; on[i, j]  <- 1                   # Force the tie ON
  off <- net; off[i, j] <- 0                   # Force the tie OFF
  count(on) - count(off)              # All else held fixed
}
f <- ~ edges + nodematch("Grade") + nodematch("Sex")
d12 <- delta_ij(mesa, 1, 2, f)
d12
#>           edges nodematch.Grade   nodematch.Sex 
#>               1               1               1

Only three lines carry the idea: on[i,j] <- 1, off[i,j] <- 0, and count(on) - count(off). The two formula lines are just machinery so summary() can see a modified network.

Student 1 and student 2 are both in grade 7 and both F, so adding the tie between them adds one edge, one same-grade edge, and one same-sex edge. The change statistic is $(1, 1, 1)$.

Now a pair who match on nothing:

k <- which(grade != grade[1] & sex != sex[1])[1]
cat("student", k, ": grade", grade[k], "sex", sex[k], "\n")
#> student 3 : grade 11 sex M
delta_ij(mesa, 1, k, f)
#>           edges nodematch.Grade   nodematch.Sex 
#>               1               0               0

$(1, 0, 0)$. It is still an edge, but it earns nothing from either homophily term.

So the difference in conditional log-odds between those two pairs is exactly $\theta_{\text{grade}} + \theta_{\text{sex}}$. That is what those coefficients mean. Not “the effect of grade on friendship” in some marginal, averaged-over-everything sense. It is the contribution to the log-odds of this tie, with the other 20,909 dyads pinned where you found them.

Stuck?: why does edges always have a change statistic of 1?

Because the edges term counts edges. Turning any tie on adds exactly one edge, turning it off removes exactly one. So $\delta_{\text{edges}} = 1$ for every dyad in the network, always.

That is why edges behaves like an intercept: it is the term whose change statistic never varies. And it is why a model with only edges is exactly an Erdős–Rényi graph, with one constant log-odds for every dyad. We will use that in B2.

Depth: the change statistic for a dependent term is not constant, and that is the whole story

edges, nodematch, nodecov all have change statistics that depend only on the two nodes involved. That is exactly the definition of dyad independence, and it is why models built only from those terms collapse to ordinary logistic regression.

triangle does not. The change statistic for triangle on dyad $ij$ is the number of nodes $k$ such that both $ik$ and $jk$ are present, which depends on the current state of other dyads. Watch it vary:

tri <- sapply(2:13, function(v) delta_ij(mesa, 1, v, ~ triangle))
setNames(tri, paste0("j=", 2:13))
#>  j=2  j=3  j=4  j=5  j=6  j=7  j=8  j=9 j=10 j=11 j=12 j=13 
#>    2    0    0    0    0    0    0    0    0    0    0    0

For most pairs, adding the tie closes zero triangles. For a few, it closes one or more. That variation is dependence, made concrete. It is also, as we will see in C1, exactly what makes the model dangerous, because a term whose change statistic can grow with the network is a term that can run away.

checkpoint(delta_1_2 = d12, label = "CHECKPOINT A3")
#> ------------------------------------------------------------------
#> CHECKPOINT A3: delta_1_2 = 1, 1, 1
#> ------------------------------------------------------------------

Advanced ERGM Work: Interpretation and Specification

How ERGM Estimation Works

For a dyad-independent ERGM, estimation reduces to familiar logistic regression. Once endogenous statistics such as shared partners enter, the normalizing constant depends on an enormous set of possible graphs, so the model uses Monte Carlo maximum likelihood.

The fitted parameter vector maximizes the graph log-likelihood

\[ \ell(\theta;y_{\mathrm{obs}}) =\theta^{\mathsf T}g(y_{\mathrm{obs}}) -\log\kappa(\theta), \qquad \kappa(\theta)=\sum_{y'\in\mathcal{Y}}\exp\{\theta^{\mathsf T}g(y')\}. \]

The first term rewards parameters that place weight on the statistics observed in our network. The normalizing term forces the model to compare that network with every other graph in the permitted sample space. We cannot enumerate those graphs in a realistic dependent model, so MCMC supplies a simulated approximation. Differentiating the log-likelihood gives the score

\[ \nabla\ell(\theta) =g(y_{\mathrm{obs}})-E_{\theta}[g(Y)]. \]

At the maximum, the observed statistics and their model-based expectations should be approximately equal. That moment-matching statement is not a separate estimator. It is the first-order condition for maximum likelihood in an ERGM.

The MCMLE loop is:

Find a starting value. A pseudolikelihood fit often supplies a computational starting point. It is not automatically the final inferential answer.
Simulate graphs at the current parameter vector. The MCMC sampler toggles ties and generates networks from the current ERGM.
Compare simulated and observed statistics. If simulated graphs have too little reciprocity or closure, the corresponding parameter needs to move upward; if they have too much, it needs to move downward.
Update the parameter vector. The likelihood approximation uses the simulated sufficient statistics to take a step toward the observed graph.
Repeat simulation and updating. The process continues until the observed statistics are plausible under graphs generated near the fitted parameter vector.
Estimate uncertainty. The variability and covariance of simulated statistics determine the curvature of the likelihood and therefore the coefficient covariance matrix.
Diagnose the chain and the model separately. MCMC diagnostics ask whether estimation worked. Goodness of fit asks whether the resulting model reproduces features relevant to the research question.

This is why ergm() can take a few seconds or many hours. Each likelihood update requires a network simulation problem, and a poorly specified model may generate graphs far from the observed one.

B1 · Reference: The Dyad-Independent Model and a Confusing Console Message

Let us fit something. Start with the two questions we can answer without any dependence at all: grade homophily and sex homophily. I will add nodecov("Grade") too, which asks whether older students simply have more friends.

m1 <- fit_or_load("m1_dyadindep",
  ergm(mesa ~ edges + nodematch("Grade") + nodematch("Sex") +
         nodecov("Grade"),
       control = control.ergm(seed = 6886)))

summary(m1)
#> Call:
#> ergm(formula = mesa ~ edges + nodematch("Grade") + nodematch("Sex") + 
#>     nodecov("Grade"), control = control.ergm(seed = 6886))
#> 
#> Maximum Likelihood Results:
#> 
#>                 Estimate Std. Error MCMC % z value Pr(>|z|)    
#> edges           -8.00333    0.50721      0 -15.779  < 1e-04 ***
#> nodematch.Grade  2.94996    0.18008      0  16.381  < 1e-04 ***
#> nodematch.Sex    0.64509    0.14952      0   4.314  < 1e-04 ***
#> nodecov.Grade    0.08908    0.02598      0   3.429 0.000605 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#>      Null Deviance: 28987  on 20910  degrees of freedom
#>  Residual Deviance:  1909  on 20906  degrees of freedom
#>  
#> AIC: 1917  BIC: 1949  (Smaller is better. MC Std. Err. = 0)

Under a second. No MCMC, no diagnostics needed, no chance of degeneracy. Why?

is.dyad.independent(m1)
#> [1] TRUE

Every term in that model has a change statistic that depends only on nodes $i$ and $j$. So the dyads are conditionally independent, the model factors into a product over dyads, and the whole thing is exactly a logistic regression on the $\binom{205}{2} = 20{,}910$ dyads.

The Thing That Confuses Everyone

When estimated from scratch, ergm reports “Starting maximum pseudolikelihood estimation (MPLE)… Finished MPLE.” The saved model object, however, labels the estimate as follows:

m1$estimate
#> [1] "MLE"

It says MLE. So which is it?

Both, and they are the same thing here. The pseudolikelihood is the product of the conditional probabilities of each dyad given the rest. When the model is dyad-independent, those conditionals are independent, so the pseudolikelihood is the likelihood, not an approximation to it, literally it. ergm runs the MPLE machinery because it is fast, and then correctly labels the result MLE because that is what it is.

When the model is dyad-dependent, that equivalence breaks and MPLE becomes an approximation with understated standard errors.

Depth: which estimator, and when: the three-line rule

This is the decision rule I would actually give you, and it has changed since the 2008 JSS paper you were assigned:

Dyad-independent model → it is exact logistic regression. Instant, no MCMC, no diagnostics. Nothing can go wrong.

Dyad-dependent, up to roughly a couple of thousand nodes → MCMC-MLE, which is the default. Let ergm choose its own sample size.

Bigger than that → MPLE, but you must change how the standard errors are computed. Schmid & Desmarais (2017) show MPLE point estimates are fine while the naive standard errors are badly understated; their parametric bootstrap gives correct coverage at about a fifth of the cost of MCMLE. This is now built in, but it is not the default:

eval(formals(control.ergm)$MPLE.covariance.method)
#> [1] "invHess"

invHess is the one that understates. If you are using MPLE on a dependent model, set MPLE.covariance.method = "bootstrap" (or "Godambe") and say in your paper that you did.

Contrastive divergence (estimate = "CD") exists, and Krivitsky (2017) introduced it precisely as a way to seed MCMLE. The literature is mixed. Current practical guidance is to use MPLE, the default seed, unless you have a specific reason to reach for CD, such as a curved model that is hard to initialize. If you use CD, check convergence carefully.

# Questions and Notes

checkpoint(coefs = round(coef(m1), 3), label = "CHECKPOINT B1")
#> ------------------------------------------------------------------
#> CHECKPOINT B1: coefs = -8.003,  2.950,  0.645,  0.089
#> ------------------------------------------------------------------

B2 · What the Numbers Mean for Students and Friendships

We have nodematch.Grade $\approx$ 2.95. The task is to turn that log-odds coefficient into a statement about two students and the chance that they are friends.

Step One: The Baseline

Start with the simplest possible reading. A model with only edges has a constant change statistic of 1 for every dyad, so the log-odds of any tie is just that one coefficient, and:

m0 <- fit_or_load("m0_edges",
  ergm(mesa ~ edges, control = control.ergm(seed = 6886)))
c(coef_edges = coef(m0)[["edges"]],
  implied_prob = plogis(coef(m0)[["edges"]]),
  observed_density = network.density(mesa))
#>       coef_edges     implied_prob observed_density 
#>     -4.625021107      0.009708274      0.009708274

plogis() of the edges coefficient recovers the density exactly. plogis() is the inverse-logit, the $\mathrm{logit}^{-1}$ that appears in the formulas above. That is a useful sanity check and it makes concrete that edges is the intercept.

Step Two: Why You Cannot Just Exponentiate and Stop

The tempting move is $e^{2.95} \approx 19$: “same-grade pairs have 19 times the odds.” That is arithmetically true and it is a conditional odds ratio, so it is defensible. But odds ratios on a network this sparse are close to useless for conveying magnitude to a reader, and there is a second problem that is worse.

The rule that matters most for reading other people’s papers

Do not compare ERGM coefficients across nested specifications.

This is Duxbury (2023, Sociological Methods & Research). An ERGM is a logistic-form model with a fixed residual variance, so, exactly as in ordinary logit, adding a term rescales all the other coefficients, even terms that are uncorrelated with it. So when a paper shows you model 1, model 2, model 3 and says “the homophily effect shrinks once we control for triadic closure,” that sentence is not supported by the numbers shown. The coefficient change blends a substantive shift with a pure rescaling.

You will see this in print constantly. Here is how you catch a broken ERGM paper, and in this literature you will, constantly: it is one of the six questions on the referee checklist at the end of today.

Step Three: Average Marginal Effects

The fix is to put everything on the probability scale, where rescaling is not an issue and where the number means something. This is one of the most useful parts of the day for communicating your results to people who do not do this.

The average marginal effect of term $k$ is

\[\mathrm{AME}_k \;=\; \frac{1}{|\mathcal{D}|}\sum_{ij \in \mathcal{D}} p_{ij}(1 - p_{ij}) \cdot \theta_k , \qquad p_{ij} = \mathrm{logit}^{-1}\!\left(\theta^{\top}\delta(y_{ij})\right)\]

This is the ordinary logit AME, computed dyad by dyad using the change statistics from A3. It is three lines of code rather than a black box. Read them once now, line by line tonight:

ame <- function(fit, term) {
  mp <- ergmMPLE(fit$formula, output = "dyadlist",
                 basis = fit$network)
  th <- coef(fit)
  X  <- mp$predictor[, names(th), drop = FALSE]  # Delta, all dyads
  p  <- plogis(as.vector(X %*% th))              # Fitted tie prob
  s  <- mean(p * (1 - p))                        # Logistic scale
  c(AME = s * th[[term]], scale = s)
}

ergmMPLE(..., output = "dyadlist") hands you the change statistic matrix for every dyad: the same $\delta(y_{ij})$ we computed by hand in A3, just for all 20,910 of them at once. That is the only non-obvious line.

Now watch what happens when we apply it to a dense network and a sparse one. First the monastery, with 18 monks and Sampson’s classic data. It is dense enough that a probability change is legible:

data(sampson)
s1 <- fit_or_load("s1_sampson",
  ergm(samplike ~ edges + nodematch("group") + mutual,
       control = control.ergm(seed = 6886)))
data.frame(estimate = round(coef(s1), 3))
#>                 estimate
#> edges             -2.271
#> nodematch.group    2.033
#> mutual             1.426
ame(s1, "nodematch.group")
#>       AME     scale 
#> 0.2681905 0.1318865

On average, two monks being in the same faction raises the probability of a nomination by about 27 percentage points. That sentence you can put in an abstract. “2.03” you cannot.

Now the same computation on Mesa:

ame(m1, "nodematch.Grade")
#>         AME       scale 
#> 0.027585167 0.009351022
c(density = network.density(mesa))
#>     density 
#> 0.009708274

Here is a puzzle, and I want you to feel it before I resolve it. Mesa’s same-grade coefficient is larger than Sampson’s: 2.95 vs 2.03, so which network has the bigger effect on the probability scale?

The answer is Sampson, and it is not close: 2.8 percentage points on Mesa against Sampson’s 27, an AME about 10 times smaller despite the larger coefficient. That is not a contradiction and it is not a bug; it is the most important practical lesson about marginal effects on networks.

Look at the scale term the function returned, $\overline{p(1-p)}$: it is 0.132 on Sampson and only 0.0094 on Mesa, a factor of about 14. In a network with density under one percent, essentially every dyad sits far out in the flat tail of the logistic curve, where even a large shift in log-odds moves the probability barely at all. A nineteen-fold change in odds against a baseline of 0.0097 is still a small absolute probability.

So report both, always

On a sparse network an AME reported alone can make a meaningful effect look trivial, while a coefficient reported alone can make it look enormous. Report the AME next to the baseline density. “Same-grade pairs are 2.8 percentage points more likely to be friends, against a baseline density of 0.97 percent” is complete and interpretable. Together those two numbers tell you the shift is about 2.8× the baseline density itself. That is substantively large while still being small on an absolute scale. Either half on its own is misleading.

Depth: two caveats on that AME, and the packaged version

First, the change statistics come from ergmMPLE evaluated at the observed network. For a dyad-independent model that is exactly right. For a dependent model it is a conditional-on-the-observed-graph marginal effect, not a marginal effect under the model’s own equilibrium distribution. It is the standard thing everybody reports, and it is what the package below reports too, but it is worth knowing what it conditions on.

Second, the formula above is the derivative marginal effect. For a binary term you could instead take the discrete difference, $\overline{\mathrm{logit}^{-1}(\eta + \theta_k) - \mathrm{logit}^{-1}(\eta)}$, where $\eta$ is the linear predictor with that term switched off.

These agree when effects are small and diverge when they are not: on Sampson the derivative version gives 0.268 and the discrete version gives about 0.38. Neither is wrong; say which you used.

The packaged version is ergMargins::ergm.AME() (Duxbury), which computes exactly the derivative formula above plus a delta-method standard error. I have checked our three lines against it and they agree to five decimal places. I am teaching the arithmetic rather than the package for two reasons: you can see what it does, and it is one less thing to install for the people. If you want the standard errors and the mediation/moderation tools, install.packages("ergMargins"); it is on CRAN, and vif.ergm() in the same package is a cheap way to diagnose the near-collinearity that produces “linear dependence detected among model statistics” warnings.

Tiered Exercise

Work through the core prompt before opening the solution tab. Everything is instant and seeded, so your numbers will match the displayed results.

Compute the change statistic for the nodecov("Grade") term, by hand, for the pair (1, 2), and confirm it against summary(). Fill in the blank:

delta_ij(mesa, 1, 2, ~ ______)

Then say in one sentence what nodecov("Grade") tests that nodematch("Grade") does not.

The nodecov("Grade") coefficient is small (0.089) but strongly significant. Compute its AME. Then explain why an AME is a slightly odd thing to report for this particular term, and what you would report instead.

delta_ij(mesa, 1, 2, ~ nodecov("Grade"))
#> nodecov.Grade 
#>            14

Core. nodecov("Grade") sums the grade values of the two endpoints, so its change statistic is $\text{grade}_i + \text{grade}_j$ = 7 + 7 = 14. nodematch asks whether similarity drives ties: do 9th graders befriend 9th graders? nodecov asks whether the level drives ties: do older students have more friends regardless of who those friends are. One is homophily, the other is activity, and a network can easily have one without the other.

Stretch.

ame(m1, "nodecov.Grade")
#>          AME        scale 
#> 0.0008330283 0.0093510216

The AME is tiny, but the deeper problem is that nodecov("Grade") is not binary; its change statistic ranges from 14 to 24. So “the average marginal effect” is the effect of a one-unit increase in the summed grades of the pair, which is not a quantity anyone has an intuition about. Better to report the effect of moving one endpoint up a full grade (multiply by 1, since one endpoint moving up one grade raises the sum by one), or better still, report predicted tie probabilities at a few substantively chosen grade combinations. This is the general lesson: AMEs are cleanest for binary terms and need thought for continuous ones.

checkpoint(ame_sampson = ame(s1, "nodematch.group")[["AME"]],
           ame_mesa    = ame(m1, "nodematch.Grade")[["AME"]],
           label = "CHECKPOINT B2")
#> ------------------------------------------------------------------
#> CHECKPOINT B2: ame_sampson = 0.268   |   ame_mesa = 0.028
#> ------------------------------------------------------------------

B3 · Endogenous Terms: The Reason You Are Here

Everything so far you could have done with a logistic regression and the right standard errors. Now we do the thing that only an ERGM does: put a structural feature of the network into the model as a parameter.

Reciprocity

The first and simplest endogenous term is mutual: if $i$ names $j$, is $j$ more likely to name $i$?

And here we hit a genuine constraint that is worth stating out loud, because it teaches something about how these terms work:

is.directed(mesa)
#> [1] FALSE

Mesa is undirected. A friendship there is a single symmetric object; there is no such thing as $i \to j$ without $j \to i$, so reciprocity is not a question you can ask. Try it and ergm will refuse. This is not an inconvenience; it is the model telling you that term availability is a property of your data, not a menu you choose from.

So reciprocity goes on Sampson, which is directed, with 18 monks each nominating the others.

round(summary(s1)$coefficients, 3)
#>                 Estimate Std. Error MCMC % z value Pr(>|z|)
#> edges             -2.271      0.228      0  -9.939    0.000
#> nodematch.group    2.033      0.304      0   6.692    0.000
#> mutual             1.426      0.478      0   2.985    0.003

mutual = 1.43, and its change statistic is 1 exactly when the reciprocal tie $ji$ is already present, 0 otherwise. So: conditional on everything else, a nomination that would reciprocate an existing one has 1.43 more log-odds than one that would not. On the probability scale:

ame(s1, "mutual")
#>       AME     scale 
#> 0.1881206 0.1318865

About 19 percentage points. Reciprocity is a big deal in this monastery, which, if you have read Sampson, is not a surprise.

The Term Vocabulary Is Searchable: Do Not Memorize It

edges, nodematch, nodecov, mutual, where do these names come from? A term database that ships inside ergm itself, several hundred entries deep across the Statnet packages. Nobody holds that in their head, and you should not try, because the package has a built-in search engine and using it is the single habit that separates people who fight this software from people who use it.

search.ergmTerms(search = "geometrically weighted edgewise")
#> Found  1  matching ergm terms:
#> dgwesp(decay, fixed=FALSE, cutoff=30, type="OTP") (binary)
#> gwesp(decay, fixed=FALSE, cutoff=30, type="OTP") (binary)
#>     Geometrically weighted edgewise shared partner distribution
search.ergmTerms(name = "nonzero")
#> No terms named 'nonzero' were found. Try searching with search='nonzero'instead.

Two lookups, two lessons. search = scans the descriptions and finds gwesp; hold that name, because it is the repair for the disaster we are going to cause on purpose in C1. And read its signature closely: fixed=FALSE is the default, which is a loaded gun we come back to in C1. name = wants an exact term name, and the second call misses, but the miss tells you the fix. (nonzero exists in the valued vocabulary, and it is the hinge of C3.)

Three more searches to run on your own time

search.ergmTerms(keywords = "triad-related"): all 17 closure-adjacent terms in one screen.
search.ergmTerms(keywords = "curved"): the geometrically weighted family; every single signature defaults to fixed = FALSE.
search.ergmConstraints("degree") and search.ergmReferences("discrete"): the same engine for the constraints = and reference = arguments, the other two vocabularies you will eventually need (the second one matters this afternoon, in C3).

Transitivity, and the Naive Way to Ask for It

Back to Mesa, and to question 3 from A1: are friendships transitive?

The obvious term is triangle, which literally counts triangles. We already know there are 62 of them. So the specification writes itself:

ergm(mesa ~ edges + triangle)

That is the natural thing to write, it is what the term is called, and it is what a reasonable person would try first.

The walkthrough leaves this specification unevaluated here. The next section shows the recorded failure from attempting the fit and uses that failure to explain degeneracy. Nothing needs to be estimated during class.

checkpoint(mutual_coef = coef(s1)[["mutual"]],
           mesa_triangles = summary(mesa ~ triangle),
           label = "CHECKPOINT B3")
#> ------------------------------------------------------------------
#> CHECKPOINT B3: mutual_coef = 1.426   |   mesa_triangles = 62
#> ------------------------------------------------------------------

The full toolbox: specification terms we did not need today

The core path uses the smallest set of terms that carries the ideas. These are the ones you will reach for first in applied work. Every chunk runs as-is; none is required for anything else today.

Descriptives before you fit. degreedist() is a quick look at the degree distribution, and isolates + degree(d) are the model-term versions of the same facts:

degreedist(mesa)
#>  degree0  degree1  degree2  degree3  degree4  degree5  degree6  degree7 
#>       57       51       30       28       18       10        2        4 
#>  degree8  degree9 degree10 degree13 
#>        1        2        1        1
summary(mesa ~ isolates + degree(0:3))
#> isolates  degree0  degree1  degree2  degree3 
#>       57       57       51       30       28

57 of 205 students are isolates. That is not a nuisance; it is a quarter of your nodes, and a model whose simulated networks have no isolates is wrong in a way gof()’s degree panel will catch (C2).

Homophily beyond nodematch. If your attribute is continuous, a referee will expect absdiff, the sum of $|x_i - x_j|$ over edges, so a negative coefficient means similarity attracts. And nodemix is the model-term version of the mixingmatrix() we use in C2: one parameter per attribute-pair cell instead of one per diagonal.

summary(mesa ~ absdiff("Grade"))
#> absdiff.Grade 
#>            79
summary(mesa ~ nodemix("Grade", levels2 = -1))
#>   mix.Grade.7.8   mix.Grade.8.8   mix.Grade.7.9   mix.Grade.8.9   mix.Grade.9.9 
#>               0              33               0               2              23 
#>  mix.Grade.7.10  mix.Grade.8.10  mix.Grade.9.10 mix.Grade.10.10  mix.Grade.7.11 
#>               1               4               7               9               1 
#>  mix.Grade.8.11  mix.Grade.9.11 mix.Grade.10.11 mix.Grade.11.11  mix.Grade.7.12 
#>               2               6               1              17               1 
#>  mix.Grade.8.12  mix.Grade.9.12 mix.Grade.10.12 mix.Grade.11.12 mix.Grade.12.12 
#>               1               4               5               5               6

Note the empty off-diagonal cells (mix.Grade.7.8 = 0, mix.Grade.7.9 = 0): fit nodemix here naively and you buy the same -Inf problem C2 demonstrates for Race: same disease, same mixingmatrix()-first cure. (levels2 = -1 drops the first cell as the reference; that default flipped in ergm 4; see Appendix A.)

Degree-based dependence. gwdegree is to popularity what gwesp is to closure, a geometrically weighted count that damps the runaway feedback, and it carries exactly the same silent fixed = TRUE trap from C1. Directed versions: gwidegree / gwodegree.

summary(mesa ~ gwdegree(0.5, fixed = TRUE))
#> gwdeg.fixed.0.5 
#>        199.5657

Dyadic covariates: the one this box exists for. Many familiar network questions, from whether trade follows alliances to whether wave-1 ties predict wave-3 ties, are dyadic covariate questions. edgecov(x) puts an entire matrix on the right-hand side: its change statistic for dyad $ij$ is simply $x_{ij}$. Sampson’s three waves make a clean demonstration. We will predict wave 3 from wave 1:

data(samplk)
ec <- fit_or_load("tb_edgecov",
  ergm(samplk3 ~ edges + mutual + edgecov(as.matrix(samplk1)),
       control = control.ergm(seed = 6886)))
round(summary(ec)$coefficients, 3)
#>                            Estimate Std. Error MCMC % z value Pr(>|z|)
#> edges                        -2.802      0.282      0  -9.921        0
#> mutual                        1.826      0.513      0   3.557        0
#> edgecov.as.matrix(samplk1)    2.676      0.379      0   7.053        0

A wave-1 nomination raises the log-odds of the same wave-3 nomination by 2.68, over and above reciprocity. That is a hand-rolled lagged-network term, and a two-line preview of C4, where the lag stops being hand-rolled and becomes the model. The sibling dyadcov(x) fits separate parameters by dyad type for directed data.

Term operators: the vocabulary became compositional in ergm 4. Appendix A lists them (B(), F(), Offset(), Symmetrize(), …); here is one actually running. F(formula, filter) evaluates a statistic on the subgraph of dyads that satisfy the filter:

summary(samplike ~ F(~edges, ~nodematch("group")) + nodematch("group"))
#> F(nodematch("group"))~edges             nodematch.group 
#>                          63                          63
f_op <- ergm(samplike ~ edges + F(~edges, ~nodematch("group")),
             control = snctrl(seed = 6886))
round(coef(f_op), 3)
#>                       edges F(nodematch("group"))~edges 
#>                      -2.001                       2.648

Both statistics are 63: “edges, restricted to within-group dyads” is the nodematch count, rebuilt from parts. So the fit is a reparametrized edges + nodematch("group") and its coefficients match that model exactly (-2.001 and 2.648). This becomes useful when no ready-made term exists. “gwesp among same-sex dyads only” is one F() away, not a feature request. One more thing in that chunk: control = snctrl(seed = 6886) is the ergm 4 Cross-function control shorthand. The same snctrl() call works as the control = for ergm(), simulate(), gof(), san(), and friends, and dispatches to the right control.*() constructor, so you can stop remembering which function wants control.gof.ergm() versus control.ergm().

Advanced Failure Diagnosis and Network Change

C1 · Diagnosing Degeneracy

The next chunk records the attempted fit. Its failure is part of the example rather than an unexpected classroom interruption.

m_tri <- ergm(mesa ~ edges + triangle,
              control = control.ergm(seed = 6886))
#> Error in `ergm.MCMLE()`:
#> ! 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.

Read that message, because it is one of the more useful error messages in R:

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.

Here is what happened. Remember step 3 of MCMC-MLE: simulate networks at the current $\theta$ and compare their statistics to the observed. ergm set out to do that, simulated from a model with a triangle term, and the simulated networks came back with more than twenty times as many edges as Mesa has. It then stopped, on purpose, rather than burn twenty minutes walking toward nothing.

Why Triangles Do This

The intuition is short and it is worth having. Adding a tie that closes a triangle makes further triangles easier to close, which makes further ties more attractive, which closes more triangles. With a positive triangle coefficient the model has a positive feedback loop in it, and the probability mass collapses onto the extremes of the sample space, the empty graph and the complete graph, with almost nothing in between. Observed social networks commonly fall between those extremes, so this specification puts essentially no probability on anything resembling the data.

Formally this is instability in the sense of Schweinberger (2011): statistics that can grow faster than linearly in the number of dyads make the log-likelihood dominated by extremal graphs. The term “degeneracy” gets used loosely; that is the precise version.

Statnet’s own framing, and it is the right one

“Unwanted degeneracy is an indicator of a poorly specified model. It is not a property of all ERGMs (nor is it unique to ERGMs!) … 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.” Source: the current Statnet ERGM workshop tutorial

This is Statnet’s position, and meaningfully softer than the 2011-era “ERGMs are fragile” story. The standing counter-view (Schweinberger 2011; Chatterjee & Diaconis 2013) is that for many natural parameterizations near-degeneracy is a structural feature, not merely user error, which is precisely why geometrically-weighted terms, and later LOLOG (see the wrap), had to be invented. Degeneracy is information about your specification. The model just told you that the process you wrote down could not have produced your data. That is a finding.

The Repair: Geometric Weighting

The problem is that each additional shared partner is worth as much as the first. The fix is a term where each additional shared partner is worth less than the last: a geometrically weighted count. That is gwesp: geometrically weighted edgewise shared partners.

m2 <- fit_or_load("m2_gwesp",
  ergm(mesa ~ edges + nodematch("Grade") + nodematch("Sex") +
         gwesp(0.25, fixed = TRUE),
       control = control.ergm(seed = 6886)))

summary(m2)
#> Call:
#> ergm(formula = mesa ~ edges + nodematch("Grade") + nodematch("Sex") + 
#>     gwesp(0.25, fixed = TRUE), control = control.ergm(seed = 6886))
#> 
#> Monte Carlo Maximum Likelihood Results:
#> 
#>                  Estimate Std. Error MCMC % z value Pr(>|z|)    
#> edges             -6.6027     0.1696      0  -38.94   <1e-04 ***
#> nodematch.Grade    1.9936     0.1661      0   12.00   <1e-04 ***
#> nodematch.Sex      0.5359     0.1288      0    4.16   <1e-04 ***
#> gwesp.fixed.0.25   1.4652     0.1051      0   13.94   <1e-04 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#>      Null Deviance: 28987  on 20910  degrees of freedom
#>  Residual Deviance:  1725  on 20906  degrees of freedom
#>  
#> AIC: 1733  BIC: 1764  (Smaller is better. MC Std. Err. = 0.71)

It converges. gwesp is 1.47, strongly positive: Mesa has substantially more triadic closure than grade and sex homophily alone can explain. That is a defensible finding about adolescent friendship, and it is exactly the statement that the naive triangle term could not deliver.

Notice also that grade homophily dropped from 2.95 to 1.99. Resist the urge to narrate that as “homophily was partly spurious closure.” That is precisely the across-specification comparison Duxbury warns about in B2. Part of that drop is rescaling. You cannot tell how much from these two numbers.

Two hard rules, both of which will bite you

1. Write fixed = TRUE unless you specifically intend a curved model. As of ergm 4.0, if you write gwesp(0.25) without fixed = TRUE, the decay value you typed is silently ignored and you get a curved model that estimates decay instead. It is not a subtle difference:

rbind(unfixed_0.25 = summary(samplike ~ gwesp(0.25))[1:5],
      unfixed_0.90 = summary(samplike ~ gwesp(0.90))[1:5])
#>              esp.OTP#1 esp.OTP#2 esp.OTP#3 esp.OTP#4 esp.OTP#5
#> unfixed_0.25        23        20        14        11         1
#> unfixed_0.90        23        20        14        11         1
c(fixed_0.25 = summary(samplike ~ gwesp(0.25, fixed = TRUE)),
  fixed_0.90 = summary(samplike ~ gwesp(0.90, fixed = TRUE)))
#> fixed_0.25.gwesp.OTP.fixed.0.25  fixed_0.90.gwesp.OTP.fixed.0.9 
#>                        80.57959                       108.08575

The two unfixed rows are byte-identical; the decay argument did nothing. The two fixed values differ, as they should. ergm does warn, but the warning scrolls past in a wall of MCMC output. I have found this bug in published replication code more than once.

2. Never load a cached ergm fit from an old version. summary.ergm compares the stored version against the current major.minor, so under 4.12.0 any object saved from any earlier version warns. There is no version to freeze on; this recurs at every release. It is also the single best argument for the format of this document: everything you are looking at was refit today.

Depth: where does 0.25 come from? It is a convention

That decay value is a convention, and the convention is under substantive scrutiny.

Stewart, Schweinberger et al. (Social Networks) show that geometrically weighted terms are nonlinear functions of the product of the base and decay parameters, so the two are “mixed up, and difficult to estimate,” and a single network “may not contain much information about decay parameters (in the statistical sense of Fisher information).” Their practical finding is uncomfortable: fixing decay at a value other than its MLE changes the estimates of all the other parameters, and can hurt both in-sample and out-of-sample performance. Their own estimated decay was 0.913, nowhere near the conventional 0 or 0.25.

So: everybody fixes decay because estimating it from one network is unstable, and the best current evidence says fixing it biases everything else. That is a live, unresolved tension and I would rather you knew about it than repeated “0.25 because that is what people use.” Their proposed remedy is multiple networks, which supply the missing information; see ergm.multi, and the ancillary document on that if you have data shaped that way.

checkpoint(gwesp = coef(m2)[["gwesp.fixed.0.25"]],
           grade = coef(m2)[["nodematch.Grade"]],
           label = "CHECKPOINT C1")
#> ------------------------------------------------------------------
#> CHECKPOINT C1: gwesp = 1.465   |   grade = 1.994
#> ------------------------------------------------------------------

C2 · Diagnostics: Separation First, Then Mixing, Then Fit

Three checks, in the order you should actually run them, followed by the one tool underneath all three. The order matters and almost nobody teaches it this way.

1. Look at Your Cells Before You Fit

We have not touched Race yet. The natural specification is nodematch("Race", diff = TRUE), which gives a separate homophily parameter per racial category. Before running it, look at the mixing matrix:

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

Look at the diagonal. Black–Black is 0 and Other–Other is 0: 6 Black students, 4 “Other” students, and none of them friends with each other. That is perfect separation: ask for a homophily parameter on a category where homophily never occurs and the maximum likelihood answer is negative infinity. Watch:

m_race <- ergm(mesa ~ edges + nodematch("Race", diff = TRUE),
               control = control.ergm(seed = 6886))
data.frame(estimate = round(coef(m_race), 3))
#>                      estimate
#> edges                  -4.826
#> nodematch.Race.Black     -Inf
#> nodematch.Race.Hisp     0.125
#> nodematch.Race.NatAm    0.944
#> nodematch.Race.Other     -Inf
#> nodematch.Race.White    1.208

Two -Inf coefficients, exactly the two empty diagonal cells. This is not a numerical failure; it is the right answer to the question you asked. The fix is to ask a better question: restrict the term to levels that can actually be estimated.

m_race2 <- ergm(
  mesa ~ edges + nodematch("Race", diff = TRUE,
                           levels = c("White", "Hisp", "NatAm")),
  control = control.ergm(seed = 6886))
data.frame(estimate = round(coef(m_race2), 3))
#>                      estimate
#> edges                  -4.828
#> nodematch.Race.White    1.210
#> nodematch.Race.Hisp     0.127
#> nodematch.Race.NatAm    0.946

Finite, interpretable, and you know exactly which groups you did and did not learn about. Checking mixingmatrix() before fitting can save you an afternoon of wondering why your model returned infinity. This generalizes well beyond nodematch; any categorical term with sparse cells has this problem.

2. Did the Chain Actually Mix?

For any dyad-dependent model, the coefficients are only as good as the MCMC that produced them. mcmc.diagnostics() shows you the sampled network statistics, centred on their observed values.

mcmc.diagnostics(m2, vars.per.page = 4)

#> Sample statistics summary:
#> 
#> Iterations = 192512:3805184
#> Thinning interval = 2048 
#> Number of chains = 1 
#> Sample size per chain = 1765 
#> 
#> 1. Empirical mean and standard deviation for each variable,
#>    plus standard error of the mean:
#> 
#>                    Mean    SD Naive SE Time-series SE
#> edges            10.551 27.97   0.6658          1.887
#> nodematch.Grade   9.793 26.25   0.6248          1.907
#> nodematch.Sex     7.388 19.67   0.4682          1.205
#> gwesp.fixed.0.25 11.388 32.75   0.7796          2.360
#> 
#> 2. Quantiles for each variable:
#> 
#>                    2.5%   25%   50%   75% 97.5%
#> edges            -41.00  -9.0 10.00 29.00 66.00
#> nodematch.Grade  -38.90  -8.0  9.00 28.00 62.90
#> nodematch.Sex    -29.00  -6.0  7.00 21.00 47.00
#> gwesp.fixed.0.25 -48.53 -11.9 10.17 33.91 79.31
#> 
#> 
#> Are sample statistics significantly different from observed?
#>                   edges nodematch.Grade nodematch.Sex gwesp.fixed.0.25
#> diff.      1.055071e+01    9.792635e+00  7.388102e+00     1.138796e+01
#> test stat. 5.590242e+00    5.134674e+00  6.129099e+00     4.825885e+00
#> P-val.     2.267530e-08    2.826331e-07  8.837825e-10     1.393830e-06
#>                  (Omni)
#> diff.                NA
#> test stat. 4.689976e+01
#> P-val.     3.299483e-09
#> 
#> Sample statistics cross-correlations:
#>                      edges nodematch.Grade nodematch.Sex gwesp.fixed.0.25
#> edges            1.0000000       0.9655049     0.9185038        0.9441814
#> nodematch.Grade  0.9655049       1.0000000     0.8823251        0.9519658
#> nodematch.Sex    0.9185038       0.8823251     1.0000000        0.8661509
#> gwesp.fixed.0.25 0.9441814       0.9519658     0.8661509        1.0000000
#> 
#> Sample statistics auto-correlation:
#> Chain 1 
#>               edges nodematch.Grade nodematch.Sex gwesp.fixed.0.25
#> Lag 0     1.0000000       1.0000000     1.0000000        1.0000000
#> Lag 2048  0.7283863       0.7692539     0.6996491        0.7435665
#> Lag 4096  0.5639869       0.6133855     0.5295487        0.5925610
#> Lag 6144  0.4542051       0.5053191     0.4104748        0.4937605
#> Lag 8192  0.3695392       0.4158056     0.3308802        0.4108785
#> Lag 10240 0.3126416       0.3534590     0.2721792        0.3495792
#> 
#> Sample statistics burn-in diagnostic (Geweke):
#> Chain 1 
#> 
#> Fraction in 1st window = 0.1
#> Fraction in 2nd window = 0.5 
#> 
#>            edges  nodematch.Grade    nodematch.Sex gwesp.fixed.0.25 
#>       -0.7754573       -0.5970973       -1.0614740       -0.7125730 
#> 
#> Individual P-values (lower = worse):
#>            edges  nodematch.Grade    nodematch.Sex gwesp.fixed.0.25 
#>        0.4380695        0.5504424        0.2884745        0.4761100 
#> Joint P-value (lower = worse):  0.5892412 
#> 
#> 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 you want: trace plots that look like fuzzy horizontal caterpillars centred near zero, and roughly symmetric densities. What you do not want: trends, long excursions, or a chain that wanders off and stays there. These traces look stable, and the Geweke check does not flag a within-chain drift, but the retained statistics are still fairly autocorrelated. Also notice that the mean-difference test printed here uses the last MCMLE simulation round, not a fresh simulation at the final estimate. The in-model gof() check below is the cleaner final test.

If the trace or Geweke check had looked poor, the moves are a longer burn-in, more separation between retained draws, and then a different specification. In practice the specification is usually the important part. Changing MCMC.interval is thinning, not evidence by itself that the underlying chain has mixed.

3. Goodness of Fit: Two Different Questions

gof() answers “if the model is true, does it generate networks that look like mine?” But there are two versions of that question and people routinely run only the second one.

In-model first. Did the model reproduce the statistics it was explicitly fitted on? This is GOF = ~model, it is fast, and it returns a table, which is often easier to inspect precisely than a plot:

set.seed(6886)
gof(m2, GOF = ~model)
#> 
#> Goodness-of-fit for model statistics 
#> 
#>                       obs       min     mean      max MC p-value
#> edges            203.0000 119.00000 196.4400 267.0000       0.64
#> nodematch.Grade  163.0000  88.00000 157.9200 223.0000       0.68
#> nodematch.Sex    132.0000  73.00000 127.6100 166.0000       0.76
#> gwesp.fixed.0.25 131.7582  60.54839 124.2299 197.0183       0.64

Each row: the observed value, the range across simulated networks, and a Monte Carlo p-value. High monte-carlo p-values mean the observed statistic sits comfortably inside the simulated distribution. Good, though note this is close to a minimum requirement rather than an achievement. A model that cannot reproduce its own sufficient statistics has not converged.

Out-of-model second, because this is the stronger check. Simulate from the fit and compare features you did not fit on. If the model captures the generative process, it should get those right too, for free.

g2 <- gof(m2)
plot(g2)

The default suite for an undirected network is degree distribution, edgewise shared partners, and geodesic distance. Black line observed, boxplots simulated: where the line sits inside the boxes, the model reproduces a feature it was never fitted to.

Here the degree and shared-partner panels are broadly compatible with the simulations. The geodesic panel finds the remaining weakness: Mesa has more connected pairs at path lengths 8 through 11 than the model usually generates. The fitted networks capture local friendship clustering but make the connected parts of the school too compact; they do not generate enough of the long chains that link students who are far apart in the network. That is a useful miss, not a reason to discard the whole fit. It tells us exactly which feature a richer specification would need to improve.

This is the same move you already know

This is exactly the posterior predictive logic from the social relations model day: simulate from the fitted model, compare to the observed data, and look for systematic discrepancy. Different implementation, identical epistemics. That is the third time this course has come back to it, and it is not an accident; simulation-based adequacy checking is close to the only general-purpose model criticism that works for network data.

4. And the Tool Underneath All of It: `simulate()`

gof() is a convenience wrapper around something you can, and should, do with your own hands: draw networks from the fitted model. The day opened with null models in A1; a fitted ERGM is a null model, tuned to your data, and simulate() is how you draw from it.

sims <- simulate(m2, nsim = 100, seed = 6886,
                 output = "stats", monitor = ~ triangle)
summary(mesa ~ triangle)
#> triangle 
#>       62
round(quantile(sims[, "triangle"], c(0.025, 0.5, 0.975)), 1)
#>  2.5%   50% 97.5% 
#>  36.4  62.0 100.5

Read the two arguments that do the work. output = "stats" says: do not keep the hundred networks, keep their statistics. monitor = ~ triangle says: also track a statistic that is not in the model: m2 has gwesp, not triangle. So this is the out-of-model GOF logic, hand-built in two lines: Mesa has 62 triangles, and the model we fitted, which was never told to match triangles, puts its median draw at 62, with the observed count sitting comfortably inside the middle 95% of 36.4 to 100.5. That is gwesp earning its keep on a statistic it was not fitted to.

Now connect this backwards: this is exactly the machine inside step 2 of MCMC-MLE from A2, and it is the machine the degeneracy guard in C1 was running when it simulated from edges + triangle, saw twenty times too many edges, and pulled the plug. One function, underneath the entire day.

Depth: what GOF cannot do, and why AIC is on thin ice here

gof() has two important limits.

It is a check, not a comparison. It tells you whether this model reproduces those statistics. It gives you no way to rank two competing models, and there is a mild circularity in that you choose the statistics to check.

And the usual model-selection tools are shakier than they look. ergm will happily print AIC and BIC. Yin, Phillips & Butts (arXiv:1908.05873) make the sharp version of the objection: those criteria “assume conditional independence of observations and/or conventional asymptotics, and their theoretical foundations are not always applicable in a network modeling context.” You have one network. There is no growing sample of independent observations. Their proposal is the HOPE scheme: held-out predictive evaluation, or cross-validation for networks.

Worth knowing: ergm computes the log-likelihood by bridge sampling and now reports a standard error on it. If you are comparing AICs that differ by less than that, you are comparing noise.

There is also a deeper problem in the background, Shalizi & Rinaldo (2013, Annals of Statistics): most network data are sampled subnetworks, and inference implicitly assumes the model is consistent under sampling (projective). Triangle-like terms violate this, so “given an ERGM for a large network, it cannot be the case that the induced model for any sub-network must be the same”; standard asymptotics do not hold. Practically, an ERGM fitted to a sampled network does not straightforwardly estimate parameters of the population network. This is not a reason to avoid ERGMs; it is a reason to be careful about what your network is.

Tiered Exercise

Read the prompt first, then open the solution tab. The fit is complete and seeded, so your numbers will match the displayed results without waiting on a new gwesp fit.

Add nodefactor("Grade") to m2 and refit. In one sentence, say what it tests that nodematch("Grade") does not.

m3 <- ergm(mesa ~ edges + nodematch("Grade") + nodematch("Sex") +
             gwesp(0.25, fixed = TRUE) + ______,
           control = control.ergm(seed = 6886))

Race has near-empty cells, as we saw. Diagnose with mixingmatrix() before fitting, then fit a model with nodematch("Race", diff = TRUE) restricted to estimable levels plus gwesp. Report the AME for the largest race-homophily effect, and state the baseline density alongside it. Then say why you would be cautious about the White estimate specifically.

m3 <- fit_or_load("m3_nodefactor",
  ergm(mesa ~ edges + nodematch("Grade") + nodematch("Sex") +
         gwesp(0.25, fixed = TRUE) + nodefactor("Grade"),
       control = control.ergm(seed = 6886)))
data.frame(estimate = round(coef(m3), 3))
#>                     estimate
#> edges                 -6.961
#> nodematch.Grade        2.166
#> nodematch.Sex          0.548
#> gwesp.fixed.0.25       1.544
#> nodefactor.Grade.8     0.096
#> nodefactor.Grade.9     0.061
#> nodefactor.Grade.10    0.215
#> nodefactor.Grade.11    0.285
#> nodefactor.Grade.12    0.571

Core. nodematch("Grade") tests homophily, whether ties form between students of the same grade. nodefactor("Grade") tests activity, whether students in a given grade form more ties overall, regardless of who the partner is. One is about sorting, the other is about volume, and they are completely different claims.

Here the nodefactor coefficients are all positive relative to the reference category (grade 7) and rise more or less monotonically, from 0.1 for grade 8 up to 0.57 for grade 12. So older students in this school are more active, not just differently sorted, and note that grade homophily is still strongly positive (2.17) once activity is controlled for. Sorting and volume are distinct features of the network.

Stretch.

m4 <- fit_or_load("m4_race",
  ergm(mesa ~ edges +
         nodematch("Race", diff = TRUE,
                   levels = c("White", "Hisp", "NatAm")) +
         gwesp(0.25, fixed = TRUE),
       control = control.ergm(seed = 6886)))
data.frame(estimate = round(coef(m4), 3))
#>                      estimate
#> edges                  -5.743
#> nodematch.Race.White    1.090
#> nodematch.Race.Hisp     0.217
#> nodematch.Race.NatAm    0.635
#> gwesp.fixed.0.25        1.848
ame(m4, "nodematch.Race.White")
#>         AME       scale 
#> 0.007611043 0.006983515
c(baseline_density = network.density(mesa))
#> baseline_density 
#>      0.009708274

White–white homophily is the largest of the three (1.09, against 0.22 for Hispanic and 0.63 for Native American). Its AME is 0.76 percentage points, against a baseline density of 0.0097, so it shifts the tie probability by roughly 0.8× the baseline while remaining well under one percentage point in absolute terms. Report both, exactly as in B2.

Why be cautious about NatAm: the mixing matrix shows 46 within-NatAm ties out of a group of 68 students, which is plenty of data; that one is fine. The fragile estimate is White, with only 4 within-group ties among 18 students. A category can be large enough to avoid -Inf and still be small enough that the estimate is driven by a handful of dyads. “Not infinite” is not the same as “well identified.”

The full toolbox: model assessment beyond the core path

Three tools the core path walks past. Each closes a loop something above opened.

anova(): nested model comparison, done properly.

anova(m0, m1)
#> Analysis of Deviance Table
#> 
#> Model 1: mesa ~ edges
#> Model 2: mesa ~ edges + nodematch("Grade") + nodematch("Sex") + nodecov("Grade")
#>      Df Deviance Resid. Df Resid. Dev Pr(>|Chisq|)    
#> NULL                 20910    28987.4                 
#> 1     1  26701.7     20909     2285.7    < 2.2e-16 ***
#> 2     3    376.6     20906     1909.1    < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

An analysis-of-deviance table across nested fits. It is fully valid here because m0 and m1 are dyad-independent, so the likelihoods are exact. For MCMC fits the log-likelihood is bridge-sampled and carries its own standard error; the AIC caution in the depth box above applies to this table too, so do not lean on small deviance differences between dependent models. One place this is especially useful: the nodematch help page itself recommends an anova() of uniform (diff = FALSE) against differential (diff = TRUE) homophily.

predict(): the packaged version of B2’s hand arithmetic.

head(predict(m1), 3)
#>   tail head           p
#> 1    1    2 0.040656162
#> 2    1    3 0.001659101
#> 3    1    4 0.001270502

One row per dyad; p is $\mathrm{logit}^{-1}(\theta^\top \delta)$, the very numbers our three-line ame() built from ergmMPLE() in B2, now with a stable API around them. The default conditional = TRUE evaluates at the observed graph, exactly as our version did; conditional = FALSE simulates from the model instead, which is the equilibrium-distribution answer to the first caveat in B2’s depth box.

simulate() for one whole network: seeing is believing. The core path drew simulated statistics from m2; you can also ask for the graph itself and put it next to the data:

one <- simulate(m2, nsim = 1, seed = 6886, output = "network")
par(mfrow = c(1, 2), mar = c(0, 0, 1, 0))
set.seed(6886)
plot(mesa, vertex.col = "Grade", vertex.cex = 0.7, main = "Observed")
plot(one,  vertex.col = "Grade", vertex.cex = 0.7, main = "One draw from m2")

The draw has 203 edges to Mesa’s 203, the same sparse, grade-sorted, locally clumpy world, different in every particular.

gof() with a suite you chose. The default panel is a choice, not a law:

g_suite <- fit_or_load("tb_gof_suite",
  gof(m2, GOF = ~ degree + esp + distance,
      control = control.gof.ergm(nsim = 200, seed = 6886)))
plot(g_suite)

GOF = names the statistics; control.gof.ergm(nsim =) sets the number of simulated networks (the default is 100; doubling it steadies the boxplot whiskers). The rule worth keeping: if a statistic matters to your substantive argument, put it in the suite and report it. Do not let the default panel decide what “fits well” means for your paper.

# Questions and Notes

set.seed(6886)
gm <- gof(m2, GOF = ~model)
checkpoint(gof_model_p = round(gm$pval.model[, "MC p-value"], 2),
           label = "CHECKPOINT C2")
#> ------------------------------------------------------------------
#> CHECKPOINT C2: gof_model_p = 0.64, 0.68, 0.76, 0.64
#> ------------------------------------------------------------------

C3 · Optional Reference: Valued ERGMs

This section is not part of the central lecture route. It remains in the handout for students whose outcomes are counts or bounded values. The class moves from advanced diagnostics directly to Temporal ERGMs.

Everything so far has assumed ties are binary. Often they are not: trade volumes, message counts, how many times one person nominated another. This is the Krivitsky (2012) reading, and we can cover it directly without leaving the monastery.

Sampson observed his monks at three time points. ergm ships all three waves, so we can build a count: how many of the three waves did $i$ nominate $j$?

suppressMessages(library(ergm.count))
data(samplk)
tot <- as.matrix(samplk1) + as.matrix(samplk2) +
       as.matrix(samplk3)
table(tot)
#> tot
#>   0   1   2   3 
#> 236  38  20  30

So most pairs never nominate each other; 30 pairs did so in all three waves. That is a valued network with values 0–3.

nw <- as.network(tot, directed = TRUE, matrix.type = "adjacency",
                 ignore.eval = FALSE, names.eval = "nom")
nw
#>  Network attributes:
#>   vertices = 18 
#>   directed = TRUE 
#>   hyper = FALSE 
#>   loops = FALSE 
#>   multiple = FALSE 
#>   bipartite = FALSE 
#>   total edges= 88 
#>     missing edges= 0 
#>     non-missing edges= 88 
#> 
#>  Vertex attribute names: 
#>     vertex.names 
#> 
#>  Edge attribute names: 
#>     nom

The Idea You Have to Get: Sufficient Statistics Are Not Enough

With binary ties, writing down $g(y)$ fully specifies the model. With valued ties it does not. If a tie can be any non-negative integer, you have to say something about the shape of the distribution before the statistics do their work. Is a value of 5 twice as surprising as a 2? Is there an upper bound?

That “something” is the reference measure. It fixes the sample space and the baseline distribution against which $\theta^\top g(y)$ tilts. ergm.count implements Poisson, binomial, geometric and discrete-uniform. Poisson is often a sensible starting point for unbounded counts, but these nominations are bounded at three because we summed exactly three waves. A binomial reference with size three respects that design. I use Poisson here only to make the excess-zero problem especially transparent, not as the final reference measure I would defend for this outcome.

sum is the valued analogue of edges: the sum of all tie values, so its coefficient controls overall volume.

v1 <- fit_or_load("v1_sum",
  ergm(nw ~ sum, response = "nom", reference = ~Poisson,
       control = control.ergm(seed = 6886)))
round(summary(v1)$coefficients, 4)
#>     Estimate Std. Error MCMC % z value Pr(>|z|)
#> sum  -0.6023     0.0731      0 -8.2374        0

sum is -0.602, which is negative. Read naively: monks nominate each other less than a Poisson baseline would predict.

That reading is wrong, and here is how you find out.

v2 <- fit_or_load("v2_sumnonzero",
  ergm(nw ~ sum + nonzero, response = "nom", reference = ~Poisson,
       control = control.ergm(seed = 6886)))
round(summary(v2)$coefficients, 4)
#>         Estimate Std. Error MCMC % z value Pr(>|z|)
#> sum       0.3778     0.1014      0  3.7270    2e-04
#> nonzero  -2.1068     0.2248      0 -9.3716    0e+00

sum has changed sign: -0.602 → 0.378.

This is the whole Krivitsky lesson in two models. A bare Poisson reference cannot accommodate how many zeros the observed network has: 218 of the 306 ordered pairs here never nominate at all (the table() above prints 236 zeros, but 18 of those are the structural diagonal, since a monk cannot nominate himself). With no way to represent that, the only lever the model has is to push overall volume down, so sum goes negative. It is not measuring low volume; it is absorbing excess zeros.

Add nonzero, a term that separately controls how many ties exist at all, and the two jobs separate. nonzero is -2.11, strongly negative, capturing the sparsity. And now sum tells you the truth: conditional on a nomination existing at all, monks nominate each other more than Poisson predicts.

The general lesson transfers directly: in valued models, zero-inflation will masquerade as a volume effect unless you give the model a way to represent zeros separately. If you fit a valued ERGM with only sum, you have almost certainly mis-read your own coefficient.

The caveat ergm prints, which you should read

Every valued fit prints this, and it is easy to scroll past:

Null model likelihood calculation is not implemented for valued ERGMs at this time. This means that all likelihood-based inference (LRT, Analysis of Deviance, AIC, BIC, etc.) Is only valid between models with the same reference distribution and constraints.

So you may compare our two models by AIC; they use the same reference and constraints. You may not compare a Poisson-reference model to a geometric-reference model that way. That comparison appears in the literature and it is not valid.

Depth: what else is in there, and where to go next

ergm.count has a full vocabulary beyond sum and nonzero: nodesqrtcovar, nodecovar, transitiveweights() with choices of two-path, combining and affecting functions, and so on. transitiveweights in particular is the valued analogue of triadic closure, and, as in the binary case, the choice among its variants matters and is under-discussed.

I am deliberately not doing that taxonomy today. The core path gives you the reference-measure concept and the zero-inflation trap, which is what you need to read a valued ERGM paper without being fooled. The taxonomy needs an hour it does not have here, and it is in the ancillary document valued_ergm_depth.qmd.

One practical note: gof() gained support for valued ERGMs in ergm 4.11, so the diagnostics from C2 apply here too. That is recent; papers from before about 2023 mostly could not do it, which is why so few valued ERGM papers report goodness of fit.

checkpoint(sum_alone = coef(v1)[["sum"]],
           sum_with_nonzero = coef(v2)[["sum"]],
           nonzero = coef(v2)[["nonzero"]],
           label = "CHECKPOINT C3")
#> ------------------------------------------------------------------
#> CHECKPOINT C3: sum_alone = -0.602   |   sum_with_nonzero = 0.378   |   nonzero = -2.107
#> ------------------------------------------------------------------

C4 · Temporal ERGMs: The Main New Model for Today

A cross-sectional ERGM represents dependence within one observed graph. A new problem begins when we observe a sequence. An alliance at time $t$ may exist because it formed during the last interval, because it survived from an earlier period, or because the actors were never eligible to form it. Those are different events. Pooling the waves into one cross-section erases the distinction before estimation begins.

Suppose we observe the same network at several waves, $Y^1, Y^2, \ldots, Y^T$: friendships each school term, alliances each decade, or legislative cosponsorship networks in successive sessions. Under a first-order Markov specification, the distribution of the current network depends on the preceding network but not on earlier waves once $Y^{t-1}$ and the included covariates are known. The joint distribution then factors into transitions:

\[\Pr(Y^2, \ldots, Y^T \mid Y^1) = \prod_{t=2}^{T} \Pr(Y^t \mid Y^{t-1}), \qquad \Pr(Y^t \mid Y^{t-1}, \theta) = \frac{\exp\{\theta^{\top} g(Y^t,\, Y^{t-1})\}}{\kappa(\theta, Y^{t-1})}\]

The statistics $g(Y^t,Y^{t-1})$ can describe configurations in the current graph, dyad states in the previous graph, or a relationship between the two. This is the central TERGM move (Robins and Pattison 2001; Hanneke, Fu, and Xing 2010). Familiar ERGM tools still matter, but the estimand, risk set, assumptions, and validation strategy are now temporal.

The Temporal Estimand

State the outcome before you choose the terms

A cross-sectional ERGM models one network conditional on its covariates, constraints, and sample space. It does not condition on a preceding network. A TERGM models a sequence of network transitions, commonly by specifying $Y^t$ conditional on $Y^{t-1}$. Its coefficients describe conditional associations with the current tie or network state under that transition distribution. They are not automatically causal effects of actors making choices.

That distinction changes the political sentence we can write. A cross-sectional mutual term says reciprocated pairs are more common than the rest of the specification predicts. A lagged reciprocity term asks whether receiving a tie in the previous interval predicts sending one back in the current interval. The first describes a configuration. The second describes an ordered sequence.

The first-order Markov assumption is substantive, not merely computational. In an alliance application, it says that once last period’s alliance network and included covariates are known, alliances from still earlier periods add no further predictive information. A longer memory may be more plausible when treaties, wars, or institutions have effects lasting beyond one observation interval. If so, include longer lags or explain why the chosen lag captures the politically relevant history.

Audit the Panel Before You Write the Formula

The likelihood only makes sense after we decide which dyads could have changed in each interval. For every wave, check the following before fitting anything:

Actor roster and eligibility. Which actors exist and are eligible to send and receive ties? State dissolution, entry into an organization, or redistricting can change the risk set. Do not silently code an ineligible dyad as an observed non-tie.
Zero, missing, or impossible. An observed zero says a tie was possible and absent. Missing says its state is unknown. A structural zero says the tie was impossible by design. Treating all three as zero changes the estimand.
Stable labels and direction. Matrix row and column positions must identify the same actors at every wave. Confirm whether $i \to j$ means that $i$ initiates, nominates, sends, or receives the relation.
Covariate timing. A time-varying predictor must be aligned with the transition it is meant to explain. Decide whether $X^{t-1}$ predicts change during $(t-1,t]$, or whether contemporaneous $X^t$ is defensible. The second choice can create simultaneity.
Comparable observation intervals. A one-month transition and a ten-year transition do not expose ties to the same opportunity for change. Unequal spacing needs an explicit argument or a model that accommodates it.
Number of transitions. $T$ observed waves provide only $T-1$ transitions. The first wave conditions the first transition but is not itself a transition outcome. Ten thousand dyads observed at three waves are still only two realizations of the temporal process.

For political science, the fifth and sixth checks are especially easy to miss. A long panel with many states does not create many temporal replications if it contains only a handful of waves. The dyads contribute information within each transition, but uncertainty about a process evolving over time cannot be learned from network size alone.

Temporal Statistics Ask Sequence Questions

Three terms illustrate why we collected a panel in the first place:

memory(type = "stability") rewards the current dyad for matching its lagged state. The covariate is $+1$ when the lagged tie is present and $-1$ when it is absent. A previously tied dyad therefore receives $+\theta$ in its current-tie log-odds, while a previously untied dyad receives $-\theta$. The contrast is $2\theta$, not $\theta$. This term combines persistence of ties and persistence of non-ties, so it should not be described as tie survival alone.
delrecip() asks whether $j \to i$ at $t-1$ is associated with $i \to j$ at $t$. That is delayed reciprocity. A current-wave mutual term asks a different question: whether the two directions co-occur at $t$. One term does not substitute automatically for the other.
A lagged two-path covariate asks whether $i \to k \to j$ at $t-1$ is associated with $i \to j$ at $t$. This is a sequence-based closure hypothesis. Current-wave gwesp instead represents dependence among ties within $Y^t$.

The practical rule is simple: write the before-and-after configuration represented by every temporal statistic. If the verbal claim mentions sequence but the statistic uses only the current graph, the model and the theory have drifted apart.

Choose the Estimation Strategy Deliberately

There is no universally best TERGM estimator. The choice should follow the estimand, the number of waves, network size, and the dependence terms in the specification.

Approach	What it estimates	What it is useful for	What needs special care
`btergm()`	MPLE with bootstrap uncertainty	Fast estimation for network panels and a familiar `ergm` term vocabulary	The bootstrap resamples temporal observations. With very few transitions, it cannot create the temporal information the panel lacks.
`mtergm()`	MCMC maximum likelihood for a temporal model represented as a constrained block-diagonal network	Likelihood-based estimation when within-wave dependence makes MPLE unattractive	More computation, MCMC diagnostics, and the same need for an adequate number of waves
Statnet `tergm(..., estimate = "CMLE")`	Conditional maximum likelihood for a network series	Direct modeling of transitions with Statnet’s temporal framework, including separable formation and dissolution specifications	MCMC behavior, constraints, risk-set definition, and model degeneracy still require diagnosis

The worked example uses btergm() because it lets us see the sequence terms quickly. That is a teaching choice, not a claim that MPLE with a bootstrap is always preferable. Leifeld, Cranmer, and Desmarais (2018) discuss btergm() and mtergm() in detail; the Statnet TERGM documentation develops conditional and separable specifications.

For the worked btergm() fit, the point estimate maximizes a pooled pseudolikelihood over the observed transitions:

\[ \widetilde\ell(\theta) =\sum_{t=2}^{T}\sum_{i\ne j} \left[ y_{ijt}\log\widetilde p_{ijt}(\theta) +(1-y_{ijt})\log\{1-\widetilde p_{ijt}(\theta)\} \right], \]

where $\widetilde p_{ijt}(\theta)$ is the conditional probability of the current tie given the lagged network, the current states of the other dyads, and the declared change statistics. This resembles a pooled logistic regression, but dependence terms make it a pseudolikelihood rather than the full transition likelihood. The temporal bootstrap repeatedly resamples observed transitions and refits this objective to describe uncertainty. It cannot recover time periods that the design never observed.

The Data: 50 Adolescents Observed at Three Waves

The teenage friends and lifestyle panel ships inside RSiena: 50 Scottish adolescents, directed friendship nominations at three waves, with alcohol use and smoking measured at every wave. In this application the outcome is whether student $i$ nominates student $j$ as a friend in wave 2 or wave 3. We first let netify verify that actor labels, direction, and matrix dimensions agree across waves. btergm then needs a list of network objects, one per wave, so to_statnet() performs that handoff before we attach each wave’s covariates:

suppressMessages({ library(RSiena); library(btergm) })  # Loaded HERE, not in setup
tw    <- list(s501, s502, s503)                 # Three 50x50 adjacency matrices
friend_panel <- netify::netify(
  tw,
  symmetric = FALSE,
  diag_to_NA = TRUE,
  missing_to_zero = FALSE
)
stopifnot(all(unlist(netify::validate_netify(friend_panel))))
tnets <- netify::to_statnet(friend_panel)
for (i in seq_along(tnets)) {
  network::set.vertex.attribute(tnets[[i]], "drink", s50a[, i])
  network::set.vertex.attribute(tnets[[i]], "smoke", s50s[, i])
}
sapply(tnets, network.edgecount)
#>   1   2   3 
#> 113 116 122

Stuck?: two things in that chunk that are deliberate

The load order. btergm gets loaded here, not at the top of the day, and that is not sloppiness: the moment its namespace loads it overrides ergm’s plot and print methods for GOF objects. Run C2’s gof() output after loading btergm and the plots break with a baffling error. If that happens to you, restart R and run the day in order. Package conflicts of exactly this kind are endemic in the Statnet ecosystem; it is also why the network:: prefix in the loop is not decoration.

The loop. In base R a for loop cannot modify a list element in place because R copies on assignment. network objects break that rule: they carry an external reference, so set.vertex.attribute() modifies in place and the loop works as written.

Fit It Twice: Without the Past, Then With It

First fit the outcome waves without temporal terms, the model a cross-sectional habit would write. We use waves two and three here so that this model and the temporal model below explain the same two outcome networks:

set.seed(6886)
t_cross <- fit_or_load("tergm_s50_cross_outcome_waves",
  btergm(tnets[-1] ~ edges + mutual + gwesp(0.25, fixed = TRUE) +
                 nodecov("drink") + absdiff("drink") + nodematch("smoke"),
         R = 200, verbose = FALSE))
round(coef(t_cross), 3)
#>                edges               mutual gwesp.OTP.fixed.0.25 
#>               -4.645                3.601                1.131 
#>        nodecov.drink        absdiff.drink      nodematch.smoke 
#>                0.043               -0.248                0.258

(Note fixed = TRUE; C1’s trap follows you here because btergm calls ergm’s term machinery.) Now hand the model last wave:

set.seed(6886)
t_stab <- fit_or_load("tergm_s50_stab",
  btergm(tnets ~ edges + mutual + gwesp(0.25, fixed = TRUE) +
                 nodecov("drink") + absdiff("drink") + nodematch("smoke") +
                 memory(type = "stability"),
         R = 200, verbose = FALSE))
round(coef(t_stab), 3)
#>                edges               mutual gwesp.OTP.fixed.0.25 
#>               -3.872                3.075                1.011 
#>        nodecov.drink        absdiff.drink      nodematch.smoke 
#>                0.048               -0.177                0.235 
#>  edgecov.memory[[i]] 
#>                1.049

Watch what one temporal term does to the structural story. mutual moves from 3.6 to 3.08, and gwesp from 1.13 to 1.01, the moment the model learns what the network looked like last wave, while memory itself comes in at 1.05 and is, as promised, doing enormous work. The movement is a warning that the cross-sectional structural story was partly entangled with dyad stability. It is not, by itself, an estimate of how much bias memory removed. The two models define different conditional distributions, and their log-odds scales are not directly comparable. The defensible next step is to compare recognizable configurations on a probability scale and then assess simulated transitions.

Two Transitions Are Not an Inferential Foundation

Point estimates are not inference, so before anyone quotes these numbers, where are the standard errors? There are none, by design; btergm’s uncertainty comes from the bootstrap, and you ask for it with confint():

round(confint(t_stab), 2)
#>                      Estimate Boot mean  2.5% 97.5%
#> edges                   -3.87     -3.79 -4.10 -3.41
#> mutual                   3.08      3.11  2.64  3.55
#> gwesp.OTP.fixed.0.25     1.01      1.02  0.96  1.08
#> nodecov.drink            0.05      0.04 -0.04  0.10
#> absdiff.drink           -0.18     -0.18 -0.20 -0.17
#> nodematch.smoke          0.24      0.23  0.07  0.43
#> edgecov.memory[[i]]      1.05      1.05  1.05  1.07
cat("unique coefficient vectors across the 200 bootstrap draws:",
    nrow(unique(round(t_stab@boot$t, 12))), "\n")
#> unique coefficient vectors across the 200 bootstrap draws: 3

Those are percentile intervals from the R = 200 pseudolikelihood bootstrap we already paid for at fit time. But this tiny panel has only two usable transitions. Resampling temporal observations produces just 3 distinct coefficient vectors across all 200 draws, so the apparently precise percentile endpoints are mostly repeated copies of three fits. Use this table to learn the software, not to make a publishable uncertainty claim. A defensible application usually needs many more waves, or a different inferential design justified for the panel at hand. A larger network does not solve this temporal replication problem. (Small quirk to expect: summary() on a btergm object prints its table twice in the current release; this is cosmetic, not two models.)

Do not turn coefficient movement into a causal result

The coefficient movement tells us that the specification matters. It does not identify which part of the movement is rescaling, changed dependence, or omitted temporal structure. For a comparison on the friendship-probability scale, use btergm::interpret() to calculate conditional tie probabilities or first differences for stated configurations (Leifeld, Cranmer, and Desmarais 2018). Report exactly which lagged and current network features are being held fixed.

Formation and Persistence Have Different Risk Sets

The stability term is useful, but it compresses two processes into one statistic: existing ties remaining present and absent ties remaining absent. Political theories often distinguish them. Common threats may help states form an alliance, while institutionalization may help an existing alliance endure. A single coefficient need not describe both.

A separable temporal ERGM, or STERGM, writes separate submodels for two risk sets:

Process	Dyads at risk	Political question
Formation	Ties absent at $t-1$	Which eligible non-allied pairs form an alliance by $t$?
Persistence	Ties present at $t-1$	Which existing alliances survive through $t$?

In Statnet syntax, Form() places terms in the formation equation, while Persist() places terms in the persistence equation. A positive Persist() coefficient means less dissolution. Diss() uses the same post-dissolution network but negates the statistics, so a positive Diss() coefficient means more dissolution. The software handles the transition networks and constraints, but the researcher still has to name the risk set and interpret the sign under the selected parameterization.

The same predictor can therefore have one association with formation and another with persistence. Regime similarity might predict which dyads begin cooperating but tell us little about whether established cooperation survives a leadership change. That is a more precise claim than saying the predictor is simply associated with “the alliance network.”

Separability is also an assumption. It treats formation and dissolution as conditionally independent components of a transition given the preceding network and model terms. It does not automatically make either equation causal, and it does not recover unobserved changes between widely spaced panel waves.

Temporal Fit Must Be Temporal

A TERGM should not be judged only by whether it reproduces each wave’s degree distribution. It also needs to reproduce the kinds of change on which the substantive claim rests. Simulate one-step-ahead networks conditional on the observed preceding wave, then compare observed and simulated quantities such as:

numbers of newly formed, persistent, and dissolved ties;
reciprocity, degree, shared partners, and mixing in the next wave;
tie duration or stability when those are part of the theory;
delayed reciprocity and lagged closure counts;
calibration of predicted tie or transition probabilities.

When the panel is long enough, hold out later waves and forecast them from earlier ones. Keep one-step-ahead validation separate from recursive forecasting, where a simulated network is fed forward as the next lag. The second task compounds model error and answers a harder question.

Sparse networks make ordinary classification accuracy nearly useless because a model can predict “no tie” for almost every dyad and look accurate. Use calibration, precision-recall, transition-specific checks, and substantively important network statistics instead. With only two transitions, the classroom data can demonstrate these ideas but cannot provide a credible held-out temporal test.

Tiered Exercise

Read the memory coefficient in t_stab off the output above. In one sentence: what does it say about these friendships, and why is it supposed to be large? (Hint: what would the waves be, if it were zero?)

Add delrecip() to t_stab and refit (seed 6886, R = 200). What happens to mutual, and what does delrecip measure that the cross-sectional mutual in B3 cannot, even in principle?

set.seed(6886)
t_temp <- fit_or_load("tergm_s50_temp",
  btergm(tnets ~ edges + mutual + gwesp(0.25, fixed = TRUE) +
                 nodecov("drink") + absdiff("drink") + nodematch("smoke") +
                 memory(type = "stability") + delrecip(),
         R = 200, verbose = FALSE))
data.frame(term          = c("mutual", "gwesp", "memory", "delrecip"),
           stability_only = round(coef(t_stab)[c("mutual", "gwesp.OTP.fixed.0.25",
                                                 "edgecov.memory[[i]]", NA)], 3),
           plus_delrecip  = round(coef(t_temp)[c("mutual", "gwesp.OTP.fixed.0.25",
                                                 "edgecov.memory[[i]]",
                                                 "edgecov.delrecip[[i]]")], 3),
           row.names = NULL)
#>       term stability_only plus_delrecip
#> 1   mutual          3.075         2.955
#> 2    gwesp          1.011         1.001
#> 3   memory          1.049         0.911
#> 4 delrecip             NA         0.528
round(confint(t_temp)["edgecov.delrecip[[i]]", ], 3)
#>  Estimate Boot mean      2.5%     97.5% 
#>     0.528     0.512     0.412     0.600
cat("unique coefficient vectors across the 200 bootstrap draws:",
    nrow(unique(round(t_temp@boot$t, 12))), "\n")
#> unique coefficient vectors across the 200 bootstrap draws: 3

Core. Dyads overwhelmingly stay in whatever state they were in last wave: ties persist, empty dyads stay empty. It is supposed to be big: its entire job is to absorb persistence so that persistence stops leaking into reciprocity and closure. If it were zero, the waves would be independent redraws of the network and the pooled fit t_cross would have been fine.

Stretch. mutual barely moves again (3.08 → 2.95); the persistence was already absorbed by memory. What delrecip adds is not a control but a new sequence hypothesis: a tie extended to you last wave is associated with a 0.53 increase in the conditional log-odds of your returning it this wave. The point estimate is positive, but the same three-vector bootstrap problem means we do not treat its interval as credible evidence that the effect excludes zero. The statistic asks which direction came first, a question no cross-sectional model of any kind, ERGM or otherwise, can make it. Cross-sectional mutual only knows that reciprocated pairs are common; it cannot see the sequence.

checkpoint(mutual_pooled   = coef(t_cross)[["mutual"]],
           mutual_temporal = coef(t_stab)[["mutual"]],
           memory          = coef(t_stab)[["edgecov.memory[[i]]"]],
           label = "CHECKPOINT C4")
#> ------------------------------------------------------------------
#> CHECKPOINT C4: mutual_pooled = 3.601   |   mutual_temporal = 3.075   |   memory = 1.049
#> ------------------------------------------------------------------

One Sentence for Tomorrow, and the Deep Dive

The TERGM is tie-oriented: the unit that changes between waves is the dyad. Tomorrow’s model, the stochastic actor-oriented model, is actor-oriented: the unit that changes is an actor’s choice. Whether those two framings are rivals or complements is a genuine, published dispute, and it is literally where Day 12 opens. You now have a fitted TERGM of your own; bring it to the argument.

Everything past this, including goodness of fit for TERGMs (with out-of-sample tie prediction), the other memory types (autoregression, loss, innovation; “stability” is one of four), timecov() for secular trends, formation versus dissolution, actors entering and leaving the panel, a larger IR application, and the full dispute reading list, is in the deep-dive companion ancillary/tergm.qmd, which starts exactly where this section stops.

Wrap: ERGMs and Latent Variable Models

I want to be careful here, because it would be easy for me to be unfair.

On Day 10 you fitted latent factor models, including the AME framework from my own work. Today you fitted ERGMs. It would be convenient for me to tell you one beats the other. That is not true and the literature does not say it.

Where they agree, which is most of the territory. Dyadic observations are not independent. Standard logit is misspecified for network data. Omitted dependence biases the coefficients, not merely the standard errors, and no robust/clustered SE repairs coefficient bias from a misspecified mean. That is distinct from the robust-variance use case: when the mean is modeled adequately and dependence is treated purely as a nuisance for inference, dyadic cluster-robust SEs (the DCR approach, which gets a full session on Day 13) are a legitimate and useful tool, not settled to be inconsistent. The stronger claim, that robust clustered SEs are inconsistent for grouped discrete-choice data, and that complex network dependence prevents the node-wise partitioning clustering requires, is Cranmer & Desmarais’s (citing Greene), and it is contested by the robust-inference camp (Aronow, Samii & Assenova; Graham). Both literatures are allies against the same enemy, and that enemy is the dyadic logit with robust standard errors that still dominates several fields. Both also assess models the same way: simulate from the fit, compare to the observed data.

Where they differ: four axes, no winner.

	ERGM	Latent factor / AME
What dependence is	A hypothesis to parameterize and test	Residual structure to represent so $\beta$ can be examined under that specification
Direction of specification	Specify the joint; conditionals follow	Specify a conditional mean; dyads independent given the latents
Nature of the structure	Deterministic function of the observed network	Unobserved random quantities
How it fails	Degeneracy; MCMC-MLE non-convergence; decay non-identification / projectivity	Choosing $R$, label switching, MCMC mixing

The subtle one, stated plainly: AME does not abolish conditional independence. It conditions on more. Dyads are still conditionally independent given $u_i^\top D v_j$.

What neither can do from one cross-section. Separate triadic closure from homophily on an unobserved attribute. That is Shalizi & Thomas (2011), which you met on Day 10. An ERGM parameterizes closure and tests it; AME absorbs it into the latent terms. Neither identifies the other, and no amount of cleverness with a single cross-sectional network will fix that. That is the actual state of the literature, and it is exactly the setup for Day 12, tomorrow. The stochastic actor-oriented model uses panel data to parameterize temporal selection and influence separately. That makes the question sharper, but it does not identify either process automatically, especially in a small panel.

One connection I find especially useful. Box-Steffensmeier, Christenson & Morgan (2018, Political Analysis) propose the frailty ERGM, which adds node-specific random intercepts, or “sociality effects,” drawn from a normal centered at zero and extended to separate sender and receiver effects for directed networks, because unmeasured node attributes otherwise cause “coefficient bias and model degeneracy.”

That is the social relations model, inside an ERGM.

Here is the thing that should bother you. On Day 9 you fit sender and receiver random effects and called it the SRM. The ERGM literature, a separate community, hit degeneracy, went looking for a fix, and landed on node-specific sender and receiver random intercepts. They rediscovered your Day 9 model from the other side of the room.

The two traditions are converging on the same structure from opposite directions.

Depth: what I am and am not claiming

Three things I want on the record, because they are easy to get wrong and colleagues will catch you:

I am not claiming “AME beats ERGM.” “Taking Dyads Seriously” mentions ERGMs in a single footnote and estimates exactly zero of them. Its comparison is against the dyadic GLM with robust standard errors, which is what all three of its replication targets used. The paper’s opponent is the naive logit, not the ERGM.

Cranmer & Desmarais are not attacking latent space models. Their sharp passage about “applying standard regression” targets the practice of computing network statistics and plugging them into a logit as covariates. The paragraph immediately after is explicitly irenic: “it is not our intention here to pit the ERGM against these alternative approaches… The quest for a universally best method is quixotic.”

And their own stated limits are worth noting, because today quietly answered them: they name ERGM’s chief limitations as binary-only and cross-sectional-only. The first was solved by Krivitsky, which you did in C3. The second by temporal ERGMs, which you fit yourself in C4. The reading list has an internal arc, and today walked all of it.

For the strongest constructive critique from inside the ERGM camp, look at LOLOG (Fellows & Handcock; Clark & Handcock 2022, JRSS-A), which models network growth under a latent node ordering and claims to avoid degeneracy. Useful fact: LOLOG and ERGM are identical under dyadic independence, so everything in blocks A and B transfers directly, and the divergence is exactly at the dependence terms. That this comes from Handcock, an ERGM founder, makes it impossible to frame as camp warfare.

The referee checklist for applied researchers

Researchers often read more ERGMs than they fit. So here is the checklist today built, and it is the deliverable:

Is dependence modelled or absorbed? And does the paper say which?
Which structural terms, and why? Justified substantively, or reverse-engineered from fit?
Was fixed = TRUE used on the geometrically weighted terms? If the decay was typed but not fixed, the reported model is not the model they think they fitted.
Is goodness of fit reported, and on which statistics? In-model only is not enough.
Are coefficients compared across nested specifications? If so, that comparison is invalid (Duxbury), and it is extremely common.
Are effects reported on a scale a reader can use, with the baseline density alongside?

After the Model Runs: Translate a Graph Parameter Into the Observed Action

An ERGM coefficient is not a standalone probability and should not be read like an ordinary logit coefficient when its change statistic varies across dyads. The clean route from the fitted model to an applied result has four steps:

Name the comparison. For a proposed tie toggle, state which network statistic changes and by how much, conditional on the rest of the graph and the other terms.
Put the result on a usable scale. Show conditional tie probabilities, first differences, or average marginal effects for recognizable configurations. A raw log-odds coefficient is rarely the final substantive result.
Simulate from the fitted model. Check degree, shared-partner, geodesic, mixing, and other statistics that matter for the argument, especially statistics not included directly in the model.
Limit the mechanism claim. A positive closure term means the observed network is more compatible with ties that close shared partners, conditional on the specification. It does not prove that actors consciously seek triangle closure.

What the Mesa result means

“Fewer than 1 in 100 possible student pairs in the Mesa network are friends. After accounting for that sparse baseline, two students being in the same grade raises their fitted friendship probability by about 2.8 percentage points on average. That change is small in absolute percentage points but large relative to the network’s baseline density of 0.97 percent. In the richer model, friendship toggles that close shared partners also receive greater conditional weight, and the simulation checks tell us whether that specification reproduces the observed distribution of shared friends. These are conditional descriptions of this school network, not evidence that students consciously count triangles.”

For the temporal comparison, add one more sentence: once prior ties enter the model, say which apparent structural relationships survive and which were largely describing persistence. With only three waves and two transitions in this demonstration, the correct result is a worked comparison, not a broad claim about a stable formation process.

What we deliberately did not show: five tools, five pointers

Today demonstrated everything you need to build, fit, diagnose, interpret and extend an ERGM. Five pieces of the machinery were left out on purpose, so that you know they exist and where to look rather than wondering whether I forgot them:

san(): simulated annealing constructs a network whose statistics hit target values. Used internally for MCMC starting points; useful to you for building synthetic test networks. ?san.
ergm.godfather(): apply a scripted sequence of tie toggles to a network and watch the statistics change: the hand-crank version of the change-statistic machinery from A3. ?ergm.godfather.
ergm.exact() / ergm.allstats(): the exact likelihood by brute-force enumeration of all $2^{\binom{n}{2}}$ graphs. Only feasible for networks of roughly eight nodes or fewer, which is exactly what makes it the right tool for checking your intuitions about A2’s normalizing constant. ?ergm.exact.
offset() / Offset(): fix a coefficient at a known value instead of estimating it: the device behind size-adjusted parameterizations (Krivitsky, Handcock & Morris 2011) and hard structural zeros. ?ergmTerm under “Offset”; the operator version is in the appendix-A list.
obs.constraints = and missing data: ergm() treats NA edges as unobserved rather than absent, and can additionally model the observation process itself. If your network has unobserved dyads, this is the principled route, not imputing zeros. ?ergmConstraint, and the observational-process section of Krivitsky, Hunter, Morris & Klumb (2023, JSS).

The other two stated omissions of the day have homes already: the full valued-ERGM taxonomy is in ancillary/valued_ergm_depth.qmd (C3), and everything temporal beyond C4’s two terms is in ancillary/tergm.qmd.

Where We Go Next: Day 12

Absorb it (Day 10) → test it (Day 11) → parameterize temporal mechanisms (Day 12). Time makes the hypotheses more distinguishable, but it does not identify selection and influence automatically.

Everything until C4 conditioned on one cross-section, and we hit a hard wall: you cannot separate closure from homophily on unobservables. C4 made the first move that waves allow: condition on the preceding wave in a tie-oriented model. Day 12 uses the same kind of panel data in a continuous-time, actor-oriented model to represent selection and influence separately and ask what the panel can actually distinguish. It opens by asking whether the model fitted here is the right way to use waves at all. Day 13 then separates dependence adjustment from causal identification, and Day 14 compares the model families on one dataset.

Appendix A: Depth Track

What changed in ergm 4, for anyone reading older tutorials

The assigned Hunter et al. (2008) JSS paper describes ergm 2. We are on 4.12.0. The things that will actually trip you up:

Curved terms ignore decay unless fixed = TRUE (4.0). Demonstrated in C1. This is the one that silently changes results.
levels2 = default flipped in nodemix and mm; it is now -1, excluding the first level combination. Reference cells in older output differ.
base = and keep = are gone → use levels = and nodes =.
fit$coef → fit$coefficients. Use coef(fit) and stop worrying.
control.ergm(SAN.control =, loglik.control =) → SAN =, loglik =.
Term operators are new in 4.0: B(), F()ilter, Offset(), Symmetrize(), Sum(), S()ubgraph, Log(), Exp(), and dyad-independent interactions via TERM1:TERM2. A pre-4 lecture teaches specification as a fixed vocabulary; ergm 4 makes it compositional. This is the largest new capability and it is barely taught anywhere; the B3 toolbox runs one (F() rebuilding nodematch from parts) so you have a working template.
snctrl() is the one control shorthand to remember: control = snctrl(seed = 6886) works across ergm(), simulate(), gof(), san() and friends, dispatching to the right control.*() constructor. Demonstrated in the B3 toolbox.
New defaults: MCMLE.termination = "confidence", adaptive MCMC.samplesize, the .triadic proposal hint on by default (models with closure mix much better than they used to), and MCMLE.density.guard, the tripwire that saved us twenty minutes in C1 and which simply did not exist before.
4.11: gof() supports valued ERGMs; samplk1/2/3 gained edge weights, which is why C3 works on data you already knew.

The current companion paper is Krivitsky, Hunter, Morris & Klumb (2023), JSS 105(6), “ergm 4: New Features.” Read that alongside the 2008 paper, not instead of it; the 2008 GOF section is still the best exposition of the logic.

Scaling up: what to do above ~1,000 nodes

Plain MCMC-MLE has a practical ceiling around one to a few thousand nodes. Above that there are three practical routes:

MPLE with bootstrap or Godambe standard errors: built into ergm since 4.5; see B1. The pragmatic default.
bigergm (CRAN): hierarchical ERGM, successor to hergm/lighthergm. It recovers blocks first, then estimates within- and between-block effects, permitting dependence only within blocks (local dependence in the Schweinberger–Handcock sense).
ergm.multi (CRAN, installed here), not one big network but many networks. This is also the practical route to estimating decay properly, per the tension in C1.

I have not run bigergm for you, so treat that row as documentation rather than a recommendation. I mention all three because many applied networks are far larger than Mesa, and it is important to know that the same ergm() workflow will not simply scale without adjustment.

A netify handoff, if your data does not start as a network object

Most of you will not receive data as a network object; you will have an edgelist and a node table. netify gets you from one to the other without hand-writing conversion code:

suppressMessages(library(netify))
data(classroom_edges); data(classroom_nodes)
cls <- netify(classroom_edges, actor1 = "from", actor2 = "to",
              symmetric = TRUE)
cls <- add_node_vars(cls, classroom_nodes, actor = "student")
cls_net <- netify_to_statnet(cls)
summary(cls_net ~ edges + nodematch("gender") +
          nodematch("grade"))
#>            edges nodematch.gender  nodematch.grade 
#>               51               34               30

30 students, 51 friendships. That is too small to fit a credible ERGM; 51 edges will not support a gwesp term, so I am not modelling it. It is here purely to show the bridge from a netify object to something ergm() accepts. Day 8’s netify session covered the rest.

Appendix B: Going Further on Your Own

These are optional and each one runs end to end. All but the first cover material the 2026 course does not teach; the first is the deep-dive companion to a section that now is taught.

Document	What it solves	Extends
`ancillary/tergm.qmd`	The C4 deep dive: TERGM GOF and tie prediction, formation versus dissolution, composition change, the TERGM-versus-SAOM dispute in full	Today, C4
`ancillary/valued_ergm_depth.qmd`	The full valued vocabulary: `transitiveweights`, `nodecovar`	Today, C3
`ancillary/ergm_applications.qmd`	Two worked applications end to end	Today, block B
`ancillary/ego_ergms.qmd`	Mixtures of ERGMs over subgroups	Today + week 1 communities
`ancillary/random_graphs.qmd`	Erdős–rényi and the null-model logic from A1	Today, A1

The referee exercise (Duque 2018, assigned today). Read it, then run the six-question checklist from the wrap against it. This is the assignment for those of you who will be reading these models rather than fitting them, which is most of you. I am deliberately not having you refit it; it is 158 actors and the substantive setting takes half an hour to explain to a non-IR audience. Reading it critically is the transferable skill.

Appendix C: Session Info

sessionInfo()
#> R version 4.3.3 (2024-02-29)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 24.04.3 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> time zone: America/New_York
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] netify_1.5.3   btergm_1.11.1  RSiena_1.5.0   ergm_4.12.0    network_1.20.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] generics_0.1.4           slam_0.1-55              robustbase_0.99-7       
#>  [4] stringi_1.8.7            lattice_0.22-5           digest_0.6.39           
#>  [7] magrittr_2.0.5           statnet.common_4.13.0    evaluate_1.0.5          
#> [10] lpSolveAPI_5.5.2.0-17.15 grid_4.3.3               RColorBrewer_1.1-3      
#> [13] fastmap_1.2.0            jsonlite_2.0.0           Matrix_1.6-5            
#> [16] ggrepel_0.9.6            ggnewscale_0.5.2         backports_1.5.0         
#> [19] rle_0.10.0               purrr_1.2.2              scales_1.4.0            
#> [22] Rdpack_2.6.4             cli_3.6.6                rlang_1.2.0             
#> [25] rbibutils_2.4            cachem_1.1.0             yaml_2.3.12             
#> [28] otel_0.2.0               tools_4.3.3              parallel_4.3.3          
#> [31] checkmate_2.3.4          memoise_2.0.1            coda_0.19-4.1           
#> [34] dplyr_1.2.1              ergm.count_4.1.3         ggplot2_4.0.3           
#> [37] boot_1.3-30              ROCR_1.0-12              vctrs_0.7.3             
#> [40] R6_2.6.1                 lifecycle_1.0.5          stringr_1.6.0           
#> [43] htmlwidgets_1.6.4        MASS_7.3-60.0.1          trust_0.1-9             
#> [46] pkgconfig_2.0.3          sna_2.8                  pillar_1.11.1           
#> [49] gtable_0.3.6             Rcpp_1.1.1-1.1           glue_1.8.1              
#> [52] tidyselect_1.2.1         DEoptimR_1.1-4           xfun_0.55               
#> [55] tibble_3.3.1             knitr_1.51               farver_2.1.2            
#> [58] htmltools_0.5.9          igraph_2.2.2             rmarkdown_2.30          
#> [61] compiler_4.3.3           S7_0.2.2

Process	Dyads at risk	Political question
Formation	Ties absent at \(t-1\)	Which eligible non-allied pairs form an alliance by \(t\)?
Persistence	Ties present at \(t-1\)	Which existing alliances survive through \(t\)?

	ERGM	Latent factor / AME
What dependence is	A hypothesis to parameterize and test	Residual structure to represent so \(\beta\) can be examined under that specification
Direction of specification	Specify the joint; conditionals follow	Specify a conditional mean; dyads independent given the latents
Nature of the structure	Deterministic function of the observed network	Unobserved random quantities
How it fails	Degeneracy; MCMC-MLE non-convergence; decay non-identification / projectivity	Choosing \(R\), label switching, MCMC mixing

Before You Reproduce the Walkthrough

Start Here: What We Are Refreshing and Building On

A Fast Refresher: Question, Contrast, Estimation, and Interpretation

1. What Is an ERGM Meant to Answer?

2. How Is This Different From What We Just Did?

3. How Is an ERGM Estimated?

4. Why Does This Matter for Applied Research?

The Applied Model Compass

The Advanced Route for Today

ERGM Reference Bridge

A1 · Where We Left Off, and What Today Is For

The Substantive Version of the Question

Two Nulls, Two Different Stories

A2 · The Joint Model, and the Constant That Ruins Your Life

\(g(y)\) Is Just Counting, and You Can Do It Right Now

And Now the Constant

A3 · Change Statistics: The Idea the Whole Day Rests On

Let Us Actually Compute One

Advanced ERGM Work: Interpretation and Specification

How ERGM Estimation Works

B1 · Reference: The Dyad-Independent Model and a Confusing Console Message

The Thing That Confuses Everyone

B2 · What the Numbers Mean for Students and Friendships

Step One: The Baseline

Step Two: Why You Cannot Just Exponentiate and Stop

Step Three: Average Marginal Effects

Tiered Exercise

B3 · Endogenous Terms: The Reason You Are Here

Reciprocity

The Term Vocabulary Is Searchable: Do Not Memorize It

Transitivity, and the Naive Way to Ask for It

Advanced Failure Diagnosis and Network Change

C1 · Diagnosing Degeneracy

Why Triangles Do This

The Repair: Geometric Weighting

C2 · Diagnostics: Separation First, Then Mixing, Then Fit

1. Look at Your Cells Before You Fit

2. Did the Chain Actually Mix?

3. Goodness of Fit: Two Different Questions

4. And the Tool Underneath All of It: simulate()

Tiered Exercise

C3 · Optional Reference: Valued ERGMs

The Idea You Have to Get: Sufficient Statistics Are Not Enough

C4 · Temporal ERGMs: The Main New Model for Today

The Temporal Estimand

Audit the Panel Before You Write the Formula

Temporal Statistics Ask Sequence Questions

Choose the Estimation Strategy Deliberately

The Data: 50 Adolescents Observed at Three Waves

Fit It Twice: Without the Past, Then With It

Two Transitions Are Not an Inferential Foundation

Formation and Persistence Have Different Risk Sets

Temporal Fit Must Be Temporal

Tiered Exercise

One Sentence for Tomorrow, and the Deep Dive

Wrap: ERGMs and Latent Variable Models

After the Model Runs: Translate a Graph Parameter Into the Observed Action

Where We Go Next: Day 12

Appendix A: Depth Track

Appendix B: Going Further on Your Own

Appendix C: Session Info

4. And the Tool Underneath All of It: `simulate()`