---
title: "ERGMs in applied political science"
engine: knitr
editor: source
subtitle: "Ancillary self-study · ICPSR Network Analysis: Advanced Topics"
---

```{r}
#| label: setup
#| include: false
source("_helpers.R")
require_pkgs(c("ergm", "network", "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). The taught session runs entirely on friendship data: adolescents and monks, so that nobody pays a translation tax on the substance while learning the method. This document is the other half of that bargain: two worked **political-science** applications, for the students who came for exactly this and who will get more out of the method when it is attached to a question they care about.

**Two applications:**

1. **Senate cosponsorship**: a small, clean, dyad-independent model. Good for seeing the whole pipeline (build the network, attach covariates, fit, interpret) with no MCMC in the way.
2. **Rebel-group cooperation in civil wars** (Gade et al. 2019), a pedagogical ERGM reanalysis of data from a published AME application. It adds endogenous structure and a real specification choice about shared partners, but it is not a replication of the published model.

**Why it is optional in 2026:** the taught session had to pick datasets that work for a room that is half non-political-science, and IR examples do not clear that bar. Here they are, done properly, for the people who want them. If you are *not* a political scientist, the method is identical to what you already saw: you can read this purely as "the same pipeline on messier data" and skip the substance.

**What you need:** `ergm`, `network`, `networkdata` (which carries both datasets: no `install_github` is executed here).
:::

# 1. Senate cosponsorship: the whole pipeline, no MCMC

## 1.1 The question and the data

Take the nine most ideologically extreme senators of the 109th Congress. Draw a directed edge from $i$ to $j$ when $i$ cosponsored at least two of $j$'s bills. The substantive question: do senators cosponsor each other because they are **ideologically similar**, or because they are **geographically close**, or both?

This maps onto three ERGM terms and nothing more, which is what makes it a good first example:

- `edges`: the baseline propensity to cosponsor (the intercept).
- `absdiff("ideol")`: ideological homophily: are cosponsorships more likely between senators with similar DW-NOMINATE scores?
- `edgecov("dist")`: geographic distance: does physical distance between home states suppress cosponsorship?

::: {.callout-note collapse="true" title="What is DW-NOMINATE, for non-Americanists"}
DW-NOMINATE is a scaling of roll-call voting that places each legislator on a liberal–conservative dimension, roughly −1 (most liberal) to +1 (most conservative). It is the standard measure of ideology in the US-politics literature. For our purposes it is just a continuous nodal covariate; nothing about the method depends on what it measures.
:::

The data ship in `networkdata` as `legnet`, which provides three objects: `net` (the adjacency matrix), `dwnom` (ideology scores), and `edist` (a matrix of geographic distances).

```{r}
#| label: senate-data
suppressMessages({
  library(ergm)
  library(network)
  library(networkdata)
})

data(legnet)
c(nodes = nrow(net), edges = sum(net), dwnom_rows = nrow(dwnom),
  labels_aligned = identical(rownames(net), dwnom$labs))
```

That last check: `labels_aligned`: is not decorative. The single most common way to ruin a network analysis is to attach a covariate vector whose order does not match the network's nodes, which silently assigns the wrong ideology to every senator. Confirm alignment before you attach anything.

::: {.callout-warning title="A trap the original version of this material fell into"}
An earlier version of this analysis built the network from an object called `el` (an edge list) that **does not exist** in `legnet`. With no `el` in scope, R quietly resolved the bare name to `methods::el()`: a function, and the whole downstream analysis cascaded from there.

The fix, which is what we do below, is to build directly from the adjacency matrix `net` that `legnet` actually provides. The general lesson: if a data object seems to appear from nowhere, check it is not a function of the same name. `exists("el", where = globalenv())` would have caught it.
:::

## 1.2 Build the network and attach covariates

```{r}
#| label: senate-build
senate <- network(as.matrix(net), directed = TRUE)

# Ideology as a vertex attribute; labels already confirmed aligned above
network::set.vertex.attribute(senate, "ideol", dwnom$dwnom)

# Geographic distance as a network (dyadic) attribute
set.network.attribute(senate, "dist", as.matrix(edist))

senate
```

## 1.3 Fit and interpret

```{r}
#| label: senate-fit
m_sen <- fit_or_load("app_senate",
  ergm(senate ~ edges + absdiff("ideol") + edgecov("dist"),
       control = control.ergm(seed = 6886)))
round(coef(m_sen), 4)
is.dyad.independent(m_sen)
```

`is.dyad.independent()` returns `TRUE`, so this fit used no MCMC at all: it is exact logistic regression on the dyads, which is why it finished instantly and needs no convergence diagnostics. Good to know before you go looking for traceplots that do not exist.

Now read the coefficients, remembering the rule from the taught session: each coefficient is a contribution to the **conditional log-odds** of a tie.

- **`absdiff("ideol")` is strongly negative.** `absdiff` grows as two senators become *more* different, so a negative coefficient means difference *suppresses* cosponsorship, which is ideological **homophily**. Similar senators cosponsor.
- **`edgecov("dist")` is negative**, so greater geographic distance lowers the odds of cosponsorship, though the effect is small.

Convert the homophily coefficient to something interpretable. Two senators one full DW-NOMINATE unit apart (about the width of the whole scale: a hard-left and a hard-right senator):

```{r}
#| label: senate-interp
b <- coef(m_sen)
# Odds ratio for a one-unit increase in ideological distance, holding geography fixed
exp(b[["absdiff.ideol"]])
```

An odds ratio around 0.04: moving from identical ideology to opposite ends of the scale cuts the odds of cosponsorship to a small fraction. In this chamber, ideology is doing most of the work, which for the most extreme senators of a polarised Congress is exactly what you would expect. That the model recovers it cleanly is the reassurance you want before trusting it on a harder question.

# 2. Rebel-group cooperation, reanalyzed with an ERGM

Gade, Gabbay, Hafez & Kelly (2019) ask why rebel groups in fragmented civil wars cooperate. Their answer, tested on Syrian rebel groups: **ideological similarity** is the primary driver, above shared sponsors or raw power. Their published regression is an additive and multiplicative effects model of a square-root-transformed count outcome, with raw-count and ordinal checks in the supplement. It is not an ERGM. We use their data for a binary ERGM teaching exercise, so the results below are a reanalysis, not a replication of Gade et al.

## 2.1 From edge list to network

The data arrive as `gadeData`, one row per ordered pair of groups, with dyadic and nodal variables mixed together in the columns. Turning that into a network plus its covariates is most of the work, and it is worth doing carefully because it is exactly the shape most real relational data comes in.

```{r}
#| label: gade-data
data(gadeData)
dim(gadeData)
names(gadeData)
```

The dependent variable, `coopActions`, is a count of cooperative acts. We binarize it as any cooperation versus none so that the example stays within the binary ERGM material taught on Day 11. This changes the estimand and discards intensity, so it cannot reproduce the article's count model.

```{r}
#| label: gade-dv
gadeData$coopBin <- as.numeric(gadeData$coopActions > 0)
table(gadeData$coopBin)
```

The reshaping. We have two kinds of covariate and they go to different places:

- **dyadic** variables (`ideol_diff.dyad`, `powerdiff.dyad`, `loc.dyad`, `spons.dyad`) become **matrices**, one value per pair;
- **nodal** variables (`averageId.node`, `size.node`, `spons_actor.node`) become **vertex attributes**, one value per group.

```{r}
#| label: gade-reshape
actors   <- sort(unique(c(gadeData$Var1, gadeData$Var2)))
n        <- length(actors)
dyadVars <- c("coopBin", "ideol_diff.dyad", "powerdiff.dyad", "loc.dyad", "spons.dyad")

# Fill an actor x actor x variable array by NAME, never by position
arr <- array(0, dim = c(n, n, length(dyadVars)),
             dimnames = list(actors, actors, dyadVars))
for (v in dyadVars) {
  for (i in seq_len(nrow(gadeData))) {
    arr[gadeData$Var1[i], gadeData$Var2[i], v] <- gadeData[i, v]
  }
}

nodeVars <- c("averageId.node", "size.node", "spons_actor.node")
nodeData <- unique(gadeData[, c("Var1", nodeVars)])
rownames(nodeData) <- nodeData$Var1
nodeData <- nodeData[actors, nodeVars]      # reorder to match `actors` exactly

c(n_actors = n, cooperating_pairs = sum(arr[, , "coopBin"]) / 2)
```

::: {.callout-tip collapse="true" title="Stuck? Why fill the array by name and not by position"}
The loop indexes with `arr[gadeData$Var1[i], gadeData$Var2[i], v]`: using the group *names* as indices into a named array. It does **not** loop over positions `1:n`.

This matters enormously. If you build the array positionally, you are trusting that the order of rows in `gadeData` matches the order of `actors`, and it does not: `actors` is sorted, the data frame is not. Positional filling would misalign every pair. Name-based indexing cannot misalign, because R looks each name up. Whenever you reshape relational data, index by name and let R do the matching. The companion [Duque SAOM document](saom_duque.html) has a longer horror story about exactly this.
:::

```{r}
#| label: gade-net
rebels <- as.network(arr[, , "coopBin"], directed = FALSE, loops = FALSE,
                     matrix.type = "adjacency")
for (v in nodeVars) network::set.vertex.attribute(rebels, v, nodeData[, v])
set.network.attribute(rebels, "loc.dyad",   arr[, , "loc.dyad"])
set.network.attribute(rebels, "spons.dyad", arr[, , "spons.dyad"])

checkpoint(groups = network.size(rebels),
           coop_ties = network.edgecount(rebels))
```

## 2.2 A dyad-independent model first

Gade et al.'s core hypotheses are all about covariates: ideological difference, power difference, shared location, shared sponsor. Start with exactly those and no structural terms.

```{r}
#| label: gade-m0
m0 <- fit_or_load("app_gade0",
  ergm(rebels ~ edges +
         nodecov("averageId.node") + nodecov("size.node") +
         absdiff("averageId.node") + absdiff("size.node") +
         edgecov("loc.dyad") + edgecov("spons.dyad"),
       control = control.ergm(seed = 6886)))
round(coef(m0), 3)
is.dyad.independent(m0)
```

Dyad-independent, so again no MCMC. The headline is `absdiff("averageId.node")`: **negative**, meaning groups that are *further apart* ideologically cooperate *less*. That is Gade et al.'s central qualitative finding, ideological proximity goes with more cooperation, echoed in one term under our different binary specification. It is not a reproduction of their coefficient. `edgecov("loc.dyad")` is strongly positive: groups operating in the same location cooperate much more, which makes sense and which you would want to control for before crediting ideology.

## 2.3 Adding endogenous structure, and a real choice

Here is where it stops being a logit. Cooperation among rebel groups plausibly has a **triadic** logic: if A cooperates with B and B with C, A and C have both a broker and a reason. That motivates a shared-partner term, which means MCMC and a decision about how to specify it. Whether that term represents closure rather than omitted actor activity is a separate question.

The naive choice is `triangle`, and the taught session showed you why that ends in degeneracy. The right choice is `gwesp`, geometrically weighted so that the tenth shared partner counts for less than the first.

```{r}
#| label: gade-m1
m1 <- fit_or_load("app_gade1",
  ergm(rebels ~ edges +
         nodecov("averageId.node") + nodecov("size.node") +
         absdiff("averageId.node") + absdiff("size.node") +
         edgecov("loc.dyad") + edgecov("spons.dyad") +
         gwesp(0.5, fixed = TRUE),
       control = control.ergm(seed = 6886)))
round(coef(m1), 3)
is.dyad.independent(m1)
```

::: {.callout-important title="`gwesp(0.5, fixed = TRUE)`: the `fixed` is load-bearing"}
Write `gwesp(decay = 0.5)` without `fixed = TRUE` and, as of `ergm` 4, the decay value is **ignored** and you silently fit a many-parameter *curved* model instead of the one-parameter model you intended. The original version of this material had `gwesp(decay = 0.25)`: unfixed: in five separate places, every one of which was fitting a curved model nobody meant to fit.

The rule from the taught session, restated because it is the single most common error in inherited ERGM code: **`fixed = TRUE` unless you are deliberately estimating a curved model.**
:::

`gwesp` comes out **positive**: conditional on this specification, ties with more shared partners have higher conditional log-odds. That is not yet evidence of a closure mechanism because hubs can mechanically create shared partners, and this model has no actor-specific activity term. And notice what happens to the covariate coefficients when you add it: they shift, because some of what looked like covariate effect in `m0` was unmodelled structure. That shift is the reason you cannot stop at the dyad-independent model if you think structure is present.

::: {.callout-warning title="But do not now compare m0 and m1 coefficients as effect sizes"}
It is tempting to say "ideology's coefficient dropped from −0.24 to −0.23 when we added closure, so most of the ideology effect is real." Resist stating it that precisely. Per Duxbury (2023), adding a term to an ERGM rescales all the others even when they are uncorrelated, so part of any across-model shift is rescaling rather than substance. The *direction* is informative; the exact magnitude of the change is not. If you need the defensible version, compare average marginal effects with `ergMargins`, not coefficients.
:::

## 2.4 Check the fit

A model with an endogenous term has to be checked: that is the whole discipline of the taught session's block C. Look at whether the model reproduces network features it was **not** fitted on.

```{r}
#| label: gade-gof
#| fig-width: 8
#| fig-height: 5.5
set.seed(6886)
gof_m1 <- fit_or_load("app_gade1_gof",
  gof(m1, GOF = ~ degree + espartners + distance - model))
par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
plot(gof_m1)
```

The observed values (the thick line) should sit inside the simulated boxplots. Where they do, the model reproduces that feature; where they stray outside, it does not. For a network this small, expect some ragged tails. Here the default panels broadly bracket the observed line, but that is compatibility with this chosen suite, not proof of fit or proof that `gwesp` identifies closure. Day 14's common-yardstick analysis adds actor heterogeneity and revisits this exact interpretation.

```{r}
#| label: gade-gof-tab
#| results: hide
# The tabular, in-model version is faster and more legible than the plots
gof_model <- fit_or_load("app_gade1_gofmodel", gof(m1, GOF = ~ model))
```

```{r}
#| label: gade-gof-tab-show
gof_model
```

The `GOF = ~model` table checks the in-model statistics directly: each observed statistic against the distribution the fitted model simulates. P-values near the middle of the range mean the model reproduces its own targets, which is the minimum you should demand before interpreting anything.

# 3. Exercises

::: {.panel-tabset}

## Core

Add `nodematch` on shared sponsorship to the Senate model: except the Senate data has no sponsorship variable, so instead add the **reciprocity** term `mutual` to the Senate model and refit.

```{r}
#| eval: false
m_sen2 <- ergm(senate ~ edges + absdiff("ideol") + edgecov("dist") + ______,
               control = control.ergm(seed = 6886))
```

Then answer: `is.dyad.independent()` was `TRUE` for the original Senate model. Predict what it returns now, and what that means for whether the fit uses MCMC.

## Stretch

In the Gade model, `absdiff("averageId.node")` tests ideological homophily as a *distance*. But Gade et al. theorise about ideological *positions*, not just distances: extremist groups might behave differently from moderate ones regardless of whom they are paired with.

Add `nodecov("averageId.node")` (already in the model) *and* think about whether an interaction or a squared term would test the "extremists behave differently" hypothesis. Fit something that tests it, and report whether the data support it. This is open-ended; the point is to translate a substantive hypothesis into a term, which is the actual skill.

## Solution

```{r}
#| label: ex-senate
m_sen2 <- fit_or_load("app_senate_mutual",
  ergm(senate ~ edges + absdiff("ideol") + edgecov("dist") + mutual,
       control = control.ergm(seed = 6886)))
round(coef(m_sen2), 4)
is.dyad.independent(m_sen2)
```

**Core.** `is.dyad.independent()` now returns `FALSE`, and the model fits by MCMC because `mutual` makes the probability of the $i \to j$ tie depend on the $j \to i$ tie, which is precisely dyadic dependence. The moment you add it, you leave exact logistic regression behind: the fit becomes stochastic, it can in principle fail to converge, and you now owe the reader `mcmc.diagnostics()`. That single term is the boundary between the two halves of the taught session.

Substantively, `mutual` comes out positive: cosponsorship is reciprocated, as you would expect , and the ideology and distance coefficients barely move, because in this small chamber reciprocity and homophily are not competing to explain the same ties.

```{r}
#| label: ex-gade
# One way to test "extremists behave differently": a squared ideology term.
# averageId.node is a group's mean ideology; its square is largest at the extremes.
gadeData2 <- gadeData
# Build the squared nodal covariate on the network
id2 <- nodeData[, "averageId.node"]^2
network::set.vertex.attribute(rebels, "id_sq", id2)

m_gade2 <- fit_or_load("app_gade_sq",
  ergm(rebels ~ edges +
         nodecov("averageId.node") + nodecov("id_sq") + nodecov("size.node") +
         absdiff("averageId.node") + absdiff("size.node") +
         edgecov("loc.dyad") + edgecov("spons.dyad") +
         gwesp(0.5, fixed = TRUE),
       control = control.ergm(seed = 6886)))
round(coef(m_gade2), 3)
```

**Stretch.** Adding `nodecov("id_sq")`: the square of a group's mean ideology, which is largest for the most extreme groups in either direction: tests whether extremists have a different baseline propensity to cooperate, separately from whom they are paired with.

Read the `nodecov.id_sq` coefficient against its standard error (`summary(m_gade2)` for the full table). If it is indistinguishable from zero, the data do not support "extremists behave differently in level" once you have already accounted for ideological *distance* via `absdiff`: which would say only that this binary ERGM does not need the squared nodal term once pairwise ideological distance is included. It would not validate the article's model by itself.

The methodological point is the transferable one: a substantive hypothesis ("extremists are different") became a specific term (`nodecov` on a squared covariate), and the term is testable. Half of applied ERGM work is exactly this translation, and it is a skill worth more than any particular result.
:::

# 4. Readings

- **Cranmer & Desmarais (2011)**, "Inferential network analysis with ERGMs," *Political Analysis* 19(1):66–86. The methodological case for using ERGMs in political science. Assigned for the taught session.
- **Gade, Gabbay, Hafez & Kelly (2019)**, "Networks of cooperation: Rebel alliances in fragmented civil wars," *Journal of Conflict Resolution*. The data and substantive argument behind the Section 2 reanalysis. The published regression is AME on a transformed count outcome, not the binary ERGM fitted here. Note this is a **day 10 (latent geometry, W2 Fri)** reading in 2026: it belongs to another session, and the latent-factor treatment of the same data is the natural companion to this ERGM one.
- **Goodreau, Kitts & Morris (2009)**, "Birds of a feather, or friend of a friend?" *Demography* 46(1):103–125. The `mesa`-style application the taught session is built on, and the template for the "covariates first, structure second" workflow used in §2.

```{r}
#| label: versions
#| cache: false
#| code-fold: true
#| code-summary: "Versions: expand if your numbers differ from mine"
pk <- c("ergm", "network", "networkdata")
data.frame(package = pk,
           version = sapply(pk, function(p) as.character(packageVersion(p))),
           row.names = NULL)
cat("R", as.character(getRversion()), "\n")
```
