---
title: "Temporal ERGMs"
subtitle: "Ancillary self-study · ICPSR Network Analysis: Advanced Topics"
author: "Shahryar Minhas"
date: last-modified
engine: knitr
editor: source
format:
  html:
    theme: [cosmo, _session.scss]
    toc: true
    toc-depth: 3
    toc-location: left
    toc-title: "Contents"
    embed-resources: true
    code-tools: true
    fig-width: 7
    fig-height: 4
    fig-align: center
knitr:
  opts_chunk:
    message: false
    warning: false
    comment: "#>"
    collapse: true
    dev: "png"
    dpi: 140
execute:
  cache: true
  freeze: false
  error: false
---

```{r}
#| label: setup
#| include: false
source("_helpers.R")
require_pkgs(c("ergm", "network", "btergm", "networkDynamic", "tsna", "RSiena", "networkdata"))
set.seed(6886)
```

{{< include _ancillary_header.qmd >}}

::: {.callout-note title="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.

::: {.callout-note collapse="true" title="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:

1. The dependent variable as a **`list` of `network` objects**, one per wave;
2. Nodal covariates as **vertex attributes on each network in that list**;
3. 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.

```{r}
#| label: data-setup
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)
sapply(nets, network.size)
```

::: {.callout-tip collapse="true" title="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,

```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.

```{r}
#| label: data-check
network::list.vertex.attributes(nets[[1]])
checkpoint(waves = length(nets),
           edges = sapply(nets, network.edgecount),
           mean_drink_w1 = mean(nets[[1]] %v% "drink"))
```

# 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.

```{r}
#| label: netdyn
suppressMessages({
  library(networkDynamic)
  library(tsna)
})

dyn <- networkDynamic(network.list = nets)
tSnaStats(dyn, "gden")     # density, wave by wave
tSnaStats(dyn, "gtrans")   # transitivity, wave by wave
```

`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.

```{r}
#| label: tergmstats
tErgmStats(dyn, "~ edges + mutual + triangle + idegree(2)")
```

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:

```{r}
#| label: snapshots
#| fig-width: 9
#| fig-height: 4.5
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.

```{r}
#| label: fit-notemporal
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)
```

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.

```{r}
#| label: fit-temporal
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)
```

## 4.1 Reading it, and the reason this document exists

Put the two side by side.

```{r}
#| label: compare
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
```

```{r}
#| label: mutual-shift
#| echo: false
mutual_cross <- coef(m_cross)[["mutual"]]
mutual_temp  <- coef(m_temp)[["mutual"]]
mutual_drop  <- mutual_cross - mutual_temp
```

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.

```{r}
#| label: bootstrap-support
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"]])
```

::: {.callout-important title="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.
:::

::: {.callout-warning title="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 `r round(mutual_drop, 1)`-point drop in `mutual` (`r round(mutual_cross, 2)` → `r round(mutual_temp, 2)`) 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.
:::

::: {.callout-note collapse="true" title="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.

```{r}
#| label: gof
#| output: false
set.seed(6886)
g <- fit_or_load("tergm_s50_gof",
  btergm::gof(m_temp, nsim = 50,
              statistics = c(btergm::esp, btergm::geodesic, btergm::rocpr)))
```

```{r}
#| label: gof-print
g
```

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:

```{r}
#| label: gof-roc
g$`Tie prediction`$auc.roc
g$`Tie prediction`$auc.pr
```

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.

::: {.callout-important title="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 `NA`s. 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.

```{r}
#| label: knecht-prep
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))
```

```{r}
#| label: knecht-nets
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)
```

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](saom_composition.html).

## 6.2 Both models, same data

```{r}
#| label: knecht-btergm
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)
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"]])
```

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.

```{r}
#| label: knecht-saom
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)
eff <- includeEffects(eff, sameX, interaction1 = "sexC")
alg <- sienaAlgorithmCreate(projname = NULL, seed = 6886, n3 = 1000)

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))
```

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

```{r}
#| label: knecht-conv
round(as.numeric(k_saom$tconv.max), 4)
```

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

```{r}
#| label: sidebyside
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)
```

**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.

::: {.callout-note collapse="true" title="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

::: {.callout-note collapse="true" title="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.

```{r}
#| label: duque-load
suppressMessages(library(networkdata))
data(duqueData)
c(waves = length(dipl_ties), dim_wave1 = dim(dipl_ties[[1]])[1],
  dim_wave8 = dim(dipl_ties[[8]])[1])
```

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.

```{r}
#| label: duque-prep
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))))
```

::: {.callout-warning title="Four data mistakes this avoids"}
The earlier construction had four ways to produce output that looked fine while changing the dependent variable:

1. **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.

2. **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.

3. **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.

4. **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.
:::

```{r}
#| label: duque-check
#| cache: false
sapply(dipl_m, function(m) sum(m != t(m), na.rm = TRUE))     # asymmetric cells: should be >> 0
table(unlist(dipl_m), useNA = "ifany")
```

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.

```{r}
#| label: duque-nets
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))
```

```{r}
#| label: duque-fit
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)
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"]])
```

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:

```{r}
#| label: duque-gof
#| output: false
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)))
```

```{r}
#| label: duque-gof-roc
#| cache: false
with(d_gof$`Tie prediction`,
     c(auc_roc = auc.roc,
       auc_pr = auc.pr,
       random_graph_pr = auc.pr.rgraph))
```

```{r}
#| label: duque-gof-geo
#| cache: false
round(d_gof$`Geodesic distances`$stats, 3)
```

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

::: {.panel-tabset}

## Core

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

```{r}
#| eval: false
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.

## Stretch

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?

## Solution

```{r}
#| label: ex-solution
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)
```

**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.

```{r}
#| label: sessioninfo
#| cache: false
#| code-fold: true
#| code-summary: "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)
cat("R", as.character(getRversion()), "\n")
```
