---
title: "Random graphs, and what a null model is for"
engine: knitr
editor: source
subtitle: "Ancillary self-study · ICPSR Network Analysis: Advanced Topics"
---

```{r}
#| label: setup
#| include: false
source("_helpers.R")
require_pkgs(c("igraph", "sna", "ergm", "network"))
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 it is the natural thing to read *before* it if you want a running start.

**The question it answers:** the ERGM session opens by saying that the whole enterprise is about building a better null model. This document is the slow version of that sentence. It builds the simplest possible null: the Erdős–Rényi random graph: by hand, compares it against real data, watches it fail, and shows why the failure is what motivates everything that comes after.

**Why it is optional in 2026:** Olga covers random graphs and the basic descriptives in week 1, so for most of you this is revision. It is here for people who want to build the null by hand once, and for anyone who joined the course with less background than the prerequisites assume.

**What you need:** `igraph`, `sna`, `ergm`, `network`. No data files.

:::

# 1. The question underneath every network paper

You have a network. You compute something about it: 51% of its ties are reciprocated, say, or it has 460 triangles. And then you want to say that number is *interesting*.

Interesting compared to what?

That is the entire problem, and it is not rhetorical. "51% reciprocity" is not a finding until you can say what reciprocity you would have expected if nothing were going on. A network with lots of ties will have lots of reciprocated pairs by sheer arithmetic. So the number on its own tells you nothing; the number *relative to a baseline* might tell you a great deal.

The baseline is called a **null model**, and choosing one is the substantive act. Get in the habit of asking, of any claim that a network is "more clustered than chance" or "more reciprocal than expected": **expected under what?** In a surprisingly large number of published papers the answer is either unstated or indefensible.

# 2. The simplest null: Erdős–Rényi

In an Erdős–Rényi random graph, every possible edge is present with the same probability $p$, independently of every other edge. That independence is the whole content of the model, and it is exactly the assumption you have spent this course learning to distrust.

Which is the point. **It is useful precisely because it is wrong in a known way.** If your data look like an ER graph on some feature, that feature needs no explanation. If they do not, the gap is what you have to explain.

## 2.1 The data

Coleman's high-school friendship network: 73 boys, directed nominations. It ships with `sna`.

```{r}
#| label: data
suppressMessages({
  library(igraph)
  library(sna)
})

data(coleman, package = "sna")
cm <- coleman[1, , ]                          # wave 1
g  <- graph_from_adjacency_matrix(cm, mode = "directed", diag = FALSE)

c(nodes      = vcount(g),
  edges      = ecount(g),
  density    = round(edge_density(g), 4),
  reciprocity = round(reciprocity(g), 4))
```

Before modelling anything, three plain questions, which is the habit worth taking away from this whole document:

1. **How big is it and how dense?** 73 boys, 243 nominations, density about 0.046. Sparse, which is normal: people have a bounded number of friends regardless of how many classmates exist.
2. **Is it reciprocated?** About 51% of nominations are returned. Is that a lot?
3. **Is it clustered?** Do friends of friends tend to be friends?

Questions 2 and 3 are the ones that need a null.

```{r}
#| label: plot
#| fig-width: 6
#| fig-height: 5
set.seed(6886)
plot(g, vertex.label = NA, vertex.color = "steelblue", vertex.size = 5,
     edge.arrow.size = 0.2, main = "Coleman friendship, wave 1")
```

## 2.2 Building the null by hand

`igraph` has `sample_gnp()` and you should use it in real work. Build one by hand first, once, because the mechanics are the argument.

$p$ is the observed number of edges over the number of possible edges. With 73 nodes and a directed graph and no self-loops, that is $73 \times 72 = 5256$ possible edges.

```{r}
#| label: build-null
n <- vcount(g)
p <- ecount(g) / (n * (n - 1))
round(p, 4)

set.seed(6886)
rnet <- matrix(rbinom(n * n, 1, p), nrow = n, ncol = n)
diag(rnet) <- 0                                # no self-loops

rg <- graph_from_adjacency_matrix(rnet, mode = "directed", diag = FALSE)
c(edges = ecount(rg), density = round(edge_density(rg), 4))
```

::: {.callout-tip collapse="true" title="Stuck? Why `rbinom` rather than a double loop"}
The classic way to write this is nested `for` loops with `runif(1) < p` inside. It works and it is transparent, and it is also about a thousand times slower, which stops being funny at 500 nodes.

`rbinom(n*n, 1, p)` draws all $n^2$ coin flips in one vectorised call and `matrix()` reshapes them. Same model, same distribution, no loops. This is the single most useful vectorisation habit in R.

One subtlety worth noticing: `diag(rnet) <- 0` sets $n$ cells to zero *after* drawing them, so the realised density comes out very slightly below $p$. With $n = 73$ that is a rounding-level difference. If it ever matters, draw only the off-diagonal cells.
:::

## 2.3 Side by side

```{r}
#| label: compare-plot
#| fig-width: 9
#| fig-height: 4.5
set.seed(6886)
par(mfrow = c(1, 2), mar = c(1, 1, 3, 1))
plot(g,  vertex.label = NA, vertex.color = "steelblue", vertex.size = 4,
     edge.arrow.size = 0.15, main = "observed")
plot(rg, vertex.label = NA, vertex.color = "grey60",  vertex.size = 4,
     edge.arrow.size = 0.15, main = "Erdős–Rényi, same density")
```

Same number of nodes, essentially the same number of edges. They do not look alike. The real network has visible clumps and a scattering of near-isolates; the random one is uniform mush.

Eyeballing is not evidence, though, so put numbers on it.

```{r}
#| label: compare-stats
suppressMessages({
  library(ergm)
  library(network)
})

obs_net <- network(cm,   directed = TRUE)
rnd_net <- network(rnet, directed = TRUE)

rbind(observed = summary(obs_net ~ edges + mutual + triangle + istar(2) + ostar(2) + idegree(6)),
      random   = summary(rnd_net ~ edges + mutual + triangle + istar(2) + ostar(2) + idegree(6)))
```

Roughly matched on `edges`, because we built it that way. Everything else is wildly off. Real friendships are reciprocated and clustered; random ones are not.

# 3. One draw is not a null distribution

Everything above compares the data against **one** random graph. That is a demonstration, not a test: a single draw could be unusual. The real procedure is to simulate many, build the sampling distribution of your statistic, and see where the observed value falls.

```{r}
#| label: simulate
#| fig-width: 8
#| fig-height: 3.4
set.seed(6886)
B <- 1000

sim_stats <- t(replicate(B, {
  r <- matrix(rbinom(n * n, 1, p), n, n)
  diag(r) <- 0
  rn <- network(r, directed = TRUE)
  summary(rn ~ mutual + triangle + istar(2))
}))

obs_stats <- summary(obs_net ~ mutual + triangle + istar(2))

par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))
for (j in 1:3) {
  hist(sim_stats[, j], breaks = 30, col = "grey85", border = "white",
       xlim = range(c(sim_stats[, j], obs_stats[j])),
       main = colnames(sim_stats)[j], xlab = "value under ER null")
  abline(v = obs_stats[j], col = "firebrick", lwd = 3)
}
```

```{r}
#| label: sim-table
data.frame(
  statistic = names(obs_stats),
  observed  = as.numeric(obs_stats),
  null_mean = round(colMeans(sim_stats), 1),
  null_max  = apply(sim_stats, 2, max),
  ratio     = round(as.numeric(obs_stats) / colMeans(sim_stats), 2),
  row.names = NULL)
```

The red line is not merely in the tail. **It is off the chart.** In a thousand draws the null never comes close. The observed reciprocity and triangle counts are multiples of what independence produces.

```{r}
#| label: checkpoint
checkpoint(mutual_obs   = obs_stats[["mutual"]],
           mutual_null  = round(mean(sim_stats[, "mutual"]), 1),
           triangle_obs = obs_stats[["triangle"]],
           triangle_null = round(mean(sim_stats[, "triangle"]), 1))
```

So we can say something definite: **friendship nominations in this school are not independent.** Whatever is generating this network, it is not coin flips.

# 4. Where this stops being enough, and the point of the whole exercise

Here is the part that matters, and it is the reason this document exists rather than just the histogram above.

We have shown the data are not Erdős–Rényi. Now try to say something more specific, and watch the problem appear.

**Claim:** "this network exhibits triadic closure: friends of friends become friends."

The evidence offered: 460 triangles observed against about 49 expected under the ER null. A factor of nine. Sounds convincing.

**The objection:** the ER null holds *nothing* fixed except density. In particular it does not preserve the degree distribution, and the real network has some very popular boys. A node with many incoming nominations mechanically generates a large number of two-paths, and two-paths close into triangles at some rate for purely combinatorial reasons. So some of that gap may not be closure at all. It may be popularity.

Notice that this is a **hypothesis**, not a conclusion. People state it as though stating it settles it. It does not: you have to check. Swap the null for one that **preserves each node's in- and out-degree exactly** and randomises only who is connected to whom.

```{r}
#| label: rewire
set.seed(6886)
B2 <- 500

# Rewire the observed simple graph directly. Every draw preserves each node's
# in- and out-degree exactly, does not introduce loops, and remains a simple graph.
din  <- igraph::degree(g, mode = "in")
dout <- igraph::degree(g, mode = "out")
rewired <- t(replicate(B2, {
  gg <- igraph::rewire(
    g,
    with = igraph::keeping_degseq(loops = FALSE, niter = 20 * igraph::ecount(g))
  )
  stopifnot(identical(igraph::degree(gg, mode = "in"), din),
            identical(igraph::degree(gg, mode = "out"), dout))
  nn <- network(as.matrix(igraph::as_adjacency_matrix(gg, sparse = FALSE)),
                directed = TRUE)
  summary(nn ~ mutual + triangle)
}))

data.frame(
  statistic     = c("mutual", "triangle"),
  observed      = as.numeric(obs_stats[c("mutual", "triangle")]),
  ER_null       = round(colMeans(sim_stats[, c("mutual", "triangle")]), 1),
  ratio_vs_ER   = round(as.numeric(obs_stats[c("mutual", "triangle")]) /
                          colMeans(sim_stats[, c("mutual", "triangle")]), 2),
  degree_null   = round(colMeans(rewired), 1),
  ratio_vs_deg  = round(as.numeric(obs_stats[c("mutual", "triangle")]) /
                          colMeans(rewired), 2),
  row.names     = NULL)
```

Now look at the two ratio columns for `triangle`, and notice that **the answer is not the dramatic one.** The excess-over-random does fall when you hold degree fixed: the direction is exactly as the objection predicted, but it falls by about one third, not by an order of magnitude. Coleman's boys really are clustering, and the degree distribution explains only a modest part of it.

That is a genuinely useful outcome, and I want to be clear that I did not arrange it. The objection in the previous paragraph is a good objection and it is often decisive. Here it turns out not to be, because this network's degree distribution is not skewed enough for the combinatorial mechanism to do much work.

```{r}
#| label: degdist
#| fig-width: 7
#| fig-height: 3.2
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
hist(igraph::degree(g, mode = "in"),  breaks = 12, col = "grey85", border = "white",
     main = "in-degree", xlab = "nominations received")
hist(igraph::degree(g, mode = "out"), breaks = 12, col = "grey85", border = "white",
     main = "out-degree", xlab = "nominations made")
```

Out-degree is nearly fixed by the survey instrument: boys were asked for a bounded number of friends, and in-degree, while it has a tail, has no hubs of the kind that would generate triangles wholesale. On a network with a genuinely heavy-tailed degree distribution, a citation network or a follower graph, the same comparison routinely wipes out most of the apparent clustering.

::: {.callout-warning title="A caveat on the `mutual` row"}
The `mutual` comparison in that table is **not** a fair test and you should not read it as one. The configuration model does not preserve reciprocity, so a degree-preserving null still expects almost no mutual dyads: the ratio stays high for a mechanical reason, not a substantive one.

If you actually wanted to test reciprocity against a null that holds degree fixed, you would need a null that also fixes the number of mutual dyads, which is `sna::rguman()` territory, or you would fit an ERGM with `mutual` in it and look at the other terms. Included here because running an inappropriate comparison and then noticing is more instructive than quietly omitting it.
:::

::: {.callout-important title="The takeaway, which is the whole document in three sentences"}
The choice of null model does the real work in any "more X than chance" claim.

Two researchers can look at the same network, compute the same statistic, and reach different conclusions purely by choosing different baselines, and neither has necessarily done anything wrong. What they have done is answer different questions while writing the same sentence.

So whenever you read or write such a claim, the question to ask is not "is the number big?" but **"what was held fixed?"**, and then, because the answer is not always the one you expect, actually run the alternative.
:::

# 5. Which is why you need ERGMs

You could keep going like this: build an ever more elaborate null by hand, preserving degree, then preserving degree and attribute mixing, then those plus something else. People did exactly that for years.

The problem is that it does not scale and it does not let you ask the question you actually care about. You do not want to know whether triangles are in excess of *some* baseline. You want to know **how much of the clustering is triadic closure once you have accounted for homophily, popularity, and reciprocity simultaneously**, and, ideally, with a standard error attached.

That is a regression question. And it is exactly the question the ERGM answers, by making the null a **fitted model** with parameters rather than a fixed baseline you chose in advance:

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

Every hand-built null above is a special case. Erdős–Rényi is `~ edges`: one parameter, and the maximum-likelihood fit reproduces the observed density, which is precisely how we chose $p$ in §2.2. Adding `mutual` gives a null that also preserves reciprocity. Adding `gwesp` lets closure be a parameter you *estimate* rather than a baseline you assume.

```{r}
#| label: ergm-er
#| output: false
m_er <- fit_or_load("rg_er",
  ergm(obs_net ~ edges, control = control.ergm(seed = 6886)))
```

```{r}
#| label: ergm-er-show
c(coef = round(coef(m_er)[["edges"]], 4),
  implied_p = round(plogis(coef(m_er)[["edges"]]), 4),
  hand_p = round(p, 4))
```

The one-parameter ERGM recovers the $p$ we computed by hand in §2.2, to four decimals. That is not a coincidence and it is worth sitting with: **the Erdős–Rényi model is an ERGM with one term.** Everything the ERGM session does from there is adding terms to this.

# 6. Exercises

::: {.panel-tabset}

## Core

Run the whole §2–3 pipeline on the **Sampson** monastery data instead of Coleman: build an ER null at matched density, simulate 1000 draws, and compare `mutual` and `triangle`.

```{r}
#| eval: false
data(sampson, package = "ergm")
sm <- as.matrix(samplike)
n_s <- nrow(sm)
p_s <- sum(sm) / (n_s * (n_s - ___))          # careful with this
```

Then: Sampson's density is about 0.29, roughly six times Coleman's. Before you run it, predict whether the observed-to-null **ratio** for triangles will be larger or smaller than Coleman's, and say why. Then check.

## Stretch

Coleman ships with **two** waves (`coleman[1,,]` and `coleman[2,,]`). Compute the reciprocity and triangle excess-over-ER for both, and then ask a question the ER null cannot answer: is the change between waves larger than you would expect from sampling variability alone?

Sketch how you would build a null for *that* question. This is genuinely harder than anything above, and it is the question the [Temporal ERGMs](tergm.html) document exists to answer properly.

## Solution

```{r}
#| label: ex-sampson
data(sampson, package = "ergm")
sm  <- as.matrix(samplike)
n_s <- nrow(sm)
p_s <- sum(sm) / (n_s * (n_s - 1))            # n(n-1) ordered pairs, no self-loops

set.seed(6886)
sim_s <- t(replicate(1000, {
  r <- matrix(rbinom(n_s * n_s, 1, p_s), n_s, n_s)
  diag(r) <- 0
  summary(network(r, directed = TRUE) ~ mutual + triangle)
}))
obs_s <- summary(samplike ~ mutual + triangle)

data.frame(
  network   = rep(c("Coleman", "Sampson"), each = 2),
  density   = round(rep(c(p, p_s), each = 2), 3),
  statistic = rep(c("mutual", "triangle"), 2),
  observed  = c(as.numeric(obs_stats[c("mutual", "triangle")]), as.numeric(obs_s)),
  null_mean = round(c(colMeans(sim_stats[, c("mutual", "triangle")]), colMeans(sim_s)), 1),
  ratio     = round(c(as.numeric(obs_stats[c("mutual", "triangle")]) /
                        colMeans(sim_stats[, c("mutual", "triangle")]),
                      as.numeric(obs_s) / colMeans(sim_s)), 2),
  row.names = NULL)
```

**Core, the blank:** `n_s - 1`. There are $n(n-1)$ ordered pairs of distinct nodes in a directed network. Using $n(n-1)/2$ is the undirected count and would double your $p$; using $n^2$ counts self-loops that cannot exist. Both mistakes are common and both are silent.

**Core, the prediction.** The triangle ratio is **much smaller** for Sampson, and the reason is worth understanding because it generalises.

Under an ER null the expected number of triangles scales roughly like $n^3 p^3$. Sampson's density is about six times Coleman's, so its *expected* triangle count is enormous: cubing a much larger $p$ dominates the smaller $n$. The null is already producing lots of triangles by chance, so there is far less room for the observed count to exceed it.

Coleman's ratio is astronomical mainly because the denominator is nearly zero: at density 0.046, random graphs on 73 nodes almost never make triangles, so *any* clustering looks miraculous.

**The general lesson, and it is a warning:** excess-over-random ratios are not comparable across networks of different density. A paper reporting "40× more clustered than random" for a sparse network and one reporting "2× more clustered than random" for a dense one may be describing identical amounts of social closure. If you need a comparable quantity, fit a model with a density term in it and compare the *closure parameter*. This is, once again, the argument for ERGMs.

**Stretch.** The honest sketch: you would need a null that preserves whatever you are willing to treat as fixed across waves: node set, and probably each wave's density, while randomising the *association between* waves. A permutation approach gets you started: repeatedly relabel nodes in wave 2, recompute the change in your statistic, and build a null distribution for "change under no wave-to-wave dependence."

That answers a narrow version of the question. It does not touch the interesting one, which is whether the *process* generating change is reciprocity-driven or closure-driven, because a permutation null has no process in it at all. For that you need a model of the transition itself, which is a TERGM or a SAOM. That is exactly why those models exist.
:::

# 7. Readings

- **Erdős & Rényi (1959)**, "On random graphs I," *Publicationes Mathematicae* 6:290–297. The original. Short, and more readable than you would expect.
- **Cranmer & Desmarais (2011)**, "Inferential network analysis with exponential random graph models," *Political Analysis* 19(1):66–86. Assigned for the ERGM session. Section 2 is the argument of §5 above, made properly.
- **Milo et al. (2002)**, "Network motifs: simple building blocks of complex networks," *Science* 298:824–827. The canonical demonstration that degree-preserving nulls and ER nulls give different answers, and that the difference is the finding.

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