Temporal ERGMs

Ancillary self-study · ICPSR Network Analysis: Advanced Topics

Author

Shahryar Minhas

Published

July 22, 2026

Optional: not covered in the 2026 course

This document is not part of any taught session. Nobody is going to lecture it at you, nobody is going to grade you on it, and you can ignore it completely without falling behind.

It exists because it is a thing you will plausibly need next, and because it is easier to hand you something that works than to hand you a citation and wish you luck. Everything here runs end to end on the same machine setup you used in class. Every number and every figure on this page was produced by a chunk in this file: there are no screenshots of R output anywhere in the ancillary package.

How to work through it alone. Get the source (click the </> Code button at the top right, then View Source, and copy it into a .qmd file: or open the .qmd you were given directly), run the chunks in order, and read the prose between them: the prose is where the argument lives. Slow fits are identified in the prose, and precomputed fits load a shipped result while still showing you the code that produced it. Folded Depth boxes are optional extras; folded Stuck? boxes are hints. Open whichever you need and ignore the other. Exercises have a Solution tab: you are working alone, so use it.

What this extends, and why it is here

Extends: the ERGM session (day 11, W3 Mon), and specifically its temporal block (C4). Everything you learned there about change statistics, specification, and goodness of fit carries over unchanged; C4 introduced time, and this document gives it the full treatment.

Where this sits in 2026: this is the day-11 deep dive. The taught session’s final block (C4) covers the TERGM core: the Markov factorization, memory and delrecip, and the s50 worked example, which day 11 promoted from this very document. What the taught block could not hold is here: the other memory flavours, GOF for temporal fits, the distinction between default in-sample checks and genuine out-of-sample prediction, the full TERGM-versus-SAOM comparison on the same data, an applied case, and the complete Leifeld–Cranmer / Block exchange. TERGM is also the single most likely thing you will actually reach for once you leave. Panel data on networks is common, and a cross-sectional ERGM on the last wave throws away most of what you paid to collect, so it gets the longest treatment on the shelf.

What you need: btergm, network, networkDynamic, tsna, RSiena (for the data), and networkdata (for the diplomacy example in §7 only). No install_github is executed anywhere in this file.

The examples are deliberately kept manageable. A self-study document that makes you wait for a fit you cannot interpret yet is a document you will close.

1. The problem, stated plainly

You have a network measured at several points in time. Friendships in a classroom in September, December, and April. Collaboration between organisations across five annual reports. Who cites whom, by year. Alliances between states, by decade.

You could fit an ERGM to the last wave and report it. People do this. The trouble is that almost every interesting question you have is about change, and a cross-section cannot see change. “Do people reciprocate ties?” is a different question from “do people reciprocate ties that were extended to them last year?” Only the second one is a claim about a process.

There is also a statistical problem hiding behind the substantive one. If you stack your waves and fit one ERGM to the pooled data, you are pretending that a tie present in wave 1 and still present in wave 2 is two independent observations. It is not. Ties are sticky. That stickiness is enormous, usually the single largest effect in the data, and if you do not model it, it leaks into everything else. You will find “reciprocity” that is really just persistence, and “triadic closure” that is really just the fact that last year’s triangles are still standing.

The temporal ERGM is the smallest sensible fix. Here is the whole idea.

1.1 The model

Start from where the ERGM session left off. For a single network,

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

Now suppose you have \(Y^1, Y^2, \ldots, Y^T\). Make a Markov assumption across waves: given \(Y^{t-1}\), the network at \(t\) is independent of everything before \(t-1\). Then the joint distribution of the whole sequence factors:

\[\Pr(Y^2, Y^3, \ldots, Y^T \mid Y^1) = \prod_{t=2}^{T} \Pr(Y^t \mid Y^{t-1}).\]

and each of those conditional pieces gets an ERGM:

\[\Pr(Y^t \mid \theta, Y^{t-1}) = \frac{\exp\{\theta^{\top} g(Y^t,\, Y^{t-1})\}}{\kappa(\theta, Y^{t-1})}.\]

That is it. That is the TERGM. The only change from an ordinary ERGM is that the sufficient statistics \(g(\cdot)\) are now allowed to depend on the previous network as well as the current one. Everything else, including change statistics, conditional log-odds, degeneracy, and the whole apparatus, is identical because it is an ERGM. The temporal part lives entirely in the statistics you are allowed to write down.

Depth: why this is really just a block-diagonal ERGM

If you take your \(T-1\) transition networks and stack them into one big block-diagonal adjacency matrix, with a constraint forbidding any edge in the off-diagonal blocks, then fitting a single ERGM to that big matrix is exactly fitting the TERGM. Each diagonal block is one \(\Pr(Y^t \mid Y^{t-1})\) term; the constraint enforces that no statistic ever spans two time points.

This is not a metaphor, it is the implementation. It is also why every ERGM warning you learned applies here verbatim: put triangle in a TERGM and you get the same degeneracy you got in the cross-section, just \(T-1\) times over.

The Markov-across-waves assumption is doing real work and is worth being uneasy about. It says that once I know last year’s network, the year before adds nothing. For annual data on slow processes that is often defensible. For three waves of adolescent friendship measured a term apart it is more of a hope than a claim. Robins & Pattison (2001) introduced the idea; Hanneke, Fu & Xing (2010) gave it the form above.

1.2 What you get for it: the two terms that matter

Two temporal statistics carry most of the value, and both are in btergm:

memory(type = "stability") is the previous network as a dyadic covariate. It asks: given what the dyad looked like last wave, how likely is it to look that way again? This is the stickiness term. It will almost always be large and it is supposed to be. Its job is to soak up persistence so that persistence stops masquerading as something more interesting. (Other flavours: "autoregression" for ties carried over, "loss" for dissolution, "innovation" for new ties among previously unconnected pairs.)
delrecip() captures delayed reciprocity. If \(j \to i\) existed last wave, is \(i \to j\) more likely this wave? This is a genuinely temporal question that a cross-section cannot ask at all. Cross-sectional mutual tells you that reciprocated pairs are common; it cannot tell you which direction came first.

There is also timecov(), which tests whether an effect strengthens or weakens over time, and which you should reach for only when you have enough waves to see a trend: five or six, not three.

2. Getting panel network data into shape

We will use the Teenage Friends and Lifestyle Study that ships with RSiena: 50 Scottish adolescents, friendship nominations at three waves, plus alcohol use and smoking measured at each wave. It needs no download and no political-science background, which matters given how many of you are not political scientists.

btergm wants three things:

The dependent variable as a list of network objects, one per wave;
Nodal covariates as vertex attributes on each network in that list;
Dyadic covariates as a list of matrices, one per wave.

That is the entire data-preparation problem. Most of the pain people report with TERGMs is this step, and most of that pain comes from doing it by hand for each wave instead of writing a loop.

suppressMessages({
  library(RSiena)
  library(network)
  library(btergm)
})

waves <- list(s501, s502, s503)                    # three 50x50 adjacency matrices
nets  <- lapply(waves, function(m) network(m, directed = TRUE))

# Nodal covariates: alcohol and smoking, measured at EACH wave, so they go on
# the matching network. This is a loop, not eight copy-pasted lines.
for (i in seq_along(nets)) {
  network::set.vertex.attribute(nets[[i]], "drink", s50a[, i])
  network::set.vertex.attribute(nets[[i]], "smoke", s50s[, i])
}

sapply(nets, network.edgecount)
#> [1] 113 116 122
sapply(nets, network.size)
#> [1] 50 50 50

Stuck? set.vertex.attribute and R’s copy semantics

If you have used R for a while, that loop should look wrong. In base R,

for (i in 1:3) f(x[[i]])

cannot modify x, because R copies on assignment. network objects break that rule: they carry an external reference, and set.vertex.attribute() modifies in place, returning nothing useful. So the loop does work, and you do not need to write nets[[i]] <- set.vertex.attribute(...).

This trips people up in both directions. It is also why you should always use the network:: prefix here. btergm and several other packages export functions with the same names, and which one you get depends on load order. Being explicit costs nine characters and saves an afternoon.

Check that the attributes actually landed. Do this every single time; a silently missing vertex attribute produces a model that fits fine and means nothing.

network::list.vertex.attributes(nets[[1]])
#> [1] "drink"        "na"           "smoke"        "vertex.names"
checkpoint(waves = length(nets),
           edges = sapply(nets, network.edgecount),
           mean_drink_w1 = mean(nets[[1]] %v% "drink"))
#> ------------------------------------------------------------------
#> CHECKPOINT: waves = 3   |   edges = 113, 116, 122   |   mean_drink_w1 = 2.88
#> ------------------------------------------------------------------

3. Looking at it before modelling it

The same rule as everywhere else in this course: look at the data before you fit anything. For panel networks, networkDynamic plus tsna gives you descriptives in about four lines.

suppressMessages({
  library(networkDynamic)
  library(tsna)
})

dyn <- networkDynamic(network.list = nets)
#> Neither start or onsets specified, assuming start=0
#> Onsets and termini not specified, assuming each network in network.list should have a discrete spell of length 1
#> Argument base.net not specified, using first element of network.list instead
#> Created net.obs.period to describe network
#>  Network observation period info:
#>   Number of observation spells: 1 
#>   Maximal time range observed: 0 until 3 
#>   Temporal mode: discrete 
#>   Time unit: step 
#>   Suggested time increment: 1
tSnaStats(dyn, "gden")     # density, wave by wave
#> Time Series:
#> Start = 0 
#> End = 3 
#> Frequency = 1 
#>        Series 1
#> [1,] 0.04612245
#> [2,] 0.04734694
#> [3,] 0.04979592
#> [4,]         NA
tSnaStats(dyn, "gtrans")   # transitivity, wave by wave
#> Time Series:
#> Start = 0 
#> End = 3 
#> Frequency = 1 
#>       Series 1
#> [1,] 0.3873874
#> [2,] 0.3728814
#> [3,] 0.4892857
#> [4,]        NA

tErgmStats() is the useful one: it computes any ERGM statistic at every wave, so you can see which structural features are actually moving before you decide to model them.

tErgmStats(dyn, "~ edges + mutual + triangle + idegree(2)")
#> Time Series:
#> Start = 0 
#> End = 3 
#> Frequency = 1 
#>   edges mutual triangle idegree2
#> 0   113     39      107       16
#> 1   116     35      116       15
#> 2   122     45      181       13
#> 3     0      0        0        0

Density and transitivity barely move; edges creep up. That is worth knowing before you fit, because it tells you the interesting variation is in which ties exist, not in how many. A model that only explains the number of edges will look fine and tell you nothing.

And the picture:

par(mfrow = c(1, 2), mar = c(1, 1, 2, 1))
p <- plot(network.extract(dyn, at = 0), main = "wave 1", vertex.cex = 1.3)
plot(network.extract(dyn, at = 2), main = "wave 3", coord = p, vertex.cex = 1.3)

Reusing p as coord puts every node in the same place in both panels, so the eye sees the edges change instead of the layout change. Always do this when comparing waves. It only works when the node set is identical across waves, which here it is. Section 7 handles a case where it is not.

4. Fitting a TERGM

btergm estimates the pooled coefficients by maximum pseudolikelihood, then resamples whole time points to describe temporal uncertainty. That second step needs a reasonably long panel. Increasing R repeats the resampling more often, but it does not manufacture new time points. There is also mtergm(), which uses MCMC-MLE and is what you would want if the model is small and you care about a likelihood. btergm() is useful for larger networks and longer panels, but its bootstrap intervals are not trustworthy in the three-wave teaching example below.

Start without the temporal terms, so you can see what they change.

set.seed(6886)
m_cross <- fit_or_load("tergm_s50_cross",
  btergm(nets ~ edges + mutual + gwesp(0.25, fixed = TRUE) +
                nodecov("drink") + absdiff("drink") + nodematch("smoke"),
         R = 200, verbose = FALSE))
round(coef(m_cross), 3)
#>                edges               mutual gwesp.OTP.fixed.0.25 
#>               -4.858                3.689                1.191 
#>        nodecov.drink        absdiff.drink      nodematch.smoke 
#>                0.061               -0.237                0.269

Note gwesp(0.25, fixed = TRUE). The fixed = TRUE is not optional and not decoration. As of ergm 4, a curved term with fixed = FALSE ignores the decay value you typed and fits a 16-parameter curved model instead. This is the most common silent error in inherited ERGM code, TERGM included, and btergm inherits it because it calls ergm’s term machinery.

Now add the two temporal terms.

set.seed(6886)
m_temp <- fit_or_load("tergm_s50_temp",
  btergm(nets ~ edges + mutual + gwesp(0.25, fixed = TRUE) +
                nodecov("drink") + absdiff("drink") + nodematch("smoke") +
                memory(type = "stability") + delrecip(),
         R = 200, verbose = FALSE))
round(coef(m_temp), 3)
#>                 edges                mutual  gwesp.OTP.fixed.0.25 
#>                -4.008                 2.955                 1.001 
#>         nodecov.drink         absdiff.drink       nodematch.smoke 
#>                 0.045                -0.161                 0.221 
#>   edgecov.memory[[i]] edgecov.delrecip[[i]] 
#>                 0.911                 0.528

4.1 Reading it, and the reason this document exists

Put the two side by side.

cmp <- merge(
  data.frame(term = names(coef(m_cross)), no_temporal = round(coef(m_cross), 3)),
  data.frame(term = names(coef(m_temp)),  with_temporal = round(coef(m_temp), 3)),
  by = "term", all = TRUE, sort = FALSE)
cmp
#>                    term no_temporal with_temporal
#> 1                 edges      -4.858        -4.008
#> 2                mutual       3.689         2.955
#> 3  gwesp.OTP.fixed.0.25       1.191         1.001
#> 4         nodecov.drink       0.061         0.045
#> 5         absdiff.drink      -0.237        -0.161
#> 6       nodematch.smoke       0.269         0.221
#> 7   edgecov.memory[[i]]          NA         0.911
#> 8 edgecov.delrecip[[i]]          NA         0.528

Look at what happens to the point estimates for mutual and gwesp. Both shrink substantially the moment you let the model know what the network looked like last wave. The point estimate for memory is also large, so dyadic stability is the biggest association in this fitted specification.

That shift is a compact reason to take temporal specification seriously. It is consistent with the cross-sectional terms partly absorbing the fact that last term’s friendships are still there, but the table alone cannot separate that story from coefficient rescaling across nonlinear models. The temporal terms change the question from a description of a photograph to a model of wave-to-wave association.

delrecip is positive in the fitted model: a tie extended to you last wave is associated with a higher conditional log-odds of extending one back this wave. That is a statement about sequence, and no cross-sectional model of any kind can make it. It is not yet an inferential claim, because two observed transitions are nowhere near enough to support the temporal bootstrap interval.

bootstrap_support <- function(fit, digits = 8L) {
  draws <- fit@boot$t
  c(requested_draws = nrow(draws),
    distinct_coefficient_vectors = nrow(unique(round(draws, digits))))
}

s50_support <- bootstrap_support(m_temp)
checkpoint(transitions = length(nets) - 1L,
           bootstrap_draws = s50_support[["requested_draws"]],
           distinct_coefficient_vectors = s50_support[["distinct_coefficient_vectors"]])
#> ------------------------------------------------------------------
#> CHECKPOINT: transitions = 2   |   bootstrap_draws = 200   |   distinct_coefficient_vectors = 3
#> ------------------------------------------------------------------

Why 200 bootstrap draws do not solve a two-transition problem

The fit reports 200 bootstrap rows, but after rounding away numerical noise there are only three distinct coefficient vectors. With two transition blocks, the temporal bootstrap has only three unordered resampling compositions available. Most of those 200 rows are repeats.

That is not a software bug. It is the data telling you that the uncertainty calculation has almost no temporal support. The pooled MPLE point estimates are still useful for walking through specification and interpretation, but a percentile interval, a standard error, or an “excludes zero” claim from this fit would look much more informative than it is. More waves, not a larger R, are what fix this problem. The same warning applies to every short-panel btergm fit later in this document.

A coefficient-comparison caveat you should already be twitchy about

I just compared coefficients across two nested specifications, which is exactly the thing the ERGM session told you not to do. Fair.

The defence is narrow and I want to be explicit about it: I am not claiming that the 0.7-point drop in mutual (3.69 → 2.95) is an effect size, and I am not testing anything. I am pointing at a large specification sensitivity. The size of the shift does not prove that rescaling is harmless, which is exactly the warning in Duxbury (2023). Read it as a qualitative prompt to investigate, not as a quantity with its own causal or inferential meaning.

If you want the defensible version, compare average marginal effects rather than coefficients, using ergMargins on an mtergm() fit.

Depth: summary() on a btergm object prints twice

summary(m_temp) prints its coefficient table two times over. This is a quirk of the S4 show method, not a sign anything went wrong. Avoid it by using coef() and confint() separately, as above, or by wrapping: invisible(summary(m_temp)).

5. Goodness of fit, including when it is bad news

Same epistemic move as ergm::gof(): simulate networks from the fitted model, compute statistics that are not in the model, and see whether the simulations look like the data. The call below is an in-sample check because it uses the fitted waves as its targets. btergm::gof() can also do genuine out-of-sample prediction, but only when you supply a held-out target and the corresponding forecasting formula explicitly. ROC and precision-recall curves from the default call are classification diagnostics for the fitted waves, not forecasts.

set.seed(6886)
g <- fit_or_load("tergm_s50_gof",
  btergm::gof(m_temp, nsim = 50,
              statistics = c(btergm::esp, btergm::geodesic, btergm::rocpr)))

g
#>   obs: mean median min max sim: mean median min max    Pr(>z)    
#> 0      42.0   42.0  34  50     39.43     39  25  53 0.8021918    
#> 1      45.0   45.0  44  46     44.45     44  15  64 0.7159375    
#> 2      28.5   28.5  22  35     15.38     13   0  43 0.2859865    
#> 3       3.5    3.5   0   7      2.05      1   0  12 0.7499871    
#> 4       0.0    0.0   0   0      0.18      0   0   3 0.0002937 ***
#> 5       0.0    0.0   0   0      0.00      0   0   0 1.0000000
#>     obs: mean median  min  max sim: mean median min  max  Pr(>z)  
#> 1       119.0  119.0  116  122    101.49  103.0  61  132 0.04982 *
#> 2       119.5  119.5  105  134    150.36  151.0  53  265 0.25273  
#> 3        96.5   96.5   67  126    184.64  183.0  41  402 0.18599  
#> 4        71.5   71.5   25  118    174.90  178.0  22  361 0.26074  
#> 5        58.5   58.5   13  104    133.92  133.0  15  296 0.34098  
#> 6        49.0   49.0    9   89     89.96   90.5   4  210 0.49089  
#> 7        33.5   33.5    4   63     53.67   53.5   0  156 0.61772  
#> 8        24.0   24.0    1   47     28.67   23.0   0   96 0.87260  
#> 9        17.0   17.0    0   34     13.67    7.5   0   98 0.87697  
#> 10       14.5   14.5    0   29      6.20    2.0   0   76 0.66893  
#> 11        5.0    5.0    0   10      2.48    0.0   0   53 0.70272  
#> 12        3.5    3.5    0    7      0.98    0.0   0   31 0.60223  
#> 13        2.5    2.5    0    5      0.37    0.0   0   14 0.55031  
#> 14        0.0    0.0    0    0      0.14    0.0   0    9 0.14967  
#> 15        0.0    0.0    0    0      0.04    0.0   0    4 0.31975  
#> 16        0.0    0.0    0    0      0.01    0.0   0    1 0.31975  
#> 17        0.0    0.0    0    0      0.00    0.0   0    0 1.00000  
#> Inf    1836.0 1836.0 1568 2104   1508.50 1473.5 768 2223 0.43399
#>   ROC model ROC random  PR model PR random
#> 1 0.8052523   0.518364 0.3827564 0.0488066

Read the tables from the right-hand column: small simulation-based p-values flag statistics for which the fitted model rarely reproduces the observed value. Most rows here are comfortable. The exceptions sit in sparse tail bins, where the observed count is zero and the simulations occasionally produce one or two. Do not automatically call those rows Monte Carlo noise, and do not turn each bin into a separate hypothesis test. Look for a persistent pattern across related bins, rerun with more simulations to check Monte Carlo stability, and remember that this diagnostic conditions on a point estimate whose short-panel uncertainty is not being propagated.

The two tie-classification summaries are:

g$`Tie prediction`$auc.roc
#> [1] 0.8052523
g$`Tie prediction`$auc.pr
#> [1] 0.3827564

An ROC around 0.81 against a chance baseline of 0.5, and a precision-recall AUC around 0.37 against a chance baseline equal to network density (about 0.05). For a sparse network, the PR number is usually more informative because ROC can look strong when a model mainly gets the many zeros right. Report both, label them as in-sample, and put the PR baseline next to the PR AUC.

Fit the model, then check it. Every time.

The original version of this material fitted a TERGM to diplomatic ties and left the gof() call commented out. When I uncommented it (§7) the model turned out to reproduce tie placement very well and to get the network’s path structure noticeably wrong. That mixed verdict that changes how you would write up the result. Nobody had looked, so nobody knew.

That is not a story about that particular model. It is a story about what happens when the check is one keystroke away and nobody presses it. A document that fits and never checks teaches a habit that will eventually cost you a paper.

6. TERGM vs SAOM on identical data

There is a live methodological argument about whether you should be doing this at all, or whether you should be fitting a stochastic actor-oriented model instead. You met SAOMs on the Thursday. The dispute is worth understanding, and the cheapest way to understand it is to fit both to the same data and look.

We will use the Knecht classroom data: 26 Dutch pupils, four waves of friendship, plus sex and a matrix of who went to primary school together. It ships with btergm.

6.1 A composition problem, and why it is not a nuisance

Knecht has pupils entering and leaving. Absent actors are coded 10 (structurally missing) and there are scattered incidental NAs. A TERGM needs the risk set at each wave to be encoded deliberately. btergm has node-adjustment and structural-zero machinery for changing rosters, but for this clean same-data comparison we use a simpler estimand and restrict to pupils present at all four waves.

suppressMessages(library(btergm))
data(knecht)

absent <- function(m) apply(m, 1, function(r) all(is.na(r) | r == 10))
gone   <- Reduce(`|`, lapply(friendship, absent))
keep   <- which(!gone)
c(total = length(gone), kept = length(keep), dropped = sum(gone))
#>   total    kept dropped 
#>      26      22       4

knets <- lapply(friendship, function(m) {
  m <- m[keep, keep]
  m[m == 10] <- NA
  diag(m) <- 0
  network(m, directed = TRUE)
})
sx <- demographics$sex[keep]
for (i in seq_along(knets)) network::set.vertex.attribute(knets[[i]], "sex", sx)
sapply(knets, network.edgecount)
#>  t1  t2  t3  t4 
#>  71  93 109 102

Dropping four of twenty-six pupils is a real cost and you should feel it. SAOMs can model composition change directly via sienaCompositionChange(), because the model is written at the level of actors taking steps rather than at the level of a graph. That is a genuine advantage of the actor-oriented formulation and not a small one. If you want the mechanics, see the companion document Composition change in SAOMs.

6.2 Both models, same data

set.seed(6886)
k_tergm <- fit_or_load("tergm_knecht",
  btergm(knets ~ edges + mutual + ttriple + nodematch("sex") +
                 memory(type = "stability") + delrecip(),
         R = 200, verbose = FALSE))
round(coef(k_tergm), 3)
#>                 edges                mutual               ttriple 
#>                -3.227                 0.899                 0.321 
#>         nodematch.sex   edgecov.memory[[i]] edgecov.delrecip[[i]] 
#>                 0.828                 0.862                 0.252
knecht_support <- bootstrap_support(k_tergm)
checkpoint(transitions = length(knets) - 1L,
           bootstrap_draws = knecht_support[["requested_draws"]],
           distinct_coefficient_vectors = knecht_support[["distinct_coefficient_vectors"]])
#> ------------------------------------------------------------------
#> CHECKPOINT: transitions = 3   |   bootstrap_draws = 200   |   distinct_coefficient_vectors = 10
#> ------------------------------------------------------------------

There are three observed transitions here. The 200 bootstrap rows therefore collapse to only ten distinct coefficient vectors after rounding, the ten unordered ways to resample three transition blocks. Use this fit to compare model structure and point-estimate signs. Do not use its bootstrap percentiles as publication-ready uncertainty.

suppressMessages(library(RSiena))

arr <- array(NA_integer_, c(length(keep), length(keep), 4))
for (i in 1:4) {
  m <- friendship[[i]][keep, keep]
  m[m == 10] <- NA
  diag(m) <- 0
  arr[, , i] <- m
}

fr    <- sienaDependent(arr)
sexC  <- coCovar(demographics$sex[keep])
primC <- coDyadCovar(primary[keep, keep])
dat   <- sienaDataCreate(fr, sexC, primC)

eff <- getEffects(dat)
eff <- includeEffects(eff, transTrip, cycle3)
#>   effectNumber effectName          shortName include fix   test  initialValue
#> 1 23           transitive triplets transTrip TRUE    FALSE FALSE          0  
#> 2 43           3-cycles            cycle3    TRUE    FALSE FALSE          0  
#>   parm
#> 1 0   
#> 2 0
eff <- includeEffects(eff, sameX, interaction1 = "sexC")
#>   effectNumber effectName shortName include fix   test  initialValue parm
#> 1 340          same sexC  sameX     TRUE    FALSE FALSE          0   0
alg <- sienaAlgorithmCreate(projname = NULL, seed = 6886, n3 = 1000)
#> If you use this algorithm object, siena07 will create/use an output file Siena.txt .

k_saom <- fit_or_load("saom_knecht",
  siena07(alg, data = dat, effects = eff,
          batch = TRUE, verbose = FALSE, silent = TRUE))

data.frame(effect = k_saom$effects$effectName,
           est    = round(k_saom$theta, 3),
           se     = round(sqrt(diag(k_saom$covtheta)), 3))
#>                effect    est    se
#> 1 outdegree (density) -2.079 0.150
#> 2         reciprocity  1.407 0.193
#> 3 transitive triplets  0.364 0.041
#> 4            3-cycles -0.398 0.073
#> 5           same sexC  0.675 0.132

Convergence first, always. A SAOM you have not checked is not a result:

round(as.numeric(k_saom$tconv.max), 4)
#> [1] 0.0572

Under 0.25, so this run clears the usual RSiena convergence screen. That check is necessary, but it does not by itself make a model publishable. Specification, uncertainty, and goodness of fit still need their own checks.

6.3 What the comparison actually shows

data.frame(
  quantity     = c("reciprocity", "transitivity", "sex homophily"),
  TERGM        = round(coef(k_tergm)[c("mutual", "ttriple", "nodematch.sex")], 3),
  SAOM         = round(k_saom$theta[match(c("reciprocity", "transitive triplets", "same sexC"),
                                          k_saom$effects$effectName)], 3),
  row.names    = NULL)
#>        quantity TERGM  SAOM
#> 1   reciprocity 0.899 1.407
#> 2  transitivity 0.321 0.364
#> 3 sex homophily 0.828 0.675

The fitted point-estimate patterns agree in sign. In both models, reciprocity, transitivity, and same-sex ties receive positive coefficients. That is a useful descriptive concordance across two model families with genuinely different assumptions. It is not an inferential agreement claim, especially given the TERGM’s three-transition bootstrap, and the coefficient magnitudes are not comparable across the families.

The interpretations do not agree, and the coefficients are not the same quantity:

The TERGM’s mutual is a contribution to the conditional log-odds of a tie, given the rest of the graph. It is defined at the level of the graph.
The SAOM’s reciprocity is a contribution to an actor’s objective function at a micro-step: how much better off actor \(i\) is, in the model’s terms, for having a reciprocated tie when it is \(i\)’s turn to act. It is defined at the level of the transition.

You cannot put those in the same table and call one bigger. Block, Stadtfeld & Snijders (2019) is the paper that makes this precise, and it is worth reading before you choose between the two families.

Depth: the Leifeld–Cranmer / Block exchange, and where it actually stands

This got sharp in print, and the sequence matters.

Leifeld & Cranmer (2019) compared the two families theoretically and empirically, using held-out-wave predictive accuracy as the criterion. Their own conclusion was deliberately non-triumphal: “Though we do not aim to make a general claim about either being superior to the other across all specifications… each model out-predicts the other one the more the specific assumptions of the respective model are met.”

Block, Hollway, Stadtfeld, Koskinen & Snijders (2022) showed the TERGM specification in that comparison used nodal attributes computed from the outcome networks’ observed degrees rather than endogenous ERGM equivalents, so 3 of 13 model terms were contaminated by the held-out wave. The out-of-sample prediction was, in part, predicting the future using information from the future. They argue further that tie-level predictive accuracy is the wrong criterion in the first place, since network science cares about endogenous processes, not about which specific dyads are filled.

Leifeld & Cranmer (2022) acknowledged the coding error, and pressed a deeper point: the SAOM hard-codes a substantive theory involving sequential one-tie-at-a-time micro-steps, one-period look-ahead, and every actor knowing the whole network. They argue that embedded theory should be falsifiable. They contend held-out prediction is the only feasible way to test it, and that Block et al. fixed the specification rather than the design.

Where it stands: the factual dispute is settled. The error was real and was conceded. The methodological dispute is not. Whether predictive fit is a legitimate criterion for comparing generative network models is genuinely open, and both positions are defensible. Anyone who tells you one side won is selling something.

7. An applied example with a larger, messier network

Diplomatic recognition (Duque 2018), optional and IR-flavoured

Everything above runs on classroom friendship because that needs no translation. This section is a political-science application, folded because it is not necessary and because the substance costs non-IR readers something to decode. If you open it, here is the one-sentence gloss you need: states establish embassies in other states, and where they choose to put them is treated as a signal of who they regard as a legitimate, high-status member of international society.

The data are eight waves (1970–2005, every five years) of who hosts an embassy from whom, which makes a directed network that grows from 134 states in 1970 to 188 in 2005. Duque (2018) argues status recognition is social: you get recognised because of your position in the network of recognition, not only because of your material attributes. That is a claim about endogenous structure, so a TERGM is one reasonable model family to consider.

suppressMessages(library(networkdata))
data(duqueData)
c(waves = length(dipl_ties), dim_wave1 = dim(dipl_ties[[1]])[1],
  dim_wave8 = dim(dipl_ties[[8]])[1])
#>     waves dim_wave1 dim_wave8 
#>         8       134       188

The node set changes across waves because states enter and leave the international system. Some data frames also retain a state’s name while coding its entire outgoing row and incoming column as missing for that wave. Initialising a pooled actor-by-actor matrix with zeros would quietly turn both kinds of absence into observed non-ties.

For this teaching fit, we use a transparent complete-roster estimand: states with observed incoming and outgoing diplomacy data in every one of the four selected waves. This avoids inventing non-ties, but it changes the target population to persistent states. A full applied analysis should model composition directly or use carefully validated structural-zero machinery, then assess how much the roster restriction changes the result.

W <- 4L                                    # first four waves: 1970-1985

present_ids <- function(df) {
  m <- as.matrix(df)
  missing <- is.na(m) | m == "."
  diag(missing) <- TRUE                    # the diagonal cannot establish presence
  ids <- names(df)
  ids[!apply(missing, 1, all) & !apply(missing, 2, all)]
}

wave_actors <- lapply(dipl_ties[1:W], present_ids)
pooled_actors <- sort(unique(unlist(lapply(dipl_ties[1:W], names))))
actors <- sort(Reduce(intersect, wave_actors))

to_mat <- function(df, actors) {
  m <- as.matrix(df[actors, actors, drop = FALSE])
  m[m == "."] <- NA                        # "." is this dataset's missing code
  storage.mode(m) <- "integer"
  diag(m) <- 0L
  if (anyNA(m)) {
    stop("The complete roster still contains missing diplomacy dyads.", call. = FALSE)
  }
  m
}

dipl_m <- lapply(dipl_ties[1:W], to_mat, actors = actors)
c(pooled_names = length(pooled_actors),
  complete_roster = length(actors),
  excluded = length(setdiff(pooled_actors, actors)),
  missing_cells_after_restriction = sum(vapply(
    dipl_m, function(m) sum(is.na(m)), integer(1))))
#>                    pooled_names                 complete_roster 
#>                             163                             132 
#>                        excluded missing_cells_after_restriction 
#>                              31                               0

Four data mistakes this avoids

The earlier construction had four ways to produce output that looked fine while changing the dependent variable:

Treating absence as a non-tie. Taking the union of state names, filling a large matrix with zeros, and writing each smaller wave into it codes every not-yet-entered, departed, or fully missing state-dyad as an observed zero. Those dyads were not at risk of a diplomatic tie in the observed wave.
Positional indexing into a name-dimensioned array. dipl[i, j, t] <- val where i and j came from looping over the wave’s rows, but dipl was dimensioned by the pooled actor list. Since wave 1 has 134 states and the pooled four-wave list has 163, everything from about actor 3 onward was written to the wrong cell. Subsetting with state names prevents this.
Silent symmetrisation. Writing both dipl[i,j,t] <- val and dipl[j,i,t] <- val forces a symmetric matrix. Diplomatic recognition is directed. That A hosts an embassy from B does not mean B hosts one from A, and the whole “who recognises whom” question depends on the asymmetry. That line destroys roughly two thousand cells per wave, which is to say it destroys the dependent variable.
Coercing a missing code to NA without noticing. as.numeric(as.character(x)) on this data frame silently turns "." into NA with a warning nobody reads. There are 588 such cells in the original 1975 frame. Here they help identify an absent state before we restrict the roster. They are never recoded as observed zeros.

None of these mistakes throws an error. All four change your answer. This is the argument for checking dim(), table(), and isSymmetric() after every reshaping step, every time.

sapply(dipl_m, function(m) sum(m != t(m), na.rm = TRUE))     # asymmetric cells: should be >> 0
#> 1970 1975 1980 1985 
#> 1988 1920 1896 1956
table(unlist(dipl_m), useNA = "ifany")
#> 
#>     0     1 
#> 51840 17856

Thousands of asymmetric cells per wave, exactly as a directed network should have. If that row had come back zeros you would know you had symmetrised something.

Now build networks and fit. The dyadic covariates, contiguity and alliance, need the same name-based restriction.

dnets <- lapply(dipl_m, network, directed = TRUE, matrix.type = "adjacency")

align <- function(df, actors) {
  m <- as.matrix(df[actors, actors, drop = FALSE])
  storage.mode(m) <- "numeric"
  if (anyNA(m)) stop("A dyadic covariate is missing on the complete roster.", call. = FALSE)
  m
}
contig_m <- lapply(contig[1:W], align, actors = actors)
allies_m <- lapply(allies[1:W], align, actors = actors)

yrs <- c(1970, 1975, 1980, 1985)
for (i in seq_len(W)) {
  dem <- polity$dem_dum[match(actors, polity$ccode[polity$year == yrs[i]])]
  if (anyNA(dem)) stop("Democracy is missing on the complete roster.", call. = FALSE)
  network::set.vertex.attribute(dnets[[i]], "dem", dem)
}
checkpoint(nodes = sapply(dnets, network.size),
           edges = sapply(dnets, network.edgecount),
           missing_edges = sapply(dnets, network.naedgecount))
#> ------------------------------------------------------------------
#> CHECKPOINT: nodes = 132, 132, 132, 132   |   edges = 3604, 4406, 4868, 4978   |   missing_edges = 0, 0, 0, 0
#> ------------------------------------------------------------------

suppressMessages(library(btergm))
set.seed(6886)
d_fit <- fit_or_load("tergm_duque_complete_roster",
  btergm(dnets ~ edges + mutual + gwesp(0.25, fixed = TRUE) +
                 nodeicov("dem") + absdiff("dem") +
                 edgecov(allies_m) + edgecov(contig_m) +
                 memory(type = "stability") + delrecip(),
         R = 200, verbose = FALSE))
round(coef(d_fit), 3)
#>                 edges                mutual  gwesp.OTP.fixed.0.25 
#>                -2.289                 2.442                 0.541 
#>          nodeicov.dem           absdiff.dem edgecov.allies_m[[i]] 
#>                -0.086                 0.127                 0.621 
#> edgecov.contig_m[[i]]   edgecov.memory[[i]] edgecov.delrecip[[i]] 
#>                 0.455                 1.654                 0.385
duque_support <- bootstrap_support(d_fit)
checkpoint(transitions = length(dnets) - 1L,
           bootstrap_draws = duque_support[["requested_draws"]],
           distinct_coefficient_vectors = duque_support[["distinct_coefficient_vectors"]])
#> ------------------------------------------------------------------
#> CHECKPOINT: transitions = 3   |   bootstrap_draws = 200   |   distinct_coefficient_vectors = 10
#> ------------------------------------------------------------------

This is the slow fit in the document.

In this complete-roster fit, memory has the largest positive point estimate, which is consistent with embassies being expensive and durable. The fitted coefficients for reciprocity, contiguity, and alliance are also positive. The democracy terms are smaller and depend on their coding and the rest of the specification. These are model-conditional point patterns, not confirmed effects and not a replication of Duque’s full analysis.

The bootstrap diagnostic shows the same short-panel problem again. Three transitions provide only ten distinct resampling compositions, however large we make R. The four-wave subset is useful for teaching data preparation and interpretation, but its bootstrap intervals cannot carry a substantive inference.

And now the check that the original never ran:

set.seed(6886)
d_gof <- fit_or_load("tergm_duque_complete_roster_gof",
  btergm::gof(d_fit, nsim = 30, statistics = c(btergm::geodesic, btergm::rocpr)))

with(d_gof$`Tie prediction`,
     c(auc_roc = auc.roc,
       auc_pr = auc.pr,
       random_graph_pr = auc.pr.rgraph))
#>         auc_roc          auc_pr random_graph_pr 
#>       0.9290458       0.8793259       0.2784350

round(d_gof$`Geodesic distances`$stats, 3)
#>     obs: mean median   min   max sim: mean median   min   max Pr(>z)
#> 1    4750.667   4868  4406  4978  4748.389   4834  4097  5261  0.991
#> 2   11860.667  11784 11778 12020 12289.722  12310 11834 12759  0.018
#> 3     375.667    398   215   514   230.400    236   106   413  0.236
#> 4       0.667      1     0     1     0.200      0     0     5  0.294
#> 5       0.000      0     0     0     0.000      0     0     0  1.000
#> Inf   304.333    131   131   651    23.289      0     0   262  0.246

Read this carefully, because it is a good example of a GOF table that needs interpreting rather than scoring.

The ROC is high in this run, but it is an in-sample classification summary. Stability and the many correctly classified non-ties can make ROC look impressive. It does not validate the temporal specification, and it is not evidence of out-of-sample forecasting performance.

The geodesic table asks a different and more structural question. Compare the observed count in each path-length row with the simulated distribution. A run of adjacent rows outside the simulation envelope is more concerning than an isolated sparse-tail row. The Inf row is especially useful because it records pairs of states with no directed path between them, so it shows whether the model makes the diplomatic system too connected or too fragmented.

There are only thirty simulations at each transition here, ninety networks in total. Raising nsim can tell you whether a GOF pattern is stable to simulation noise, but it cannot add transitions or propagate the severe parameter uncertainty from a three-transition bootstrap. The honest workflow is to show these model-conditional diagnostics beside the point estimates, describe persistent mismatches, and then refit on a longer panel before making substantive inferential claims.

8. Exercises

Refit the s50 model from §4, replacing memory(type = "stability") with memory(type = "autoregression").

set.seed(6886)
m_auto <- btergm(nets ~ edges + mutual + gwesp(0.25, fixed = TRUE) +
                        nodecov("drink") + absdiff("drink") + nodematch("smoke") +
                        memory(type = "__________") + delrecip(),
                 R = 200, verbose = FALSE)

Then answer in one sentence: "stability" codes a dyad 1 if it is in the same state as last wave (tie→tie or non-tie→non-tie), while "autoregression" codes a dyad 1 only if there was a tie last wave. Which one do you get a bigger coefficient for? Work out the answer from the coding before you look at the number.

Add timecov() to the s50 model and interpret the result. Then explain, without fitting anything, why you should not trust it much here.

Hint for the second part: timecov() estimates a trend, and a TERGM on \(T\) waves has \(T-1\) transitions to learn it from. Count the transitions you have.

Second stretch, harder: the §4 comparison showed mutual shrinking when temporal terms were added. Design a check that distinguishes “the temporal terms absorbed real persistence” from “the coefficients rescaled for the reason Duxbury describes.” What would you have to compute?

set.seed(6886)
m_auto <- fit_or_load("tergm_s50_auto",
  btergm(nets ~ edges + mutual + gwesp(0.25, fixed = TRUE) +
                nodecov("drink") + absdiff("drink") + nodematch("smoke") +
                memory(type = "autoregression") + delrecip(),
         R = 200, verbose = FALSE))

data.frame(
  term        = c("memory", "mutual", "gwesp"),
  stability   = round(coef(m_temp)[c("edgecov.memory[[i]]", "mutual", "gwesp.OTP.fixed.0.25")], 3),
  autoregress = round(coef(m_auto)[c("edgecov.memory[[i]]", "mutual", "gwesp.OTP.fixed.0.25")], 3),
  row.names   = NULL)
#>     term stability autoregress
#> 1 memory     0.911       1.821
#> 2 mutual     2.955       2.955
#> 3  gwesp     1.001       1.001

Core. "autoregression" gives the larger coefficient, roughly double. If you predicted "stability" because it “covers more dyads,” you reasoned about coverage when you should have reasoned about contrast.

Here is why. This network’s density is about 5%, so 95% of dyads are non-ties. The "stability" covariate is 1 for almost every dyad in the data, because almost every dyad was a non-tie last wave and is a non-tie now. A covariate that is nearly constant carries very little information. Worse, edges is already soaking up the “most dyads are empty” baseline, so "stability" is left estimating a modest increment over a baseline that has already been accounted for.

"autoregression" is 1 only for the 5% of dyads that actually had a tie. That is a sharp contrast against the 95% that did not, and its coefficient measures a fitted difference in conditional log-odds between dyads with and without a prior tie.

Both are legitimate; they are answering different questions. The practical consequences are that their magnitudes are not comparable. Never report one as if it were the other. On sparse networks, "autoregression" is usually the more interpretable of the two, while "stability" is the more complete control precisely because it also constrains the non-ties.

Note that the mutual and gwesp point estimates come out identical to the third decimal under both. That is worth pausing on: swapping the memory flavour changes what memory itself measures, but does not disturb the rest of the model. Compare that with §4.1, where adding memory at all moved mutual by about three-quarters of a point. In this example, adding a temporal term creates much more point-estimate movement than swapping these two memory codings. That comparison is descriptive, not an uncertainty statement.

Stretch, part 1. With three waves you have two transitions. timecov() is fitting a time trend to two points. It will estimate, it will report a standard error, and the number is close to meaningless. Two points determine a line exactly, with nothing left over to say whether a line was the right shape. You want five or six waves before a time trend is a claim rather than an interpolation.

Stretch, part 2. You would need a quantity that is not on the log-odds scale, so that rescaling cannot move it. Average marginal effects are the standard answer: fit both specifications with mtergm() (which returns an ergm-compatible object) and compute ergMargins::ergm.AME() for mutual under each. If the AME in percentage points drops as much as the coefficient did, the shift is substantive. If the coefficient moves and the AME does not, you were watching residual rescaling.

Report the baseline density alongside the AME either way. On a network this sparse, absolute marginal effects are small numbers and mean nothing without the base rate next to them.

9. Where you would actually reach for the neighbours

The TERGM is one point in a space, and the terminal slides of the original deck listed the neighbours as bare URLs. Here is what each one buys you, so you can tell whether you want it.

tergm (the statnet package), and STERGM. Where btergm fits one equation to the transition, tergm() can fit separate formation and dissolution equations. The things that make ties appear need not be the things that make ties persist, and often are not. Reach for it when that separation is your actual question, e.g. “homophily drives who becomes friends, but proximity drives who stays friends.” One caution: stergm() is deprecated in tergm 4.2.2 in favour of tergm(), of which it is a special case. Older tutorials will tell you to call stergm(); do not.

Bergm: Bayesian ERGMs. Caimo & Friel (2011). Gives you full posteriors instead of point estimates plus asymptotic errors, which matters most exactly where ERGMs are hardest: small networks, near-degenerate specifications, and any situation where the likelihood surface is not remotely quadratic. Reach for it when your standard errors look too good to be true, or when you want to be honest about a curved term’s decay parameter instead of fixing it.

ergm.ego: inference from egocentric samples. Most survey data is not a complete network; it is a sample of respondents each reporting a handful of alters. ergm.ego estimates a population ERGM from that, with proper design weights. Reach for it whenever your data came from a survey rather than a bounded roster. For many of you this is the most relevant item on this list, and it is the one people most often do not know exists.

Hierarchical / block-structured ERGMs (bigergm). Permits dependence within blocks and not across them. This is “local dependence” in the Schweinberger–Handcock sense, which is what makes large networks tractable. Reach for it above roughly a thousand nodes, where plain MCMC-MLE stops finishing.

Relational event models (relevent, rem, goldfish, remstimate). If your data are timestamped events, such as emails, calls, citations, or meetings, rather than panel waves, you are not in TERGM/SAOM territory at all. Those models assume the between-wave process is unobserved and must be inferred. When you see the event sequence directly, you do not need to infer it, and most of the inferential burden disappears. Bianchi, Filippi-Mazzola, Lomi & Wit (2024) is the current overview; Butts, Lomi, Snijders & Stadtfeld (2023) frame why the event level matters. Reach for these when your natural unit is “a thing happened at 14:32” rather than “a tie existed in 2019.”

10. Readings

The 2026 syllabus assigns no dedicated TERGM reading. Day 11’s temporal block (C4) points here instead. If you are working through the deep dive on your own, these four are the set:

Leifeld, Cranmer & Desmarais (2018), “Temporal Exponential Random Graph Models with btergm,” Journal of Statistical Software. The software paper for everything in §4–5. Start here.
Hanneke, Fu & Xing (2010), “Discrete temporal models of social networks,” Electronic Journal of Statistics. The model in §1.1, done properly.
Block, Stadtfeld & Snijders (2019), “Forms of dependence: comparing SAOMs and ERGMs from basic principles,” Sociological Methods & Research. The graph-level vs transition-level distinction in §6.3. This is also the reference “Taking Dyads Seriously” leans on for its ERGM caveat, so it closes a loop back to the Thursday of week 2.
The Leifeld–Cranmer / Block exchange (2019, 2022). Read all three pieces in order or none of them; reading one is worse than reading none.

Cranmer & Desmarais (2011), which is on the syllabus for the ERGM day, names cross-sectionality as one of two chief limitations of the ERGM. Day 11’s C4 opened the fix for that one, and this document finishes it. Krivitsky (2012), also on the syllabus for the ERGM day, is the fix for the other.

Versions: expand if your numbers differ from mine

pk <- c("ergm", "network", "btergm", "networkDynamic", "tsna", "RSiena", "networkdata")
data.frame(package = pk,
           version = sapply(pk, function(p) as.character(packageVersion(p))),
           row.names = NULL)
#>          package version
#> 1           ergm  4.12.0
#> 2        network  1.20.0
#> 3         btergm  1.11.1
#> 4 networkDynamic  0.12.0
#> 5           tsna   0.3.6
#> 6         RSiena   1.5.0
#> 7    networkdata     0.1
cat("R", as.character(getRversion()), "\n")
#> R 4.3.3