Checking a community phylogenetics analysis

community phylogenetics

model checking

ggplot2

ecology tutorial

Three checks that decide whether an NRI result means anything: sweep the species pool, calibrate your own pipeline on null data, and break the branch lengths.

Author

Tidy Ecology

Published

2026-07-15

The previous three posts built the machinery: NRI and NTI, the null models, the beta diversity metrics. All of it runs. None of it tells you whether your result survives the decisions you made before the analysis started.

Three checks do. One of them will probably kill your result, and it is better that it happens here.

Setup

A data set with a real, strong signal: forty sites, all the same habitat, filtering hard on a conserved trait. If any analysis deserves to report clustering, this one does.

library(ape)
library(ggplot2)

set.seed(292)
S  <- 200
tr <- rcoal(S, tip.label = sprintf("sp%03d", 1:S))
tr$edge.length <- tr$edge.length / max(node.depth.edgelength(tr))
pool  <- tr$tip.label
D     <- cophenetic(tr)[pool, pool]
trait <- rTraitCont(tr, model = "BM", sigma = 1)[pool]
trait <- (trait - mean(trait)) / sd(trait)

mpd <- function(sp, Dm = D) { d <- Dm[sp, sp]; mean(d[lower.tri(d)]) }

nri <- function(sp, pl, Dm = D, R = 199) {
  o   <- mpd(sp, Dm)
  nul <- replicate(R, mpd(sample(pl, length(sp)), Dm))
  -(o - mean(nul)) / sd(nul)
}

set.seed(43);  optC <- rnorm(40, 1.2, 0.25)
M <- matrix(0L, 40, S, dimnames = list(NULL, pool))
set.seed(431)
for (i in 1:40) M[i, sample(S, 15, prob = exp(-(trait - optC[i])^2 / (2 * 0.5^2)))] <- 1L
comms <- lapply(1:40, function(i) pool[M[i, ] == 1])

Check 1: sweep the pool

NRI compares your community against a draw from a species pool. The pool is not in your data. You chose it, usually by default, usually without saying so in the methods. Four choices, all defensible, all used in the literature:

obs_sp    <- pool[colSums(M) > 0]                                    # seen somewhere in the matrix
hab       <- pool[trait >= min(trait[colSums(M) > 0]) &
                  trait <= max(trait[colSums(M) > 0])]               # inside the observed trait range
plausible <- pool[abs(trait - 1.2) <= 1.2]                           # could tolerate this habitat

pools <- list(`whole tree`              = pool,
              `observed anywhere`       = obs_sp,
              `trait range of the data` = hab,
              `habitat-plausible`       = plausible)

sweep <- do.call(rbind, lapply(names(pools), function(nm) {
  set.seed(11)
  v <- sapply(comms, nri, pl = pools[[nm]])
  data.frame(pool = nm, size = length(pools[[nm]]),
             mean_NRI = round(mean(v), 2), pct_sig = round(100 * mean(abs(v) > 1.96)))
}))
sweep

                     pool size mean_NRI pct_sig
1              whole tree  200     3.62      78
2       observed anywhere  122     0.56      10
3 trait range of the data  136     1.03      22
4       habitat-plausible  124     0.71      12

Boxplots of NRI for four pool definitions. The whole-tree pool sits far above the significance threshold; the three narrower pools sit near zero. — Figure 1: NRI across the same 40 sites under four species pool definitions. The data never changed; only the answer to ‘random with respect to what’ did.

The whole-tree pool gives mean NRI +3.62 and flags 78 per cent of sites. Restricting the pool to species seen somewhere in the matrix, a narrowing most reviewers would call conservative and sensible, gives +0.56 and 10 per cent. The other two narrowings land in the same place.

Note what did not happen: the data did not change, the tree did not change, the filtering process is exactly as real under every row. What changed is the question. Against the whole tree, these communities are strongly clustered. Against the species that actually occur in this landscape, they are ordinary. Both statements are true, and only one of them is usually reported.

The sweep is three lines. Run it, put the table in your supplement, and if the sign or the significance moves, that is your result: the answer depends on a pool you cannot observe.

Check 2: run your own pipeline on nothing

Before blaming the statistic for check 1, find out whether the statistic works. Feed the pipeline data with no signal at all and count how often it cries wolf. If the standardised effect size behaves, five per cent of null communities should cross the threshold.

set.seed(12); fp15 <- replicate(400, abs(nri(sample(pool, 15), pool)) > 1.96)
set.seed(13); fp5  <- replicate(400, abs(nri(sample(pool,  5), pool)) > 1.96)
set.seed(14); z15  <- replicate(400, nri(sample(pool, 15), pool))

c(false_positive_k15 = mean(fp15), false_positive_k5 = mean(fp5),
  null_mean = mean(z15), null_sd = sd(z15))

false_positive_k15  false_positive_k5          null_mean            null_sd 
        0.05250000         0.06250000         0.01538244         0.97700793

Histogram of null NRI values closely overlaid by a standard normal density curve, centred on zero with unit width. — Figure 2: NRI for 400 communities of 15 species drawn at random from the pool, against the standard normal the 1.96 threshold assumes.

5.2 per cent at fifteen species, 6.2 per cent at five, against a nominal five. The null distribution has mean +0.02 and standard deviation 0.98. The arithmetic is fine. Whatever went wrong in check 1 was not the standardised effect size failing; it was the pool answering a different question, exactly as instructed.

This check costs almost nothing and it is worth doing on your own data rather than trusting this post, because your pipeline is not this pipeline. Anything above about seven per cent means a bug, usually a null draw that does not match the observed richness.

Check 3: break the branch lengths

Most community phylogenies are not estimated from your species. They are a published backbone with your taxa grafted on, polytomies resolved by hand, and branch lengths supplied by a rule. Ask what that costs by running the analysis on a tree that keeps the topology and throws the branch lengths away.

tr2 <- tr
tr2$edge.length <- rep(1, nrow(tr2$edge))
D2  <- cophenetic(tr2)[pool, pool]

set.seed(15); v_true <- sapply(comms, nri, pl = pool)
set.seed(15); v_eq   <- sapply(comms, nri, pl = pool, Dm = D2)

c(true_mean = mean(v_true), equal_mean = mean(v_eq),
  true_pct = mean(abs(v_true) > 1.96), equal_pct = mean(abs(v_eq) > 1.96),
  correlation = cor(v_true, v_eq))

  true_mean  equal_mean    true_pct   equal_pct correlation 
  3.5744651   1.3774118   0.7750000   0.2750000   0.4814676

Scatter plot of NRI under equal branch lengths against NRI under the real tree. Points sit well below the one-to-one line and scatter widely around it. — Figure 3: Per-site NRI from the real tree against NRI from the same topology with every branch length set to one. Points on the dashed line would mean branch lengths did not matter.

Mean NRI drops from +3.57 to +1.38, and the share of sites flagged falls from 78 to 28 per cent. Site by site, the two versions correlate at 0.48, so knowing a site’s NRI from the topology-only tree tells you roughly a quarter of what you want to know about its NRI from the real one.

Equal branch lengths are an extreme, but they are the honest bound on a tree whose lengths came from a rule rather than from data. If your tree is dated at the family level and interpolated below, your effective resolution is somewhere between these two columns, and nothing in the output will tell you where.

What to report

The three checks fit in one supplementary table and four sentences of methods:

Which pool, why that pool, and what the sweep across alternatives does to the sign and the significance.
The false positive rate of your own pipeline on communities you generated with no signal.
The correlation between your result and the same analysis on a degraded tree.

None of this makes a weak result strong. It makes a strong result believable, and it tells you early when a result is only an artefact of a pool you picked without noticing.

The honest limit

These checks bound the damage from choices you made. They cannot bound the damage from the thing community phylogenetics is actually built on: the assumption that phylogenetic distance stands in for ecological similarity. When that assumption fails, and it fails whenever the traits that matter are convergent or labile, every number in every one of these four posts is precise and about nothing. No amount of checking recovers from that, because the tree cannot tell you it happened. Measuring the traits can.

References

Cavender-Bares, Kozak, Fine, Kembel 2009 Ecology Letters 12(7):693-715 (10.1111/j.1461-0248.2009.01314.x)
Swenson 2009 PLoS ONE 4(2):e4390 (10.1371/journal.pone.0004390)
Kembel, Cowan, Helmus, Cornwell, Morlon, Ackerly, Blomberg, Webb 2010 Bioinformatics 26(11):1463-1464 (10.1093/bioinformatics/btq166)
Lessard, Borregaard, Fordyce, Rahbek, Weiser, Dunn, Sanders 2012 Proceedings of the Royal Society B 279(1727):266-274 (10.1098/rspb.2011.0552)

--- title: "Checking a community phylogenetics analysis" description: "Three checks that decide whether an NRI result means anything: sweep the species pool, calibrate your own pipeline on null data, and break the branch lengths." date: "2026-07-15 09:00" categories: [R, community phylogenetics, model checking, ggplot2, ecology tutorial] image: thumbnail.png image-alt: "Phylogenetic clustering scores across four species pool definitions, strong under the widest pool and absent under every narrower one." --- The previous three posts built the machinery: [NRI and NTI](../phylogenetic-community-structure-nri-nti/), [the null models](../null-models-for-phylogenetic-structure/), [the beta diversity metrics](../phylogenetic-beta-diversity/). All of it runs. None of it tells you whether your result survives the decisions you made before the analysis started. Three checks do. One of them will probably kill your result, and it is better that it happens here. ## Setup A data set with a real, strong signal: forty sites, all the same habitat, filtering hard on a conserved trait. If any analysis deserves to report clustering, this one does. ```{r setup} #| message: false #| warning: false library(ape) library(ggplot2) set.seed(292) S <- 200 tr <- rcoal(S, tip.label = sprintf("sp%03d", 1:S)) tr$edge.length <- tr$edge.length / max(node.depth.edgelength(tr)) pool <- tr$tip.label D <- cophenetic(tr)[pool, pool] trait <- rTraitCont(tr, model = "BM", sigma = 1)[pool] trait <- (trait - mean(trait)) / sd(trait) mpd <- function(sp, Dm = D) { d <- Dm[sp, sp]; mean(d[lower.tri(d)]) } nri <- function(sp, pl, Dm = D, R = 199) { o <- mpd(sp, Dm) nul <- replicate(R, mpd(sample(pl, length(sp)), Dm)) -(o - mean(nul)) / sd(nul) } set.seed(43); optC <- rnorm(40, 1.2, 0.25) M <- matrix(0L, 40, S, dimnames = list(NULL, pool)) set.seed(431) for (i in 1:40) M[i, sample(S, 15, prob = exp(-(trait - optC[i])^2 / (2 * 0.5^2)))] <- 1L comms <- lapply(1:40, function(i) pool[M[i, ] == 1]) ``` ## Check 1: sweep the pool NRI compares your community against a draw from a species pool. The pool is not in your data. You chose it, usually by default, usually without saying so in the methods. Four choices, all defensible, all used in the literature: ```{r pools} obs_sp <- pool[colSums(M) > 0] # seen somewhere in the matrix hab <- pool[trait >= min(trait[colSums(M) > 0]) & trait <= max(trait[colSums(M) > 0])] # inside the observed trait range plausible <- pool[abs(trait - 1.2) <= 1.2] # could tolerate this habitat pools <- list(`whole tree` = pool, `observed anywhere` = obs_sp, `trait range of the data` = hab, `habitat-plausible` = plausible) sweep <- do.call(rbind, lapply(names(pools), function(nm) { set.seed(11) v <- sapply(comms, nri, pl = pools[[nm]]) data.frame(pool = nm, size = length(pools[[nm]]), mean_NRI = round(mean(v), 2), pct_sig = round(100 * mean(abs(v) > 1.96))) })) sweep ``` ```{r fig-pool} #| echo: false #| fig-cap: "NRI across the same 40 sites under four species pool definitions. The data never changed; only the answer to 'random with respect to what' did." #| fig-alt: "Boxplots of NRI for four pool definitions. The whole-tree pool sits far above the significance threshold; the three narrower pools sit near zero." #| fig-width: 7 #| fig-height: 3.8 dfp <- do.call(rbind, lapply(names(pools), function(nm) { set.seed(11) data.frame(pool = nm, nri = sapply(comms, nri, pl = pools[[nm]])) })) dfp$pool <- factor(dfp$pool, levels = names(pools)) ggplot(dfp, aes(pool, nri, fill = pool)) + geom_hline(yintercept = 0, colour = "#16241d", linewidth = 0.3) + geom_hline(yintercept = 1.96, linetype = 2, colour = "#b5534e", linewidth = 0.4) + geom_boxplot(width = 0.6, colour = "#16241d", linewidth = 0.35, outlier.size = 0.7, outlier.colour = "#5d6b61") + scale_fill_manual(values = c("#275139", "#93a87f", "#93a87f", "#93a87f"), guide = "none") + labs(x = NULL, y = "NRI", subtitle = "dashed line: the usual significance threshold") + theme_minimal(base_size = 11) + theme(panel.grid.minor = element_blank(), axis.text.x = element_text(size = 9), plot.subtitle = element_text(colour = "#5d6b61", size = 9)) ``` The whole-tree pool gives mean NRI `r sprintf("%+.2f", sweep$mean_NRI[1])` and flags `r sweep$pct_sig[1]` per cent of sites. Restricting the pool to species seen somewhere in the matrix, a narrowing most reviewers would call conservative and sensible, gives `r sprintf("%+.2f", sweep$mean_NRI[2])` and `r sweep$pct_sig[2]` per cent. The other two narrowings land in the same place. Note what did not happen: the data did not change, the tree did not change, the filtering process is exactly as real under every row. What changed is the question. Against the whole tree, these communities are strongly clustered. Against the species that actually occur in this landscape, they are ordinary. Both statements are true, and only one of them is usually reported. The sweep is three lines. Run it, put the table in your supplement, and if the sign or the significance moves, that is your result: the answer depends on a pool you cannot observe. ## Check 2: run your own pipeline on nothing Before blaming the statistic for check 1, find out whether the statistic works. Feed the pipeline data with no signal at all and count how often it cries wolf. If the standardised effect size behaves, five per cent of null communities should cross the threshold. ```{r calib} set.seed(12); fp15 <- replicate(400, abs(nri(sample(pool, 15), pool)) > 1.96) set.seed(13); fp5 <- replicate(400, abs(nri(sample(pool, 5), pool)) > 1.96) set.seed(14); z15 <- replicate(400, nri(sample(pool, 15), pool)) c(false_positive_k15 = mean(fp15), false_positive_k5 = mean(fp5), null_mean = mean(z15), null_sd = sd(z15)) ``` ```{r fig-calib} #| echo: false #| fig-cap: "NRI for 400 communities of 15 species drawn at random from the pool, against the standard normal the 1.96 threshold assumes." #| fig-alt: "Histogram of null NRI values closely overlaid by a standard normal density curve, centred on zero with unit width." #| fig-width: 5.6 #| fig-height: 3.4 ggplot(data.frame(z = z15), aes(z)) + geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "#93a87f", colour = "#f5f4ee", linewidth = 0.2) + stat_function(fun = dnorm, colour = "#b5534e", linewidth = 1) + geom_vline(xintercept = c(-1.96, 1.96), linetype = 2, colour = "#5d6b61", linewidth = 0.4) + labs(x = "NRI under no signal", y = "density", subtitle = "red: standard normal") + theme_minimal(base_size = 11) + theme(panel.grid.minor = element_blank(), plot.subtitle = element_text(colour = "#5d6b61", size = 9)) ``` `r sprintf("%.1f", 100 * mean(fp15))` per cent at fifteen species, `r sprintf("%.1f", 100 * mean(fp5))` per cent at five, against a nominal five. The null distribution has mean `r sprintf("%+.2f", mean(z15))` and standard deviation `r sprintf("%.2f", sd(z15))`. The arithmetic is fine. Whatever went wrong in check 1 was not the standardised effect size failing; it was the pool answering a different question, exactly as instructed. This check costs almost nothing and it is worth doing on your own data rather than trusting this post, because your pipeline is not this pipeline. Anything above about seven per cent means a bug, usually a null draw that does not match the observed richness. ## Check 3: break the branch lengths Most community phylogenies are not estimated from your species. They are a published backbone with your taxa grafted on, polytomies resolved by hand, and branch lengths supplied by a rule. Ask what that costs by running the analysis on a tree that keeps the topology and throws the branch lengths away. ```{r tree-break} tr2 <- tr tr2$edge.length <- rep(1, nrow(tr2$edge)) D2 <- cophenetic(tr2)[pool, pool] set.seed(15); v_true <- sapply(comms, nri, pl = pool) set.seed(15); v_eq <- sapply(comms, nri, pl = pool, Dm = D2) c(true_mean = mean(v_true), equal_mean = mean(v_eq), true_pct = mean(abs(v_true) > 1.96), equal_pct = mean(abs(v_eq) > 1.96), correlation = cor(v_true, v_eq)) ``` ```{r fig-tree} #| echo: false #| fig-cap: "Per-site NRI from the real tree against NRI from the same topology with every branch length set to one. Points on the dashed line would mean branch lengths did not matter." #| fig-alt: "Scatter plot of NRI under equal branch lengths against NRI under the real tree. Points sit well below the one-to-one line and scatter widely around it." #| fig-width: 5 #| fig-height: 4 ggplot(data.frame(t = v_true, e = v_eq), aes(t, e)) + geom_abline(linetype = 2, colour = "#5d6b61", linewidth = 0.4) + geom_hline(yintercept = 1.96, linetype = 3, colour = "#b5534e", linewidth = 0.4) + geom_point(colour = "#275139", size = 2, alpha = 0.75) + labs(x = "NRI, real branch lengths", y = "NRI, all branches equal") + theme_minimal(base_size = 11) + theme(panel.grid.minor = element_blank()) ``` Mean NRI drops from `r sprintf("%+.2f", mean(v_true))` to `r sprintf("%+.2f", mean(v_eq))`, and the share of sites flagged falls from `r sprintf("%.0f", 100 * mean(abs(v_true) > 1.96))` to `r sprintf("%.0f", 100 * mean(abs(v_eq) > 1.96))` per cent. Site by site, the two versions correlate at `r sprintf("%.2f", cor(v_true, v_eq))`, so knowing a site's NRI from the topology-only tree tells you roughly a quarter of what you want to know about its NRI from the real one. Equal branch lengths are an extreme, but they are the honest bound on a tree whose lengths came from a rule rather than from data. If your tree is dated at the family level and interpolated below, your effective resolution is somewhere between these two columns, and nothing in the output will tell you where. ## What to report The three checks fit in one supplementary table and four sentences of methods: - Which pool, why that pool, and what the sweep across alternatives does to the sign and the significance. - The false positive rate of your own pipeline on communities you generated with no signal. - The correlation between your result and the same analysis on a degraded tree. None of this makes a weak result strong. It makes a strong result believable, and it tells you early when a result is only an artefact of a pool you picked without noticing. ## The honest limit These checks bound the damage from choices you made. They cannot bound the damage from the thing community phylogenetics is actually built on: the assumption that phylogenetic distance stands in for ecological similarity. When that assumption fails, and it fails whenever the traits that matter are convergent or labile, every number in every one of these four posts is precise and about nothing. No amount of checking recovers from that, because the tree cannot tell you it happened. Measuring the traits can. ## Related tutorials - [Phylogenetic community structure: NRI and NTI in R](../phylogenetic-community-structure-nri-nti/) - [Null models for phylogenetic community structure](../null-models-for-phylogenetic-structure/) - [Phylogenetic beta diversity: UniFrac and PhyloSor](../phylogenetic-beta-diversity/) - [Phylogenetic signal: Blomberg K and Pagel lambda](../phylogenetic-signal-k-lambda/) ## References - Cavender-Bares, Kozak, Fine, Kembel 2009 Ecology Letters 12(7):693-715 (10.1111/j.1461-0248.2009.01314.x) - Swenson 2009 PLoS ONE 4(2):e4390 (10.1371/journal.pone.0004390) - Kembel, Cowan, Helmus, Cornwell, Morlon, Ackerly, Blomberg, Webb 2010 Bioinformatics 26(11):1463-1464 (10.1093/bioinformatics/btq166) - Lessard, Borregaard, Fordyce, Rahbek, Weiser, Dunn, Sanders 2012 Proceedings of the Royal Society B 279(1727):266-274 (10.1098/rspb.2011.0552)