Checking a diversity estimate

A richness estimate is only as sound as its rarest counts. Check sample coverage, singleton sensitivity and interval coverage in base R, and see what the number can promise.

The earlier posts in this cluster produced richness and diversity estimates that beat the naive observed counts. This one turns them over and asks how far to trust them. The theme running through the whole cluster is that the unseen species are not recoverable from the seen without an assumption, so the honest question is not “what is the estimate” but “what does the estimate depend on, and how much”. Three checks answer that: how far you are extrapolating, how much the estimate rides on a fragile count, and whether its confidence interval means what it seems to.

library(ggplot2); library(dplyr); library(tidyr)
te <- c(forest = "#2f5d50", moss = "#7a9b76", rust = "#b5651d", gold = "#c9a227",
        slate = "#3d4b53", sand = "#d9cbb2", sky = "#5b8aa6", brick = "#8c3b2e")
theme_te <- function(base_size = 12) {
  theme_minimal(base_size = base_size) +
    theme(panel.grid.minor = element_blank(),
          panel.grid.major = element_line(colour = "#e7e1d5", linewidth = 0.3),
          axis.title = element_text(colour = "#3d4b53"),
          axis.text  = element_text(colour = "#5c6670"),
          plot.title = element_text(colour = "#2f3a40", face = "bold", size = base_size + 1),
          plot.subtitle = element_text(colour = "#5c6670", size = base_size - 1),
          legend.position = "bottom",
          legend.title = element_text(colour = "#3d4b53"),
          plot.background  = element_rect(fill = "white", colour = NA),
          panel.background = element_rect(fill = "white", colour = NA))
}
covCJ <- function(x) { n <- sum(x); f1 <- sum(x == 1); f2 <- sum(x == 2)
  if (f2 == 0) f2 <- f1 * (f1 - 1) / 2
  1 - (f1 / n) * ((n - 1) * f1 / ((n - 1) * f1 + 2 * f2)) }
chao1 <- function(x) { n <- sum(x); f1 <- sum(x == 1); f2 <- sum(x == 2)
  length(x) + (n - 1) / n * f1 * (f1 - 1) / (2 * (f2 + 1)) }

set.seed(4193)
S_true <- 180
rel <- exp(rnorm(S_true, 0, 1.6)); rel <- rel / sum(rel)

Check one: how far are you extrapolating?

Sample coverage is the plainest reliability check. It says what fraction of the community’s individuals belong to species you have already seen, so the coverage deficit is roughly how much of the estimate is projection rather than observation. Compare a light effort with a heavy one on the same community.

n_low <- 250; n_high <- 1500
x_lo <- as.vector(rmultinom(1, n_low,  rel)); x_lo <- x_lo[x_lo > 0]
x_hi <- as.vector(rmultinom(1, n_high, rel)); x_hi <- x_hi[x_hi > 0]
data.frame(effort = c("low (n = 250)", "high (n = 1500)"),
           S_obs   = c(length(x_lo), length(x_hi)),
           coverage = round(c(covCJ(x_lo), covCJ(x_hi)), 4),
           Chao1    = round(c(chao1(x_lo), chao1(x_hi)), 1),
           unseen   = round(c(chao1(x_lo) - length(x_lo), chao1(x_hi) - length(x_hi)), 1))

           effort S_obs coverage Chao1 unseen
1   low (n = 250)    73   0.8764 103.9   30.9
2 high (n = 1500)   134   0.9834 147.0   13.0

The light sample reaches 0.876 coverage and asks Chao1 to invent 30.9 unseen species on top of the 73 observed. The heavy sample reaches 0.983 coverage, so its estimate leans on only 13 projected species and sits much closer to the true 180. Low coverage is a warning that the estimate is mostly extrapolation, and extrapolation is where assumptions do the work.

Check two: does the estimate ride on the singletons?

Chao1 is built from singletons and doubletons, so the singleton count is load-bearing. That matters because singletons are exactly the counts most easily corrupted: a misidentification, a contaminant, a sequencing error in a metabarcoding run all masquerade as a species seen once. Vary the singleton count and watch the estimate move.

x <- x_lo; So <- length(x); f1 <- sum(x == 1); f2 <- sum(x == 2); nn <- sum(x)
chao1_f1 <- function(f1v) So + (nn - 1) / nn * f1v * (f1v - 1) / (2 * (f2 + 1))
c(f1_observed = f1, f2 = f2,
  Chao1_observed = round(chao1_f1(f1), 1),
  Chao1_if_20pct_spurious = round(chao1_f1(round(f1 * 0.8)), 1),
  Chao1_if_20pct_more     = round(chao1_f1(round(f1 * 1.2)), 1),
  slope_per_singleton = round(chao1_f1(f1 + 1) - chao1_f1(f1), 2))

            f1_observed                      f2          Chao1_observed 
                  31.00                   14.00                  103.90 
Chao1_if_20pct_spurious     Chao1_if_20pct_more     slope_per_singleton 
                  92.90                  117.20                    2.06 

With 31 singletons the estimate is 103.9. Treat a fifth of them as spurious and it drops to 92.9; add a fifth and it climbs to 117.2. Near the observed value each singleton is worth about 2.06 species of estimated richness. An estimate this sensitive to a fragile count deserves a hard look at where the singletons came from; corrected estimators that down-weight spurious singletons exist for exactly this reason.

f1grid <- seq(max(2, round(f1 * 0.6)), round(f1 * 1.4), by = 1)
fb <- data.frame(f1 = f1grid, chao1 = sapply(f1grid, chao1_f1))
marks <- data.frame(f1 = c(round(f1 * 0.8), f1, round(f1 * 1.2)),
                    chao1 = chao1_f1(c(round(f1 * 0.8), f1, round(f1 * 1.2))))
ggplot(fb, aes(f1, chao1)) +
  geom_line(colour = te["forest"], linewidth = 0.9) +
  geom_point(data = marks, colour = te["rust"], size = 3) +
  geom_text(data = marks, aes(label = round(chao1)), vjust = -1.1, size = 3, colour = te["slate"]) +
  labs(title = "Chao1 rides on the singleton count",
       subtitle = "A 20% change in singletons (miscounts, contamination) moves the estimate by about 12 species",
       x = "number of singletons (species seen once)", y = "Chao1 richness estimate") +
  theme_te()

Figure 1: Chao1 as a function of the singleton count, holding the other counts fixed, with the observed value and plus or minus twenty per cent marked.

Check three: what does the confidence interval cover?

Chao’s variance formula gives an interval, but the symmetric version misbehaves because the estimator’s sampling distribution is right-skewed. The log-transformed interval fixes the shape and keeps the lower limit above the observed richness. The deeper point is what either interval is an interval for.

chao1_ci <- function(x, level = 0.95) { n <- sum(x); So <- length(x)
  f1 <- sum(x == 1); f2 <- sum(x == 2); Sc <- chao1(x); D <- Sc - So
  r <- f1 / f2; v <- f2 * (0.5 * r^2 + r^3 + 0.25 * r^4)     # Chao 1987 variance
  z <- qnorm(1 - (1 - level) / 2)
  wald <- c(Sc - z * sqrt(v), Sc + z * sqrt(v))
  K <- exp(z * sqrt(log(1 + v / D^2)))                        # log-transformed
  list(Sc = Sc, wald = wald, logci = c(So + D / K, So + D * K)) }
ci <- chao1_ci(x_lo)
c(Chao1 = round(ci$Sc, 1),
  Wald_lo = round(ci$wald[1], 1), Wald_hi = round(ci$wald[2], 1),
  log_lo = round(ci$logci[1], 1), log_hi = round(ci$logci[2], 1))

  Chao1 Wald_lo Wald_hi  log_lo  log_hi 
  103.9    71.6   136.1    84.6   155.2 

Now check both intervals by simulation, against two targets: the estimator’s own expectation, and the true richness.

set.seed(5193); R <- 2000
Sc <- wl <- wh <- ll <- lh <- numeric(R)
for (r in 1:R) { xx <- as.vector(rmultinom(1, n_low, rel)); xx <- xx[xx > 0]
  ci <- chao1_ci(xx); Sc[r] <- ci$Sc
  wl[r] <- ci$wald[1]; wh[r] <- ci$wald[2]; ll[r] <- ci$logci[1]; lh[r] <- ci$logci[2] }
E_chao1 <- mean(Sc)
data.frame(target = c("E[Chao1] (estimator)", "true richness"),
  Wald = round(c(mean(wl <= E_chao1 & wh >= E_chao1), mean(wl <= S_true & wh >= S_true)), 3),
  log  = round(c(mean(ll <= E_chao1 & lh >= E_chao1), mean(ll <= S_true & lh >= S_true)), 3))

                target  Wald   log
1 E[Chao1] (estimator) 0.899 0.962
2        true richness 0.230 0.427

fc <- data.frame(
  method = factor(rep(c("symmetric Wald", "log-transformed"), each = 2),
                  levels = c("symmetric Wald", "log-transformed")),
  target = rep(c("E[Chao1] (estimator target)", "true richness"), 2),
  coverage = c(mean(wl <= E_chao1 & wh >= E_chao1), mean(wl <= S_true & wh >= S_true),
               mean(ll <= E_chao1 & lh >= E_chao1), mean(ll <= S_true & lh >= S_true)))
ggplot(fc, aes(method, coverage, fill = target)) +
  geom_hline(yintercept = 0.95, linetype = "dashed", colour = te["brick"], linewidth = 0.5) +
  geom_col(position = position_dodge(width = 0.7), width = 0.6) +
  geom_text(aes(label = sprintf("%.2f", coverage)), position = position_dodge(width = 0.7),
            vjust = -0.4, size = 3, colour = te["slate"]) +
  scale_fill_manual(values = c("E[Chao1] (estimator target)" = unname(te["forest"]),
                               "true richness" = unname(te["sky"])), name = NULL) +
  scale_y_continuous(limits = c(0, 1.05)) +
  labs(title = "The interval covers the lower bound, not the truth",
       subtitle = "n = 250, 2000 samples; dashed line 0.95: log CI reaches the estimator target, neither reaches true richness 180",
       x = NULL, y = "interval coverage") +
  theme_te()

Figure 2: Interval coverage of the estimator target and of the true richness, for the symmetric and log-transformed intervals.

The mean Chao1 across samples is 116.8, and the log-transformed interval covers that at about 0.96, close to the nominal 0.95, while the symmetric interval falls to 0.9. So the log interval is the right construction. But watch the second target: neither interval covers the true richness of 180, at 0.23 and 0.43. That is not a broken interval. It is Chao1 being a lower bound: the interval brackets the estimator, and the estimator sits below the truth by construction.

What the number can promise

Read the three checks together and the estimate stops being a single figure and becomes a claim with conditions attached. It is a lower bound, its interval is a confidence statement about that lower bound, and both depend on sample coverage and on a singleton count that field or laboratory error can distort. This is the same shape as the honest limits elsewhere on the blog: the discovery set that depends on the chosen error rate, the effect whose sign depends on an unmeasured reference frame. The estimate is a function of what you assume about the unseen, not a fact about the community, and the checks make the assumptions visible rather than pretending they are not there.

References

Chao A 1987 Biometrics 43(4):783-791 (10.2307/2531532)

Colwell R K, Chao A, Gotelli N J, Lin S-Y, Mao C X, Chazdon R L, Longino J T 2012 Journal of Plant Ecology 5(1):3-21 (10.1093/jpe/rtr044)

Chao A, Jost L 2012 Ecology 93(12):2533-2547 (10.1890/11-1952.1)

Walther B A, Moore J L 2005 Ecography 28(6):815-829 (10.1111/j.2005.0906-7590.04112.x)

Chiu C-H, Chao A 2016 PeerJ 4:e1634 (10.7717/peerj.1634)

Related tutorials

• Estimating species richness beyond your sample
• Coverage-based rarefaction and extrapolation
• Estimating diversity with Hill numbers
• Rarefaction and accumulation curves