Checking a network analysis

R
ecological networks
species interactions
ggplot2
ecology tutorial
How far sampling effort, network size and unseen links move connectance, H2’ and nestedness, what a null model rescues, and what the matrix can never tell you.
Author

Tidy Ecology

Published

2026-06-10

Modified

2026-09-30

Corrected 30 September 2026: the null-model section said the web has no nestedness beyond its degree sequence, that the one significant z-score was unremarkable, and that the z-score does not move with effort; the trait-matching web is less nested than its degrees imply, repeated samples show the z-score falling as effort grows, and the post now says the z-score does not track the raw index but still depends on effort through power. The size-dependence section also read the simulated fall in connectance as the biological scaling of connectance with richness; it now says the fall here is a sampling effect, because the same effort per species is spread over more cells.

The three earlier posts in this series built the toolkit: metrics for describing a web, modularity for its compartments, and null models for deciding whether any of it beats chance. Each ended on the same unfinished thought: the numbers move when the sampling moves. This post makes that concrete by building a web whose true interaction probabilities are known, then sampling it at different intensities and watching what happens. Some of the damage a null model repairs. One part of it nothing repairs, and that is the honest limit worth carrying into your own data.

Setup: packages, colours and plot theme
library(ggplot2)
library(grid)
te_ink <- "#16241d"; te_body <- "#2c3a31"; te_forest <- "#275139"; te_label <- "#46604a"
te_sage <- "#93a87f"; te_paper <- "#f5f4ee"; te_line <- "#dad9ca"; te_gold <- "#c9b458"
te_rust <- "#b5534e"
theme_te <- function(base = 12) theme_minimal(base_size = base) +
  theme(plot.background = element_rect(fill = te_paper, colour = NA),
        panel.background = element_rect(fill = te_paper, colour = NA),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_line(colour = te_line, linewidth = 0.3),
        axis.text = element_text(colour = te_body), axis.title = element_text(colour = te_label),
        plot.title = element_text(colour = te_ink, face = "bold", size = rel(1.02)),
        plot.subtitle = element_text(colour = te_label, size = rel(0.9)),
        legend.text = element_text(colour = te_body), legend.title = element_text(colour = te_label))
# H2' and NODF as coded in the earlier posts of this series
H2prime <- function(W) {
  m <- sum(W); if (m == 0) return(NA_real_)
  pm <- W[W > 0] / m; H2 <- -sum(pm * log(pm))
  r <- rowSums(W) / m; cc <- colSums(W) / m
  H2max <- -sum(r[r > 0] * log(r[r > 0])) - sum(cc[cc > 0] * log(cc[cc > 0]))
  rr <- rowSums(W); ccx <- colSums(W); ent <- 0
  repeat {
    if (sum(rr) <= 0) break
    best <- -1; bi <- bj <- NA
    for (i in which(rr > 0)) for (j in which(ccx > 0)) { v <- min(rr[i], ccx[j]); if (v > best) { best <- v; bi <- i; bj <- j } }
    p <- best / m; ent <- ent - p * log(p); rr[bi] <- rr[bi] - best; ccx[bj] <- ccx[bj] - best
  }
  if (H2max <= ent) return(0)
  (H2max - H2) / (H2max - ent)
}
nodf <- function(M) {
  pc <- function(deg, mat, mar) {
    idx <- order(deg, decreasing = TRUE); s <- 0; np <- 0
    for (a in seq_along(idx)) for (b in seq_along(idx)) if (a < b) {
      i <- idx[a]; j <- idx[b]; di <- deg[i]; dj <- deg[j]
      if (di > dj && dj > 0) { ov <- if (mar == "r") sum(mat[i, ] & mat[j, ]) / dj else sum(mat[, i] & mat[, j]) / dj; s <- s + ov * 100 }
      np <- np + 1
    }
    c(s = s, np = np)
  }
  r <- pc(rowSums(M), M, "r"); cc <- pc(colSums(M), M, "c")
  as.numeric((r["s"] + cc["s"]) / (r["np"] + cc["np"]))
}
null_swap <- function(M, nsw = 2000) {
  for (s in seq_len(nsw)) {
    ii <- sample(nrow(M), 2); jj <- sample(ncol(M), 2); sub <- M[ii, jj]
    if (all(sub == matrix(c(1, 0, 0, 1), 2)))      M[ii, jj] <- matrix(c(0, 1, 1, 0), 2)
    else if (all(sub == matrix(c(0, 1, 1, 0), 2))) M[ii, jj] <- matrix(c(1, 0, 0, 1), 2)
  }
  M
}

A web whose truth you know

Every interaction matrix in the wild is a sample: the observer watched for a fixed number of hours and wrote down what happened. To see what that sampling does, define the true interaction probabilities first, then draw from them.

set.seed(4263)
A <- 20L; P <- 16L
a_row <- exp(rnorm(A, 0, 0.8)); b_col <- exp(rnorm(P, 0, 0.7))
mu <- runif(A); pos <- seq(0.05, 0.95, length.out = P); niche_w <- 0.22
Ptrue <- outer(a_row, b_col) * exp(-(outer(mu, pos, "-")^2) / (2 * niche_w^2))
Ptrue <- Ptrue / sum(Ptrue)
p_min <- 1e-6
n_possible <- sum(Ptrue > p_min)
c(cells = A * P, ecologically_possible = n_possible)
                cells ecologically_possible 
                  320                   315 

Of the 320 cells, 315 carry a non-negligible probability: those are the links that exist in this system. The rest are trait mismatches that will never be observed no matter how long you watch. Sampling n interactions is a multinomial draw over the cells.

The first diagnostic: how much of the web have you seen?

Before any index, ask how complete the sample is. The interaction accumulation curve counts distinct links discovered as interactions are recorded, the network analogue of a species accumulation curve. A curve still climbing at the end of the data is a web with links left to find.

set.seed(6)
grid_n <- c(seq(20, 400, 20), seq(500, 3000, 250))
pool <- sample(A * P, max(grid_n), replace = TRUE, prob = as.vector(Ptrue))
accum <- vapply(grid_n, function(k) length(unique(pool[seq_len(k)])), numeric(1))
acc_df <- data.frame(sampled = grid_n, unique_links = accum)
acc_end <- accum[length(accum)]
c(links_found_at_3000 = acc_end, links_that_exist = n_possible,
  fraction = round(acc_end / n_possible, 2))
links_found_at_3000    links_that_exist            fraction 
             237.00              315.00                0.75 

After 3000 recorded interactions the curve has found 237 of the 315 links that exist, 75 per cent, and is still rising. Real pollination datasets are usually far short of that (Chacoff et al 2012; Jordano 2016), which matters because the missing links are not a random subset: they are the rare ones.

Every metric drifts with effort

Now sample the same true web at increasing intensity and recompute the indices from the earlier posts.

samp_metrics <- function(n) {
  W <- matrix(rmultinom(1, n, as.vector(Ptrue)), A, P)
  Wt <- W[rowSums(W) > 0, colSums(W) > 0, drop = FALSE]
  c(effort = n, connectance = sum(W > 0) / (A * P),
    H2prime = H2prime(Wt), NODF = nodf((Wt > 0) * 1))
}
set.seed(5)
efforts <- c(50, 100, 200, 400, 800, 1600, 3200)
drift <- as.data.frame(t(vapply(efforts, samp_metrics, numeric(4))))
conn_lo <- drift$connectance[1]; conn_hi <- drift$connectance[nrow(drift)]
h2_lo <- drift$H2prime[1]; h2_hi <- drift$H2prime[nrow(drift)]
nodf_lo <- drift$NODF[1]; nodf_hi <- drift$NODF[nrow(drift)]
round(drift, 3)
  effort connectance H2prime   NODF
1     50       0.122   0.607 17.820
2    100       0.219   0.404 38.474
3    200       0.328   0.336 42.527
4    400       0.456   0.249 54.896
5    800       0.594   0.185 71.032
6   1600       0.681   0.147 76.533
7   3200       0.766   0.142 81.365

Nothing about the ecology changed between these rows. Only the watching did, and every number moved a long way. Connectance climbed from 0.12 to 0.77, because more effort finds more rare links. Nestedness rose from 18 to 81. Most striking, the specialisation index H2’ fell from 0.61 to 0.14: at low effort the web looks highly specialised, and the apparent specialisation dissolves as sampling improves. H2’ was introduced in the first post as the index that travels best because it standardises against the marginal totals. It still drifts by a factor of four here. A thin sample makes generalists look like specialists, which is the same confound that made poorly sampled pollinators score high on d’ (Bluthgen et al 2008; Frund et al 2016).

An accumulation curve rising steeply then bending but not levelling, and three metric trajectories against effort, two rising and the specialisation index falling.
Figure 1: Left: the interaction accumulation curve, still climbing after 3000 recorded interactions. Right: connectance, H2’ and nestedness against sampling effort on a log scale; all three move a long way although the underlying web never changes.

Size dependence: bigger webs look sparser

The second confound is dimension. Generate webs from the same process at different species richness, with effort scaled to richness so that per-species sampling is constant.

gen_web <- function(S, eff_per_sp = 32) {
  a <- exp(rnorm(S, 0, 0.7)); b <- exp(rnorm(S, 0, 0.6))
  mu_s <- runif(S); pos_s <- seq(0.05, 0.95, length.out = S)
  Pt <- outer(a, b) * exp(-(outer(mu_s, pos_s, "-")^2) / (2 * 0.22^2)); Pt <- Pt / sum(Pt)
  W <- matrix(rmultinom(1, eff_per_sp * S, as.vector(Pt)), S, S)
  W[rowSums(W) > 0, colSums(W) > 0, drop = FALSE]
}
set.seed(77)
sizes <- c(6, 9, 12, 15, 18, 21, 24)
size_res <- as.data.frame(t(vapply(sizes, function(S) {
  reps <- replicate(8, { W <- gen_web(S); c(conn = sum(W > 0) / (nrow(W) * ncol(W)), nd = nodf((W > 0) * 1)) })
  c(size = S, connectance = mean(reps["conn", ]), NODF = mean(reps["nd", ]))
}, numeric(3))))
size_conn_lo <- size_res$connectance[nrow(size_res)]; size_conn_hi <- size_res$connectance[1]
round(size_res, 3)
  size connectance   NODF
1    6       0.733 67.181
2    9       0.657 67.755
3   12       0.582 65.539
4   15       0.577 64.056
5   18       0.554 63.599
6   21       0.516 62.539
7   24       0.473 59.972

The niche width, the abundance spread and the effort per species are identical across these webs. Only the number of species differs, and connectance falls from 0.73 in the smallest web to 0.47 in the largest. This looks like the well-known fall of connectance with richness (Olesen and Jordano 2002; Bluthgen 2010), but here it is a sampling effect: the extra species fill the same trait axis more densely, so each species’ true partners grow with the web and nearly every cell is a possible link at every size, while the same number of interactions per species is spread over more cells. Holding the effort per cell constant instead leaves connectance flat (a check not shown here). Comparing raw connectance between a 20-species web and a 60-species web says more about their dimensions than their ecology.

What the null model rescues

Given that the raw numbers move, does the null-model machinery from the previous post fare better? Recompute nestedness across effort, but read each value against a fixed-degree null rather than on its own.

set.seed(41)
zt <- as.data.frame(t(vapply(c(100, 200, 400, 800, 1600, 3200), function(n) {
  W <- matrix(rmultinom(1, n, as.vector(Ptrue)), A, P)
  B <- (W > 0) * 1; B <- B[rowSums(B) > 0, colSums(B) > 0, drop = FALSE]
  obs <- nodf(B); nd <- replicate(99, nodf(null_swap(B)))
  c(effort = n, raw_NODF = obs, z = (obs - mean(nd)) / sd(nd))
}, numeric(3))))
z_raw_lo <- zt$raw_NODF[1]; z_raw_hi <- zt$raw_NODF[nrow(zt)]
z_range <- range(zt$z)
round(zt, 2)
  effort raw_NODF     z
1    100    32.35 -0.27
2    200    47.59 -0.64
3    400    61.14 -1.74
4    800    70.80  0.16
5   1600    80.09  0.46
6   3200    78.74 -3.09

The raw nestedness climbs from 32 to 79 across this range of effort, a factor of 2.4. The z-score does not follow it: in this run it wanders between -3.1 and 0.5 with no clear trend in effort. Five of the six efforts land within two standard deviations of their null. The sixth, at the highest effort, reaches -3.1, pointing at less nestedness than the degrees already imply, and that is probably not a chance excursion. The web was built by trait matching, which is structure beyond the degree sequence, and repeating the sampling many times at each effort (a check not shown here) moves the z-score downwards as effort grows, beyond minus two in most samples at the highest effort. What the null buys you is that the verdict no longer tracks the raw index: the raw index more than doubles across this range from fill-in alone, while the z-score moves only as far as a thicker sample gives the test power to see real structure. So compare z-scores only between webs sampled with similar effort. This is the practical case for reporting departures from a null rather than raw values, and it holds for modularity and specialisation as much as nestedness.

The z-score is noisy, though, and reading a single one as a precise effect size would overstate it. It answers a yes or no question about structure beyond the degrees; it does not restore the raw metric’s comparability.

A downward curve of connectance against richness, and a two-line panel where the raw nestedness climbs while the z-score does not follow it, wandering around zero at the lower efforts and dipping well below zero at the highest effort.
Figure 2: Left: connectance against species richness, with per-species effort held constant. Right: raw nestedness rises steeply with effort while its z-score against a fixed-degree null does not follow that drift in this run.

The honest limit: a zero has two meanings

Here is what no null model, no standardised index and no amount of care with the analysis can fix. Take a realistic sample and look at the empty cells.

set.seed(42)
W_obs <- matrix(rmultinom(1, 600, as.vector(Ptrue)), A, P)
zero_cells <- which(W_obs == 0)
thr <- p_min   # the same cutoff that defined the links that exist
n_forbidden <- sum(Ptrue[zero_cells] < thr)
n_missing   <- sum(Ptrue[zero_cells] >= thr)
exp_visits  <- 600 * Ptrue[zero_cells][Ptrue[zero_cells] >= thr]
c(zeros = length(zero_cells), forbidden = n_forbidden, missing = n_missing,
  pct_missing = round(100 * n_missing / length(zero_cells)),
  median_expected_visits = round(median(exp_visits), 2))
                 zeros              forbidden                missing 
                149.00                   5.00                 144.00 
           pct_missing median_expected_visits 
                 97.00                   0.17 

The sample has 149 empty cells. Of those, 5 are genuine: trait mismatches, phenological non-overlap, links that cannot happen. The other 144, which is 97 per cent of the zeros, are real interactions that simply were not seen, most of them expected to occur less than once at this effort (median 0.17 expected visits). In the matrix these two kinds of zero are the same character. Nothing in the data distinguishes a forbidden link from a missing one, and yet every index in this series treats them identically: connectance counts them as absences, NODF reads them as gaps in a nested pattern, H2’ scores them as evidence of choosiness.

Only outside information separates them. Trait measurements can rule a link impossible, phenological records can show two species never overlap, and more effort turns a missing link into an observed one (Olesen et al 2011; Jordano 2016). None of it comes from the matrix. This is the same shape as the honest limits elsewhere on this blog: the unseen species behind a richness estimate, the untestable exchangeability behind a conformal interval, the assumption of no unmeasured confounding behind a causal claim. The checking constrains the reading; it does not license certainty.

Where this leaves you

Three checks are worth running on any interaction network before the indices are reported. Plot the interaction accumulation curve and say plainly how complete the sample is. Refuse to compare raw connectance, linkage density or nestedness across webs of different size or effort, because those numbers scale with both. Report departures from an explicit null instead of raw values, since the z-score does not follow a raw index that more than doubles with effort, a drift that would otherwise be read as ecology; but compare z-scores only between webs of similar effort, because a thicker sample gives the test more power.

None of that reaches the deepest problem. Most of the zeros in a typical web are links that exist and were missed, and the matrix cannot tell you which. Network analysis is a description of a sample, and the further your conclusions travel from that sample, the more they lean on evidence that is not in it.

References

Bluthgen N 2010 Basic and Applied Ecology 11(3):185-195 (10.1016/j.baae.2010.01.001)

Bluthgen N, Frund J, Vazquez DP, Menzel F 2008 Ecology 89(12):3387-3399 (10.1890/07-2121.1)

Chacoff NP, Vazquez DP, Lomascolo SB, Stevani EL, Dorado J, Padron B 2012 Journal of Animal Ecology 81(1):190-200 (10.1111/j.1365-2656.2011.01883.x)

Frund J, McCann KS, Williams NM 2016 Oikos 125(4):502-513 (10.1111/oik.02256)

Jordano P 2016 Functional Ecology 30(12):1883-1893 (10.1111/1365-2435.12763)

Olesen JM, Jordano P 2002 Ecology 83(9):2416-2424 (10.1890/0012-9658(2002)083[2416:GPIPPM]2.0.CO;2)

Olesen JM, Bascompte J, Dupont YL, Elberling H, Rasmussen C, Jordano P 2011 Proceedings of the Royal Society B: Biological Sciences 278(1706):725-732 (10.1098/rspb.2010.1371)

Newsletter

Get updates by email

An occasional email when tutorials are added or substantially corrected. No spam; unsubscribe anytime.

By subscribing you agree to receive these emails and confirm your address once. See the privacy policy.