Random graphs, and what a null model is for

Ancillary self-study · ICPSR Network Analysis: Advanced Topics

Author

Shahryar Minhas

Published

July 22, 2026

Optional: not covered in the 2026 course

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

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

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

What this extends, and why it is here

Extends: the ERGM session (day 11, W3 Mon), and 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.

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))
#>       nodes       edges     density reciprocity 
#>     73.0000    243.0000      0.0462      0.5103

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

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.
Is it reciprocated? About 51% of nominations are returned. Is that a lot?
Is it clustered? Do friends of friends tend to be friends?

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

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.

n <- vcount(g)
p <- ecount(g) / (n * (n - 1))
round(p, 4)
#> [1] 0.0462

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))
#>    edges  density 
#> 255.0000   0.0485

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

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.

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)))
#>          edges mutual triangle istar2 ostar2 idegree6
#> observed   243     62      460    542    383        5
#> random     255      4       60    495    421        7

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.

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

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)
#>   statistic observed null_mean null_max ratio
#> 1    mutual       62       5.6       16 11.02
#> 2  triangle      460      49.2       97  9.36
#> 3    istar2      542     398.0      624  1.36

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.

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))
#> ------------------------------------------------------------------
#> CHECKPOINT: mutual_obs = 62   |   mutual_null = 5.6   |   triangle_obs = 460   |   triangle_null = 49.2
#> ------------------------------------------------------------------

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.

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)
#>   statistic observed ER_null ratio_vs_ER degree_null ratio_vs_deg
#> 1    mutual       62     5.6       11.02         7.2         8.64
#> 2  triangle      460    49.2        9.36        73.4         6.27

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.

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.

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.

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.

m_er <- fit_or_load("rg_er",
  ergm(obs_net ~ edges, control = control.ergm(seed = 6886)))

c(coef = round(coef(m_er)[["edges"]], 4),
  implied_p = round(plogis(coef(m_er)[["edges"]]), 4),
  hand_p = round(p, 4))
#>      coef implied_p    hand_p 
#>   -3.0267    0.0462    0.0462

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

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.

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.

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 document exists to answer properly.

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)
#>   network density statistic observed null_mean ratio
#> 1 Coleman   0.046    mutual       62       5.6 11.02
#> 2 Coleman   0.046  triangle      460      49.2  9.36
#> 3 Sampson   0.288    mutual       28      12.8  2.19
#> 4 Sampson   0.288  triangle      193     156.3  1.24

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.

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)
#>   package version
#> 1  igraph   2.2.2
#> 2     sna     2.8
#> 3    ergm  4.12.0
#> 4 network  1.20.0
cat("R", as.character(getRversion()), "\n")
#> R 4.3.3