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))
#> nodes edges dwnom_rows labels_aligned
#> 9 27 9 1ERGMs in applied political science
Ancillary self-study · ICPSR Network Analysis: Advanced Topics
ImportantOptional: 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.
NoteWhat 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:
- 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.
- 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?
NoteWhat 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).
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.
WarningA 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
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
#> Network attributes:
#> vertices = 9
#> directed = TRUE
#> hyper = FALSE
#> loops = FALSE
#> multiple = FALSE
#> bipartite = FALSE
#> dist: 9x9 matrix
#> total edges= 27
#> missing edges= 0
#> non-missing edges= 27
#>
#> Vertex attribute names:
#> ideol vertex.names
#>
#> No edge attributes1.3 Fit and interpret
m_sen <- fit_or_load("app_senate",
ergm(senate ~ edges + absdiff("ideol") + edgecov("dist"),
control = control.ergm(seed = 6886)))
round(coef(m_sen), 4)
#> edges absdiff.ideol edgecov.dist
#> 1.1349 -3.1330 -0.0167
is.dyad.independent(m_sen)
#> [1] TRUEis.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.absdiffgrows 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):
b <- coef(m_sen)
# Odds ratio for a one-unit increase in ideological distance, holding geography fixed
exp(b[["absdiff.ideol"]])
#> [1] 0.0435869An 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.
data(gadeData)
dim(gadeData)
#> [1] 930 11
names(gadeData)
#> [1] "Var1" "Var2" "coopActions" "id"
#> [5] "ideol_diff.dyad" "powerdiff.dyad" "loc.dyad" "spons.dyad"
#> [9] "averageId.node" "size.node" "spons_actor.node"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.
gadeData$coopBin <- as.numeric(gadeData$coopActions > 0)
table(gadeData$coopBin)
#>
#> 0 1
#> 758 172The 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.
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)
#> n_actors cooperating_pairs
#> 31 86
TipStuck? 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 has a longer horror story about exactly this.
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))
#> ------------------------------------------------------------------
#> CHECKPOINT: groups = 31 | coop_ties = 86
#> ------------------------------------------------------------------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.
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)
#> edges nodecov.averageId.node nodecov.size.node
#> -6.286 0.293 0.127
#> absdiff.averageId.node absdiff.size.node edgecov.loc.dyad
#> -0.241 -0.113 3.168
#> edgecov.spons.dyad
#> -0.129
is.dyad.independent(m0)
#> [1] TRUEDyad-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.
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)
#> edges nodecov.averageId.node nodecov.size.node
#> -7.175 0.202 0.074
#> absdiff.averageId.node absdiff.size.node edgecov.loc.dyad
#> -0.233 -0.081 2.647
#> edgecov.spons.dyad gwesp.fixed.0.5
#> -0.024 1.295
is.dyad.independent(m1)
#> [1] FALSE
Important
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.
WarningBut 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.
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.
# 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))gof_model
#>
#> Goodness-of-fit for model statistics
#>
#> obs min mean max MC p-value
#> edges 86.0000 59.00000 84.5700 112.0000 0.90
#> nodecov.averageId.node 460.5623 319.55720 456.9320 604.5088 1.00
#> nodecov.size.node 1386.1000 1022.90000 1364.9455 1719.8000 0.92
#> absdiff.averageId.node 129.8897 88.44333 126.5079 180.6118 0.84
#> absdiff.size.node 639.9000 413.90000 629.9765 885.3000 0.88
#> edgecov.loc.dyad 85.0000 59.00000 83.5800 110.0000 0.90
#> edgecov.spons.dyad 14.0000 7.00000 13.6000 21.0000 1.00
#> gwesp.fixed.0.5 120.2135 75.07166 117.8774 164.8288 0.88The 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
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.
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.
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.
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)
#> edges absdiff.ideol edgecov.dist mutual
#> 1.2702 -3.1879 -0.0163 -0.2453
is.dyad.independent(m_sen2)
#> [1] FALSECore. 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.
# 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)
#> edges nodecov.averageId.node nodecov.id_sq
#> -10.543 1.376 -0.205
#> nodecov.size.node absdiff.averageId.node absdiff.size.node
#> 0.074 -0.014 -0.079
#> edgecov.loc.dyad edgecov.spons.dyad gwesp.fixed.0.5
#> 3.243 0.080 1.205Stretch. 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.
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)
#> package version
#> 1 ergm 4.12.0
#> 2 network 1.20.0
#> 3 networkdata 0.1
cat("R", as.character(getRversion()), "\n")
#> R 4.3.3