Phylogenetic community structure: NRI and NTI in R

community phylogenetics

biodiversity

ggplot2

ecology tutorial

Raw MPD and MNTD are confounded with species richness in a plot. Build standardised effect sizes, NRI and NTI, by hand in R and learn to read their signs.

Author

Tidy Ecology

Published

2026-06-14

Modified

2026-10-08

Corrected 8 October 2026: the post explained the wide spread of the filtered NRI by trait optima landing in crowded or empty parts of the pool, which the runs do not bear out, and said that the spacing rule rejects candidates close to a resident; it now puts the spread down to where each optimum falls and says that the rule down-weights close candidates.

Corrected 1 October 2026: the post called the index calibrated on the strength of its mean and spread, pointed to the checking post for a case it does not contain, and said that limiting similarity showing up mainly in NTI follows from the metrics alone; it now notes that the skewed MNTD null puts the dashed cut-offs off their nominal rates and that the spacing rule here acts directly on distances in the tree.

You have a species list for a plot and a phylogeny for the regional flora. Are the species in that plot closer relatives than you would expect, and if so, does that say anything about how the community was assembled?

The metric part of this question is easy. The interpretive part is where the whole field lives, and it starts with a fact that is easy to miss: the raw number you compute from a plot is not a property of the plot. For MNTD it is mostly a property of how many species the plot has, and even for MPD how unusual a given value is depends on the species count.

This post builds the standard fix from scratch: the standardised effect size, and its two sign-flipped aliases, the Net Relatedness Index (NRI) and the Nearest Taxon Index (NTI). Everything is base R plus ape for the tree.

Setup

A regional pool of 200 species on a coalescent tree, scaled so the root sits at depth 1 (the largest possible distance between two tips is therefore 2). One trait evolves on the tree under Brownian motion, so close relatives resemble each other: that is the assumption the whole method rests on, and it is worth naming out loud.

library(ape)
library(ggplot2)
te_canvas <- theme(plot.background  = element_rect(fill = "#f5f4ee", colour = NA),
                   panel.background = element_rect(fill = "#f5f4ee", colour = NA))

set.seed(292)
S <- 200
pool_tree <- rcoal(S, tip.label = sprintf("sp%03d", 1:S))
pool_tree$edge.length <- pool_tree$edge.length / max(node.depth.edgelength(pool_tree))
pool <- pool_tree$tip.label

D <- cophenetic(pool_tree)[pool, pool]
trait <- rTraitCont(pool_tree, model = "BM", sigma = 1)[pool]
trait <- (trait - mean(trait)) / sd(trait)

cophenetic() gives the patristic distance matrix: the summed branch length on the path between each pair of tips. Everything below reads from that matrix and nothing else.

Two distances, two depths

Two summaries dominate the literature, and they are two lines each.

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

mntd <- function(sp) {
  d <- D[sp, sp]
  diag(d) <- NA
  mean(apply(d, 1, min, na.rm = TRUE))
}

Mean pairwise distance (MPD) averages over every pair in the community. Mean nearest taxon distance (MNTD) averages, over each species, the distance to its single closest relative in the community.

They are not two flavours of the same thing. MPD is dominated by the deep splits: swap one species for a member of a distant clade and every pair involving it changes. MNTD only looks at each species’ nearest neighbour, so it is a measurement of the terminal structure. A community can be built from two distant clades (high MPD) with tight species clusters inside each (low MNTD). Keep that picture; it explains most of the confusing results later.

The raw number depends on richness

Here is the experiment nobody runs before reporting an MNTD. Draw communities at random from the pool at five richness levels, with no assembly rule whatsoever, and look at what the raw metrics do.

neutral <- function(k) sample(pool, k)

set.seed(1)
kg <- rep(c(4, 8, 15, 25, 40), each = 200)
raw <- t(sapply(kg, function(k) { s <- neutral(k); c(mpd(s), mntd(s)) }))
tab <- data.frame(
  k      = sort(unique(kg)),
  mpd_m  = tapply(raw[, 1], kg, mean), mpd_s  = tapply(raw[, 1], kg, sd),
  mntd_m = tapply(raw[, 2], kg, mean), mntd_s = tapply(raw[, 2], kg, sd)
)
round(tab, 3)

    k mpd_m mpd_s mntd_m mntd_s
4   4 1.129 0.287  0.634  0.243
8   8 1.129 0.183  0.342  0.135
15 15 1.144 0.091  0.204  0.074
25 25 1.142 0.064  0.129  0.044
40 40 1.150 0.046  0.084  0.026

Two panels against species richness. The MPD panel shows a flat cloud whose mean stays near 1.14 while the spread narrows sharply. The MNTD panel shows a cloud whose mean falls steeply from about 0.63 at richness 4 to about 0.08 at richness 40. — Figure 1: Raw MPD and MNTD under purely random draws from the pool, 200 replicates per richness level. Points are replicates, the line joins the means. No assembly process is present in any of these communities.

The MNTD mean falls by a factor of 7.57 between richness 4 and richness 40. Nothing ecological happened. Add species to a fixed pool and each species’ nearest relative simply gets nearer, because there are more candidates. Any comparison of raw MNTD between a species-poor plot and a species-rich plot is a comparison of richness wearing a phylogenetic costume.

MPD looks safer at first: the mean barely moves (1.13 to 1.15). But look at the second column. The standard deviation shrinks by a factor of 6.21 over the same range. A four-species plot sitting 0.2 below the pool mean is unremarkable; a forty-species plot sitting 0.2 below it is extraordinary. The raw deviation means different things at different richness, which is the same problem in a quieter voice.

Standardising against a null

The fix is the obvious one. Ask what the metric would look like if this community, with exactly this richness, had been drawn at random from the pool, and express the observation in units of that distribution.

ses <- function(sp, R = 999) {
  k   <- length(sp)
  obs <- c(mpd(sp), mntd(sp))
  nul <- replicate(R, { s <- sample(pool, k); c(mpd(s), mntd(s)) })
  z   <- (obs - rowMeans(nul)) / apply(nul, 1, sd)
  p   <- (rowSums(nul <= obs) + 1) / (R + 1)
  list(obs = obs, null = nul, z = z, p = p, nri = -z[1], nti = -z[2])
}

Three conventions are buried in those few lines, and each one has cost somebody a paragraph in review.

The sign is flipped. NRI is defined as minus the standardised MPD, and NTI as minus the standardised MNTD. So positive NRI means clustering: the community is more closely related than chance, which is a low MPD. If you report ses.mpd from picante and then discuss it as if it were NRI, every sign in your discussion is backwards. This is the single most common slip in the method.

The p-value is a rank, not a z-test. The p line counts how many null draws fall at or below the observation. That works whatever shape the null distribution has, which matters more than it sounds: the MNTD null in the figure below is skewed to the right, and for a skewed null a z-score cut-off and the rank p-value can disagree.

The null is a decision, not a default. Here the null draws species uniformly from the pool. That is one hypothesis out of several, and it is the subject of the next post.

Three assembly rules, three signatures

Now build communities under three rules and see what the index reports. Habitat filtering picks species whose trait is near a local optimum. Limiting similarity picks species that are far from whatever has already been picked. The neutral rule picks at random.

filtered <- function(k, opt, w = 0.3) {
  sample(pool, k, prob = exp(-(trait - opt)^2 / (2 * w^2)))
}

spaced <- function(k) {
  sel <- sample(pool, 1)
  while (length(sel) < k) {
    cand <- setdiff(pool, sel)
    sel  <- c(sel, sample(cand, 1, prob = apply(D[cand, sel, drop = FALSE], 1, min)^8))
  }
  sel
}

set.seed(10); cf <- filtered(15, 1.2)
set.seed(11); cn <- neutral(15)
set.seed(12); co <- spaced(15)

set.seed(20); sf <- ses(cf)
set.seed(21); sn <- ses(cn)
set.seed(22); so <- ses(co)

data.frame(
  rule = c("filtered", "neutral", "spaced"),
  MPD  = round(c(sf$obs[1], sn$obs[1], so$obs[1]), 3),
  NRI  = round(c(sf$nri, sn$nri, so$nri), 2),
  MNTD = round(c(sf$obs[2], sn$obs[2], so$obs[2]), 3),
  NTI  = round(c(sf$nti, sn$nti, so$nti), 2)
)

      rule   MPD   NRI  MNTD   NTI
1 filtered 0.724  4.67 0.029  2.29
2  neutral 1.123  0.20 0.140  0.85
3   spaced 1.294 -1.62 0.433 -3.25

Two histograms of null values in sage green. The MPD null is centred near 1.15 and the observed value 0.72 sits well outside its left tail. The MNTD null is centred near 0.20 with a right skew, and the observed value 0.03 sits at its left edge. — Figure 2: Null distributions for the filtered community (15 species, 999 draws) with the observed value marked. The observation sits far into the left tail of both distributions, which is what phylogenetic clustering looks like before it is turned into a z-score.

The filtered community, which really was assembled by filtering, lands at NRI +4.67 with a rank p of 0.001. The neutral community lands at NRI +0.20, p 0.424, against a truth of no structure at all. The spaced community lands at NTI -3.25: overdispersed at the tips.

Repeat over 100 communities per rule to see how reliable that is.

rep_ses <- function(gen, n = 100, R = 199) {
  t(replicate(n, {
    sp  <- gen(15)
    obs <- c(mpd(sp), mntd(sp))
    nul <- replicate(R, { s <- sample(pool, 15); c(mpd(s), mntd(s)) })
    -(obs - rowMeans(nul)) / apply(nul, 1, sd)
  }))
}

set.seed(30); A <- rep_ses(function(k) filtered(k, rnorm(1)))
set.seed(31); B <- rep_ses(neutral)
set.seed(32); C <- rep_ses(spaced)

data.frame(
  rule    = c("filtered", "neutral", "spaced"),
  NRI     = round(c(mean(A[, 1]), mean(B[, 1]), mean(C[, 1])), 2),
  NRI_sd  = round(c(sd(A[, 1]), sd(B[, 1]), sd(C[, 1])), 2),
  NTI     = round(c(mean(A[, 2]), mean(B[, 2]), mean(C[, 2])), 2),
  NTI_sd  = round(c(sd(A[, 2]), sd(B[, 2]), sd(C[, 2])), 2)
)

      rule   NRI NRI_sd   NTI NTI_sd
1 filtered  3.39   3.76  1.33   0.90
2  neutral  0.06   1.09  0.17   1.04
3   spaced -1.57   0.13 -3.09   0.18

Boxplots of NRI and NTI for three rules. Filtered communities sit well above zero on NRI with a long upper tail and modestly above zero on NTI. Neutral communities straddle zero on both. Spaced communities sit below zero on both, much further below on NTI. — Figure 3: NRI and NTI over 100 communities per assembly rule, 15 species each. A solid line marks zero and dashed lines mark the conventional two-tailed cut-offs. Positive is clustered, negative is overdispersed.

The neutral row is the one to read first. Its NRI averages +0.06 with a spread of 1.09, and its NTI averages +0.17 with a spread of 1.04. Centred on zero, spread of about one: that is what a standardised effect size is supposed to do when the data really were generated by the null. Centre and spread are not the whole story, though: the MNTD null is skewed to the right (the figure above), so under the null NTI almost never reaches the upper dashed line and passes the lower one a little more often than a normal curve would. The rank p-value, which assumes no shape, is the one to test with.

Why NRI and NTI disagree

Now the asymmetry. Habitat filtering shows up mainly in NRI (+3.39, against +1.33 for NTI). Limiting similarity shows up mainly in NTI (-3.09, against -1.57 for NRI).

Part of that is built into these simulations. A conserved trait varies most between distant clades, so filtering on it selects whole clades and collapses the deep structure, which is exactly what MPD sees. The spacing rule here down-weights each candidate the closer it sits in the tree to its nearest resident, so it acts on nearest neighbours, which is exactly what MNTD sees; a rule that spaced species on the trait instead would reach the tree only through the trait, and need not leave the same mark on NTI. The two indices are not redundant, and a paper reporting only one of them has answered only one of the two questions.

Two more things in that table deserve a mention. The filtered NRI has a spread of 3.76, far wider than the other rules, because each filtered community has its own randomly placed trait optimum, and how clustered the community comes out depends heavily on where that optimum falls. And the spaced rule produces a very tight distribution, because the draw that builds it is weighted so heavily towards species far from every resident that the index barely moves, even though the species picked differ from one community to the next. Effect size and its variability are separate facts about a process.

The honest limit

The number that came out of ses() is not “how clustered this community is”. It is “how clustered this community is relative to a uniform draw from these 200 species”. Change the 200 and the answer changes, and there is no data-internal test that tells you which 200 is right. Everything downstream inherits this. The tree inherits it too: swap the tree and the distances move.

The method is worth using anyway, because the alternative is reading raw MNTDs across plots of unequal richness, which measures nothing at all. But the standardised number is a statement about a comparison you chose, not a discovery about nature. The checking post makes that concrete by moving the pool and watching the conclusion move with it.

Where to go next

The obvious next question is the one this post kept deferring: which null? Drawing uniformly from the pool is not the only option, and the alternatives ask genuinely different questions rather than being better or worse approximations of one question. That comes next. After that, the same tree and the same distance matrix answer a between-community question instead of a within-community one.

References

Webb 2000 American Naturalist 156(2):145-155 (10.1086/303378)
Webb, Ackerly, McPeek, Donoghue 2002 Annual Review of Ecology and Systematics 33:475-505 (10.1146/annurev.ecolsys.33.010802.150448)
Webb, Ackerly, Kembel 2008 Bioinformatics 24:2098-2100 (10.1093/bioinformatics/btn358)
Kembel, Cowan, Helmus, Cornwell, Morlon, Ackerly, Blomberg, Webb 2010 Bioinformatics 26(11):1463-1464 (10.1093/bioinformatics/btq166)

--- title: "Phylogenetic community structure: NRI and NTI in R" description: "Raw MPD and MNTD are confounded with species richness in a plot. Build standardised effect sizes, NRI and NTI, by hand in R and learn to read their signs." date: "2026-06-14 13:00" date-modified: "2026-10-08" categories: [R, community phylogenetics, biodiversity, ggplot2, ecology tutorial] image: thumbnail.png image-alt: "Boxplots of NRI and NTI for three community assembly rules, showing clustering, neutral and overdispersed signatures." --- *Corrected 8 October 2026: the post explained the wide spread of the filtered NRI by trait optima landing in crowded or empty parts of the pool, which the runs do not bear out, and said that the spacing rule rejects candidates close to a resident; it now puts the spread down to where each optimum falls and says that the rule down-weights close candidates.* *Corrected 1 October 2026: the post called the index calibrated on the strength of its mean and spread, pointed to the checking post for a case it does not contain, and said that limiting similarity showing up mainly in NTI follows from the metrics alone; it now notes that the skewed MNTD null puts the dashed cut-offs off their nominal rates and that the spacing rule here acts directly on distances in the tree.* You have a species list for a plot and a phylogeny for the regional flora. Are the species in that plot closer relatives than you would expect, and if so, does that say anything about how the community was assembled? The metric part of this question is easy. The interpretive part is where the whole field lives, and it starts with a fact that is easy to miss: the raw number you compute from a plot is not a property of the plot. For MNTD it is mostly a property of how many species the plot has, and even for MPD how unusual a given value is depends on the species count. This post builds the standard fix from scratch: the standardised effect size, and its two sign-flipped aliases, the Net Relatedness Index (NRI) and the Nearest Taxon Index (NTI). Everything is base R plus `ape` for the tree. ## Setup A regional pool of 200 species on a coalescent tree, scaled so the root sits at depth 1 (the largest possible distance between two tips is therefore 2). One trait evolves on the tree under Brownian motion, so close relatives resemble each other: that is the assumption the whole method rests on, and it is worth naming out loud. ```{r setup} #| message: false #| warning: false library(ape) library(ggplot2) te_canvas <- theme(plot.background = element_rect(fill = "#f5f4ee", colour = NA), panel.background = element_rect(fill = "#f5f4ee", colour = NA)) set.seed(292) S <- 200 pool_tree <- rcoal(S, tip.label = sprintf("sp%03d", 1:S)) pool_tree$edge.length <- pool_tree$edge.length / max(node.depth.edgelength(pool_tree)) pool <- pool_tree$tip.label D <- cophenetic(pool_tree)[pool, pool] trait <- rTraitCont(pool_tree, model = "BM", sigma = 1)[pool] trait <- (trait - mean(trait)) / sd(trait) ``` `cophenetic()` gives the patristic distance matrix: the summed branch length on the path between each pair of tips. Everything below reads from that matrix and nothing else. ## Two distances, two depths Two summaries dominate the literature, and they are two lines each. ```{r metrics} mpd <- function(sp) { d <- D[sp, sp]; mean(d[lower.tri(d)]) } mntd <- function(sp) { d <- D[sp, sp] diag(d) <- NA mean(apply(d, 1, min, na.rm = TRUE)) } ``` Mean pairwise distance (MPD) averages over every pair in the community. Mean nearest taxon distance (MNTD) averages, over each species, the distance to its single closest relative in the community. They are not two flavours of the same thing. MPD is dominated by the deep splits: swap one species for a member of a distant clade and every pair involving it changes. MNTD only looks at each species' nearest neighbour, so it is a measurement of the terminal structure. A community can be built from two distant clades (high MPD) with tight species clusters inside each (low MNTD). Keep that picture; it explains most of the confusing results later. ## The raw number depends on richness Here is the experiment nobody runs before reporting an MNTD. Draw communities at random from the pool at five richness levels, with no assembly rule whatsoever, and look at what the raw metrics do. ```{r confounding} neutral <- function(k) sample(pool, k) set.seed(1) kg <- rep(c(4, 8, 15, 25, 40), each = 200) raw <- t(sapply(kg, function(k) { s <- neutral(k); c(mpd(s), mntd(s)) })) tab <- data.frame( k = sort(unique(kg)), mpd_m = tapply(raw[, 1], kg, mean), mpd_s = tapply(raw[, 1], kg, sd), mntd_m = tapply(raw[, 2], kg, mean), mntd_s = tapply(raw[, 2], kg, sd) ) round(tab, 3) ``` ```{r fig-confounding} #| echo: false #| fig-cap: "Raw MPD and MNTD under purely random draws from the pool, 200 replicates per richness level. Points are replicates, the line joins the means. No assembly process is present in any of these communities." #| fig-alt: "Two panels against species richness. The MPD panel shows a flat cloud whose mean stays near 1.14 while the spread narrows sharply. The MNTD panel shows a cloud whose mean falls steeply from about 0.63 at richness 4 to about 0.08 at richness 40." #| fig-width: 8 #| fig-height: 3.6 dfc <- data.frame( k = rep(kg, 2), value = c(raw[, 1], raw[, 2]), metric = rep(c("MPD", "MNTD"), each = length(kg)) ) dfc$metric <- factor(dfc$metric, levels = c("MPD", "MNTD")) ggplot(dfc, aes(k, value)) + geom_point(colour = "#93a87f", alpha = 0.35, size = 0.9, position = position_jitter(width = 0.8, height = 0)) + stat_summary(fun = mean, geom = "line", colour = "#275139", linewidth = 0.9) + stat_summary(fun = mean, geom = "point", colour = "#275139", size = 2) + facet_wrap(~metric, scales = "free_y") + labs(x = "species richness", y = "raw value") + theme_minimal(base_size = 11) + te_canvas + theme(panel.grid.minor = element_blank(), strip.text = element_text(face = "bold", colour = "#16241d")) ``` The MNTD mean falls by a factor of `r round(tab$mntd_m[1] / tab$mntd_m[5], 2)` between richness 4 and richness 40. Nothing ecological happened. Add species to a fixed pool and each species' nearest relative simply gets nearer, because there are more candidates. Any comparison of raw MNTD between a species-poor plot and a species-rich plot is a comparison of richness wearing a phylogenetic costume. MPD looks safer at first: the mean barely moves (`r round(tab$mpd_m[1], 2)` to `r round(tab$mpd_m[5], 2)`). But look at the second column. The standard deviation shrinks by a factor of `r round(tab$mpd_s[1] / tab$mpd_s[5], 2)` over the same range. A four-species plot sitting 0.2 below the pool mean is unremarkable; a forty-species plot sitting 0.2 below it is extraordinary. The raw deviation means different things at different richness, which is the same problem in a quieter voice. ## Standardising against a null The fix is the obvious one. Ask what the metric would look like if this community, with exactly this richness, had been drawn at random from the pool, and express the observation in units of that distribution. ```{r ses} ses <- function(sp, R = 999) { k <- length(sp) obs <- c(mpd(sp), mntd(sp)) nul <- replicate(R, { s <- sample(pool, k); c(mpd(s), mntd(s)) }) z <- (obs - rowMeans(nul)) / apply(nul, 1, sd) p <- (rowSums(nul <= obs) + 1) / (R + 1) list(obs = obs, null = nul, z = z, p = p, nri = -z[1], nti = -z[2]) } ``` Three conventions are buried in those few lines, and each one has cost somebody a paragraph in review. **The sign is flipped.** NRI is defined as minus the standardised MPD, and NTI as minus the standardised MNTD. So *positive NRI means clustering*: the community is more closely related than chance, which is a *low* MPD. If you report `ses.mpd` from `picante` and then discuss it as if it were NRI, every sign in your discussion is backwards. This is the single most common slip in the method. **The p-value is a rank, not a z-test.** The `p` line counts how many null draws fall at or below the observation. That works whatever shape the null distribution has, which matters more than it sounds: the MNTD null in the figure below is skewed to the right, and for a skewed null a z-score cut-off and the rank p-value can disagree. **The null is a decision, not a default.** Here the null draws species uniformly from the pool. That is one hypothesis out of several, and it is the subject of [the next post](../null-models-for-phylogenetic-structure/). ## Three assembly rules, three signatures Now build communities under three rules and see what the index reports. Habitat filtering picks species whose trait is near a local optimum. Limiting similarity picks species that are far from whatever has already been picked. The neutral rule picks at random. ```{r rules} filtered <- function(k, opt, w = 0.3) { sample(pool, k, prob = exp(-(trait - opt)^2 / (2 * w^2))) } spaced <- function(k) { sel <- sample(pool, 1) while (length(sel) < k) { cand <- setdiff(pool, sel) sel <- c(sel, sample(cand, 1, prob = apply(D[cand, sel, drop = FALSE], 1, min)^8)) } sel } set.seed(10); cf <- filtered(15, 1.2) set.seed(11); cn <- neutral(15) set.seed(12); co <- spaced(15) set.seed(20); sf <- ses(cf) set.seed(21); sn <- ses(cn) set.seed(22); so <- ses(co) data.frame( rule = c("filtered", "neutral", "spaced"), MPD = round(c(sf$obs[1], sn$obs[1], so$obs[1]), 3), NRI = round(c(sf$nri, sn$nri, so$nri), 2), MNTD = round(c(sf$obs[2], sn$obs[2], so$obs[2]), 3), NTI = round(c(sf$nti, sn$nti, so$nti), 2) ) ``` ```{r fig-null} #| echo: false #| fig-cap: "Null distributions for the filtered community (15 species, 999 draws) with the observed value marked. The observation sits far into the left tail of both distributions, which is what phylogenetic clustering looks like before it is turned into a z-score." #| fig-alt: "Two histograms of null values in sage green. The MPD null is centred near 1.15 and the observed value 0.72 sits well outside its left tail. The MNTD null is centred near 0.20 with a right skew, and the observed value 0.03 sits at its left edge." #| fig-width: 8 #| fig-height: 3.4 dfn <- data.frame( value = c(sf$null[1, ], sf$null[2, ]), metric = rep(c("MPD", "MNTD"), each = ncol(sf$null)) ) dfn$metric <- factor(dfn$metric, levels = c("MPD", "MNTD")) dfo <- data.frame(metric = factor(c("MPD", "MNTD"), levels = c("MPD", "MNTD")), obs = sf$obs) ggplot(dfn, aes(value)) + geom_histogram(bins = 40, fill = "#93a87f", colour = "white", linewidth = 0.15) + geom_vline(data = dfo, aes(xintercept = obs), colour = "#b5534e", linewidth = 0.9) + facet_wrap(~metric, scales = "free") + labs(x = "value under the null", y = "null draws") + theme_minimal(base_size = 11) + te_canvas + theme(panel.grid.minor = element_blank(), strip.text = element_text(face = "bold", colour = "#16241d")) ``` The filtered community, which really was assembled by filtering, lands at NRI `r sprintf("%+.2f", sf$nri)` with a rank p of `r sprintf("%.3f", sf$p[1])`. The neutral community lands at NRI `r sprintf("%+.2f", sn$nri)`, p `r sprintf("%.3f", sn$p[1])`, against a truth of no structure at all. The spaced community lands at NTI `r sprintf("%+.2f", so$nti)`: overdispersed at the tips. Repeat over 100 communities per rule to see how reliable that is. ```{r replicate} rep_ses <- function(gen, n = 100, R = 199) { t(replicate(n, { sp <- gen(15) obs <- c(mpd(sp), mntd(sp)) nul <- replicate(R, { s <- sample(pool, 15); c(mpd(s), mntd(s)) }) -(obs - rowMeans(nul)) / apply(nul, 1, sd) })) } set.seed(30); A <- rep_ses(function(k) filtered(k, rnorm(1))) set.seed(31); B <- rep_ses(neutral) set.seed(32); C <- rep_ses(spaced) data.frame( rule = c("filtered", "neutral", "spaced"), NRI = round(c(mean(A[, 1]), mean(B[, 1]), mean(C[, 1])), 2), NRI_sd = round(c(sd(A[, 1]), sd(B[, 1]), sd(C[, 1])), 2), NTI = round(c(mean(A[, 2]), mean(B[, 2]), mean(C[, 2])), 2), NTI_sd = round(c(sd(A[, 2]), sd(B[, 2]), sd(C[, 2])), 2) ) ``` ```{r fig-rules} #| echo: false #| fig-cap: "NRI and NTI over 100 communities per assembly rule, 15 species each. A solid line marks zero and dashed lines mark the conventional two-tailed cut-offs. Positive is clustered, negative is overdispersed." #| fig-alt: "Boxplots of NRI and NTI for three rules. Filtered communities sit well above zero on NRI with a long upper tail and modestly above zero on NTI. Neutral communities straddle zero on both. Spaced communities sit below zero on both, much further below on NTI." #| fig-width: 8 #| fig-height: 3.8 dfr <- data.frame( value = c(A[, 1], B[, 1], C[, 1], A[, 2], B[, 2], C[, 2]), rule = rep(rep(c("filtered", "neutral", "spaced"), each = nrow(A)), 2), index = rep(c("NRI (from MPD)", "NTI (from MNTD)"), each = 3 * nrow(A)) ) dfr$rule <- factor(dfr$rule, levels = c("filtered", "neutral", "spaced")) ggplot(dfr, aes(rule, value, fill = rule)) + geom_hline(yintercept = 0, colour = "#16241d", linewidth = 0.3) + geom_hline(yintercept = c(-1.96, 1.96), linetype = 2, colour = "#5d6b61", linewidth = 0.3) + geom_boxplot(outlier.size = 0.7, outlier.colour = "#5d6b61", width = 0.6, colour = "#16241d", linewidth = 0.35) + scale_fill_manual(values = c(filtered = "#275139", neutral = "#93a87f", spaced = "#cda23f"), guide = "none") + facet_wrap(~index) + labs(x = NULL, y = "index value") + theme_minimal(base_size = 11) + te_canvas + theme(panel.grid.minor = element_blank(), strip.text = element_text(face = "bold", colour = "#16241d")) ``` The neutral row is the one to read first. Its NRI averages `r sprintf("%+.2f", mean(B[, 1]))` with a spread of `r sprintf("%.2f", sd(B[, 1]))`, and its NTI averages `r sprintf("%+.2f", mean(B[, 2]))` with a spread of `r sprintf("%.2f", sd(B[, 2]))`. Centred on zero, spread of about one: that is what a standardised effect size is supposed to do when the data really were generated by the null. Centre and spread are not the whole story, though: the MNTD null is skewed to the right (the figure above), so under the null NTI almost never reaches the upper dashed line and passes the lower one a little more often than a normal curve would. The rank p-value, which assumes no shape, is the one to test with. ## Why NRI and NTI disagree Now the asymmetry. Habitat filtering shows up mainly in NRI (`r sprintf("%+.2f", mean(A[, 1]))`, against `r sprintf("%+.2f", mean(A[, 2]))` for NTI). Limiting similarity shows up mainly in NTI (`r sprintf("%+.2f", mean(C[, 2]))`, against `r sprintf("%+.2f", mean(C[, 1]))` for NRI). Part of that is built into these simulations. A conserved trait varies most between distant clades, so filtering on it selects whole clades and collapses the deep structure, which is exactly what MPD sees. The spacing rule here down-weights each candidate the closer it sits in the tree to its nearest resident, so it acts on nearest neighbours, which is exactly what MNTD sees; a rule that spaced species on the trait instead would reach the tree only through the trait, and need not leave the same mark on NTI. The two indices are not redundant, and a paper reporting only one of them has answered only one of the two questions. Two more things in that table deserve a mention. The filtered NRI has a spread of `r sprintf("%.2f", sd(A[, 1]))`, far wider than the other rules, because each filtered community has its own randomly placed trait optimum, and how clustered the community comes out depends heavily on where that optimum falls. And the spaced rule produces a very tight distribution, because the draw that builds it is weighted so heavily towards species far from every resident that the index barely moves, even though the species picked differ from one community to the next. Effect size and its variability are separate facts about a process. ## The honest limit The number that came out of `ses()` is not "how clustered this community is". It is "how clustered this community is *relative to a uniform draw from these 200 species*". Change the 200 and the answer changes, and there is no data-internal test that tells you which 200 is right. Everything downstream inherits this. The tree inherits it too: swap the tree and the distances move. The method is worth using anyway, because the alternative is reading raw MNTDs across plots of unequal richness, which measures nothing at all. But the standardised number is a statement about a comparison you chose, not a discovery about nature. The [checking post](../checking-a-community-phylogenetics-analysis/) makes that concrete by moving the pool and watching the conclusion move with it. ## Where to go next The obvious next question is the one this post kept deferring: which null? Drawing uniformly from the pool is not the only option, and the alternatives ask genuinely different questions rather than being better or worse approximations of one question. That comes next. After that, the same tree and the same distance matrix answer a between-community question instead of a within-community one. ## References - Webb 2000 American Naturalist 156(2):145-155 (10.1086/303378) - Webb, Ackerly, McPeek, Donoghue 2002 Annual Review of Ecology and Systematics 33:475-505 (10.1146/annurev.ecolsys.33.010802.150448) - Webb, Ackerly, Kembel 2008 Bioinformatics 24:2098-2100 (10.1093/bioinformatics/btn358) - Kembel, Cowan, Helmus, Cornwell, Morlon, Ackerly, Blomberg, Webb 2010 Bioinformatics 26(11):1463-1464 (10.1093/bioinformatics/btq166) ## Related tutorials - [Null models for phylogenetic structure](../null-models-for-phylogenetic-structure/) - [Phylogenetic beta diversity: UniFrac and PhyloSor](../phylogenetic-beta-diversity/) - [Checking a community phylogenetics analysis](../checking-a-community-phylogenetics-analysis/) - [Phylogenetic diversity with picante](../phylogenetic-diversity-picante/)