library(vegan)
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
prune <- function(M){
M <- M[, colSums(M) > 0 & colSums(M) < nrow(M), drop = FALSE]
M <- M[rowSums(M) > 0, , drop = FALSE]
M
}
# Segregated: two guilds with opposite responses to a site gradient.
make_segregated <- function(){
set.seed(11); nsite <- 30; nsp <- 24
grad <- seq(-1.5, 1.5, length.out = nsite)
guild <- rep(c(-1, 1), length.out = nsp)
M <- matrix(0L, nsite, nsp)
for (j in seq_len(nsp)) {
p <- plogis(-0.2 + 3.0 * guild[j] * grad + rnorm(nsite, 0, 0.3))
M[, j] <- rbinom(nsite, 1, p)
}
prune(M)
}
# Random: independent Bernoulli presences, no structure.
make_random <- function(){
set.seed(11); nsite <- 30; nsp <- 24
prune(matrix(rbinom(nsite * nsp, 1, 0.45), nsite, nsp))
}
# Aggregated: every species tracks one shared habitat axis.
make_aggregated <- function(){
set.seed(7); nsite <- 30; nsp <- 24
hab <- seq(-1.4, 1.4, length.out = nsite)
M <- matrix(0L, nsite, nsp)
for (j in seq_len(nsp)) {
a <- runif(1, -0.6, 0.2)
p <- plogis(a + 1.8 * hab + rnorm(nsite, 0, 0.35))
M[, j] <- rbinom(nsite, 1, p)
}
prune(M)
}
Seg <- make_segregated()
Ran <- make_random()
Agg <- make_aggregated()Species co-occurrence null models in R
R
community ecology
biodiversity
vegan
ecology tutorial
Test for non-random species co-occurrence in R with vegan: the C-score, oecosimu, the fixed-fixed (quasiswap) null model, and standardised effect sizes.
Corrected 4 November 2026: the post said Connor and Simberloff (1979) made the point that the same pattern can come from dispersal limitation, habitat filtering or historical accident; their abstract supports that Diamond’s assembly rules are circular, trivial or what random colonisation would also produce, and that showing competition needs a null hypothesis of chance to be rejected, and the text now says that and gives the list of other mechanisms as the post’s own.
Corrected 14 October 2026: the introduction gave the fixed-fixed model’s conservatism as the reason it does not flag the aggregated community; it now gives the reason the later section gives, that the aggregation built into that community lies in the margins the model holds fixed.
Corrected 1 October 2026: the post called the fixed-fixed model’s silence on the aggregated community a Type II error; it now says that aggregation is made of site and species effects alone, which the fixed margins already contain, so that null is essentially true of the community and the aggregation is invisible to it by design.
Corrected 30 September 2026: the post said Gotelli showed that freeing site richness inflates false positives, making the free-richness result partly an artefact, and that he documented the C-score missing aggregation under the fixed-fixed model; it now says his equivalent model was among those invulnerable to false positives, that freeing richness changes the question instead, and that his Type II warnings concern other indices. It also said he found the fixed-fixed model the safest against false positives; it now says every model that holds species totals fixed shared that safety.
Look at a site-by-species table and ask a simple question: do these species avoid each other, ignore each other, or travel together? The pattern that started the argument is the checkerboard, where two species never share a site, as if each excludes the other (Diamond 1975). The trouble is that a raw count of checkerboards means nothing on its own. A common species and a rare one will rarely co-occur whatever the ecology, and a table where every site holds the same number of species constrains co-occurrence before any process acts. To say a pattern is non-random you need a model of what random looks like, and that model is a choice with consequences (Connor & Simberloff 1979).
This post builds three synthetic communities, scores their co-occurrence with the C-score, and tests each against a null model with vegan::oecosimu. Two of the three behave as you would hope. The third does not, and the reason it does not is the most useful thing here: what the fixed-fixed null model can detect depends on which margins of the table it holds fixed, and the aggregation built into the third community lies in those margins.
Three communities
Each community is a sites-by-species incidence matrix (1 for present, 0 for absent). The segregated one has two guilds that respond to a gradient in opposite directions, so they tend to occupy different sites. The random one has independent presences. The aggregated one has every species respond to a single shared habitat axis, so they pile into the same good sites.
The C-score is the engine. For a pair of species it counts the checkerboard units, (Ri - S)(Rj - S), where Ri and Rj are the number of sites each occupies and S is the number of sites they share (Stone & Roberts 1990). Two species that never co-occur reach the maximum for their occurrences; two that always co-occur score zero. vegan::nestedchecker averages this over every species pair, so a high C-score means segregation and a low one means aggregation. The seriated incidence matrices below show the three structures: a clean two-block split, a speckle, and a filled gradient.
ink <- "#16241d"; forest_col <- "#275139"; gold <- "#cda23f"; red <- "#b5534e"
faint <- "#5d6b61"; line_col <- "#dad9ca"; paper <- "#f5f4ee"; absent <- "#e8e7dc"
te <- theme_minimal(base_size = 12) + theme(
panel.grid.minor = element_blank(),
panel.grid.major = element_line(color = line_col, linewidth = 0.3),
plot.background = element_rect(fill = paper, color = NA),
panel.background = element_rect(fill = paper, color = NA),
strip.text = element_text(face = "bold", color = ink, hjust = 0),
axis.title = element_text(color = "#46604a"),
legend.position = "top", legend.title = element_blank())
serorder <- function(M){
sv <- svd(scale(M, center = TRUE, scale = FALSE))
list(r = order(sv$u[, 1]), c = order(sv$v[, 1]))
}
tidymat <- function(M, lab){
o <- serorder(M); M2 <- M[o$r, o$c, drop = FALSE]
df <- expand.grid(site = seq_len(nrow(M2)), species = seq_len(ncol(M2)))
df$present <- as.vector(M2)
df$community <- lab; df
}
mats <- bind_rows(tidymat(Seg, "segregated"),
tidymat(Ran, "random"),
tidymat(Agg, "aggregated"))
mats$community <- factor(mats$community, levels = c("segregated", "random", "aggregated"))
ggplot(mats, aes(species, site, fill = factor(present))) +
geom_tile(color = paper, linewidth = 0.25) +
facet_wrap(~community) +
scale_fill_manual(values = c(`0` = absent, `1` = forest_col), guide = "none") +
scale_y_reverse() +
labs(x = "species (seriated)", y = "site (seriated)") + te +
theme(axis.text = element_blank(), panel.grid = element_blank())
Scoring against a null model
The raw C-score cannot be read directly, because its scale depends on how common the species are and how rich the sites are. oecosimu solves this by shuffling the matrix many times under a null model, recomputing the C-score on each shuffle, and reporting where the observed value falls in that null distribution. The default tool here is the fixed-fixed null, generated by the quasiswap algorithm: it keeps both the row sums (site richness) and the column sums (species occurrence) exactly as observed, and randomises only which cells are filled within those constraints. Gotelli (2000) compared nine such algorithms and found that those holding each species’ occurrence total fixed, the fixed-fixed model included, were all invulnerable to false positives, which is one reason it is the common default.
ff <- function(M, alt){
set.seed(42)
oecosimu(M, nestedchecker, method = "quasiswap", nsimul = 499, alternative = alt)
}
os_seg <- ff(Seg, "greater") # test for segregation: more checkerboards than null
os_ran <- ff(Ran, "two.sided") # random control
os_agg <- ff(Agg, "less") # test for aggregation: fewer checkerboards than null
ses <- function(os) as.numeric(os$oecosimu$z[1])
pval <- function(os) as.numeric(os$oecosimu$pval[1])
# effect size and p-value for each community under the fixed-fixed null
data.frame(community = c("segregated", "random", "aggregated"),
SES = round(c(ses(os_seg), ses(os_ran), ses(os_agg)), 1),
p = round(c(pval(os_seg), pval(os_ran), pval(os_agg)), 3)) community SES p
1 segregated 82.4 0.002
2 random 0.8 0.446
3 aggregated -0.8 0.210The standardised effect size (SES) is the observed C-score minus the null mean, divided by the null standard deviation, so it reads like a z-score: large positive means segregation, large negative means aggregation, near zero means the data sit inside the cloud of what randomness produces.
The segregated community lands far in the upper tail, with SES 82.4 and p = 0.002 (the smallest value 499 shuffles can return), so its checkerboard structure is well outside the null. The random community is the control, and it gives SES 0.8 and p = 0.446, which is what the test returns on a matrix built with no co-occurrence structure at all. The aggregated community, tested for fewer checkerboards than expected, gives SES -0.8 and p = 0.210. That is not significant. A community built so that every species shares the same sites is not flagged as aggregated.
meannull <- function(os, M) as.numeric(os$oecosimu$simulated[1, ]) / choose(ncol(M), 2)
obsmean <- function(M) nestedchecker(M)$C.score
nulldf <- bind_rows(
data.frame(C = meannull(os_seg, Seg), community = "segregated"),
data.frame(C = meannull(os_agg, Agg), community = "aggregated"))
nulldf$community <- factor(nulldf$community, levels = c("segregated", "aggregated"))
obsdf <- data.frame(
community = factor(c("segregated", "aggregated"), levels = c("segregated", "aggregated")),
obs = c(obsmean(Seg), obsmean(Agg)),
lab = c(sprintf("SES = %.1f\np = %.3f", ses(os_seg), pval(os_seg)),
sprintf("SES = %.1f\np = %.3f", ses(os_agg), pval(os_agg))),
lab_x = c(-Inf, Inf), lab_h = c(-0.15, 1.1)) # segregated: top left; aggregated: top right, clear of the red line
ggplot(nulldf, aes(C)) +
geom_histogram(bins = 28, fill = line_col, color = paper, linewidth = 0.2) +
geom_vline(data = obsdf, aes(xintercept = obs), color = red, linewidth = 0.9) +
geom_text(data = obsdf, aes(x = obs, y = Inf, label = "observed"), color = red,
angle = 90, hjust = 1.1, vjust = -0.4, size = 3.0) +
geom_text(data = obsdf, aes(x = lab_x, y = Inf, label = lab, hjust = lab_h), color = ink,
vjust = 1.5, size = 3.2, lineheight = 0.95) +
facet_wrap(~community, scales = "free") +
labs(x = "C-score under fixed-fixed null (mean checkerboard units per pair)",
y = "simulations") +
scale_x_continuous(expand = expansion(mult = 0.06)) +
scale_y_continuous(expand = expansion(mult = c(0.05, 0.3))) + # headroom so the labels clear the bars
te + theme(legend.position = "none")
Why aggregation escapes, and what the null model decides
The non-result is not a coding accident. Aggregation, by piling species into the same good sites, makes site richness uneven: the good sites are full, the poor sites are empty. The fixed-fixed null holds those row sums fixed, so it preserves the very richness pattern that aggregation produced. Held against a null that already contains its own footprint, the observed matrix looks unremarkable. That is not a Type II error, in which a false null goes unrejected: every species responds to the shared habitat axis with the same slope and to none of the other species, so the pattern is made of site and species effects alone, which the fixed row and column sums already contain; the fixed-fixed null is essentially true of this community, and real aggregation of this kind is invisible to it by design. Gotelli (2000) warned that two other indices, the number of checkerboard combinations and the number of species combinations, may miss patterns in noisy data, and found that the C-score used here had good power for detecting species pairs that rarely co-occur; the next figure shows that the blind spot in this example comes from the null model, not the index.
The fix and its cost become clear when you relax which margins are held. Below, the same aggregated community is tested against four null models, from the strict fixed-fixed to ones that free a margin.
methlab <- c(quasiswap = "fixed-fixed\n(quasiswap)",
r1 = "rows fixed,\ncols proportional (r1)",
r0 = "rows fixed,\ncols equal (r0)",
c0 = "cols fixed,\nrows free (c0)")
zz <- c(); pp <- c()
for (m in c("quasiswap", "r1", "r0", "c0")){
set.seed(42)
os <- oecosimu(Agg, nestedchecker, method = m, nsimul = 499, alternative = "less")
zz[m] <- as.numeric(os$oecosimu$z[1]); pp[m] <- as.numeric(os$oecosimu$pval[1])
}
fig3 <- data.frame(method = factor(methlab[names(zz)],
levels = methlab[c("quasiswap", "r1", "r0", "c0")]),
SES = as.numeric(zz), p = as.numeric(pp))
fig3$sig <- ifelse(fig3$p < 0.05, "p < 0.05", "n.s.")
ggplot(fig3, aes(SES, method)) +
geom_vline(xintercept = 0, color = faint, linewidth = 0.4) +
geom_segment(aes(x = 0, xend = SES, yend = method), color = line_col, linewidth = 1.1) +
geom_point(aes(color = sig), size = 4) +
geom_text(aes(label = sprintf("p = %.3f", p), hjust = ifelse(SES >= 0, -0.08, 1.08)),
nudge_y = 0.30, color = ink, size = 3.2) +
scale_color_manual(values = c(`p < 0.05` = forest_col, `n.s.` = faint)) +
scale_x_continuous(limits = c(-23, 4), breaks = seq(-20, 0, 5)) +
labs(x = "standardised effect size (SES) for aggregation", y = NULL) + te
The fixed-fixed model returns SES -0.8 (p = 0.210), no signal. The r1 model, which keeps site richness and makes species occurrence proportional to the observed totals, lands on the other side of zero at SES 0.3 (p = 0.558), no signal either. The r0 model, which fixes site richness but lets species occurrence vary freely, edges toward it at SES -1.4 (p = 0.074). The c0 model, which fixes species occurrence but lets site richness vary, declares strong aggregation at SES -19.2 (p = 0.002). Same matrix, opposite conclusions.
The temptation is to read c0 as the model that finally got it right. Resist it, though not because it is prone to false positives: in Gotelli’s (2000) tests the three algorithms that fix each species’ occurrence total were invulnerable to Type I errors, and his fixed-equiprobable model, which is what c0 does, is one of them. What c0 changes is the question. With site richness no longer held fixed, it asks whether species share sites more than they would if every site were equally able to hold them, and the uneven richness that the shared habitat axis produced counts as co-occurrence under that null. The c0 result says the species pile into the same sites; it does not say whether they respond to each other or to the habitat. The lesson is not “use c0 to find aggregation.” It is that the null model is the hypothesis. The fixed-fixed model trades power for safety; relaxing it buys power, at the risk of false positives when the species totals are freed and at the cost of a different question when site richness is freed; and the analyst, not the software, decides which trade is acceptable for the question.
How far to read these tests
A significant C-score quantifies segregation; it does not name a cause. Diamond’s original reading was competition, but Connor & Simberloff (1979) argued that his assembly rules are circular, trivial, or what random colonisation would produce anyway, and that showing competition means rejecting a null hypothesis of independent, chance colonisation. The same pattern can also come from dispersal limitation, shared habitat filtering, or historical accident, so pattern alone cannot separate them. Stone & Roberts (1992) showed something sharper still: strong checkerboardedness can be produced by forces driving aggregation just as well as by avoidance, so even the sign of the departure does not map cleanly onto a mechanism. Across many real matrices, segregation is the common finding (Gotelli & McCabe 2002), but reading it as evidence of competition requires more than a low p-value.
A few practical notes. Set a seed before each oecosimu call, since the p-value comes from a finite number of shuffles and is otherwise irreproducible. With 499 simulations the smallest p-value is about 0.002, so report effect sizes alongside p-values rather than leaning on significance alone. oecosimu tracks the null distribution for the total checkerboard count and reports the per-pair C-score as a derived value; the two differ only by a constant, so the SES is identical either way. And the C-score is one of several co-occurrence indices: the number of species combinations and the variance ratio answer related but distinct questions, and the index you pick, like the null model you pick, is part of the hypothesis you are testing.
References
Diamond JM 1975 In: Cody ML, Diamond JM (eds), Ecology and Evolution of Communities. Harvard University Press, pp 342-444
Connor EF, Simberloff D 1979 Ecology 60(6):1132-1140 (10.2307/1936961)
Stone L, Roberts A 1990 Oecologia 85(1):74-79 (10.1007/BF00317345)
Stone L, Roberts A 1992 Oecologia 91(3):419-424 (10.1007/BF00317632)
Gotelli NJ 2000 Ecology 81(9):2606-2621 (10.1890/0012-9658(2000)081[2606:NMAOSC]2.0.CO;2)
Gotelli NJ, McCabe DJ 2002 Ecology 83(8):2091-2096 (10.1890/0012-9658(2002)083[2091:SCOAMA]2.0.CO;2)
Oksanen J, Simpson G, Blanchet F, et al. 2022 vegan: Community Ecology Package, R package version 2.6-4