How Many Samples Before Your Species Count Means Anything

The hundredth species is a rumour

Here is the number that started this. Walk into a plot holding exactly one hundred species, identify five hundred individuals with all the care you have, and you will still leave the plot with about eighty-six names on the sheet. The fault is not carelessness. Fourteen of those species sit rare enough that five hundred draws went past them without one individual turning up, and no amount of care at the frame changes that.

Ecology has known this since the 1940s. A small industry of corrections has grown up around the problem since [13][14], and the most cited of them, Chao1, is nearly forty years old [9] and turns up in a remarkable share of the biodiversity papers published in any given month. We wanted its delivered performance rather than its advertised one, which is a different question and a harder one to ask. No field site can settle that question. Field sites do not come labelled.

So we made the label first.

What we did, and where we have never been

None of us has ever run a transect. The club meets in a room with computers in it, we have no field site, no quadrat frame, no pooter, and the nearest thing we own to a taxonomic key is a library card. We put that at the top rather than the bottom. It decides what the paper may claim, and the honest limit is narrower than a reader might assume.

What we can do is build a community whose species list we wrote ourselves, sample it the way a fieldworker would, and watch what the estimators say. The community is a vector of one hundred relative abundances. A survey is ten quadrats, each an independent multinomial draw of the same size summing to the total effort \(n\), and summing across those quadrats gives the abundance counts. Thresholding at zero gives incidence counts, which is what Chao2 and the incidence jackknives need. Both come off the same draw, so every estimator sees the same survey.

Two families of abundance distribution. Log-normal, which is what Preston found when he plotted moth catches on a logarithmic axis in 1948 [2], with the spread parameter \(\sigma\) running from 0.5 to 2.0. And MacArthur's broken stick [3], which you get by snapping a unit stick at ninety-nine uniformly placed points, and which sits close to the most even arrangement anyone has seriously proposed for a real assemblage. Evenness across the five configurations: 0.973, 0.909, 0.895, 0.780 and 0.652 on Pielou's \(J\).

The computation is the experiment. Nothing in this article was observed. When we write that Chao1 was short by 4.58 species, we mean that a program we wrote printed that number, and the program is in the repository with its output beside it.

The arithmetic

Sample \(n\) individuals from a community with relative abundances \(p_1, \ldots, p_S\). A given species is missed with probability \((1-p_i)^n\). The complement gives the chance of finding it. Sum those complements across all one hundred species and you have the expected count.

The curve has no free parameters. Put \(p_i = 0.002\) and \(n = 500\) into it. A single such species is found with probability \(1 - 0.998^{500} = 0.632\). Drop the abundance to \(p_i = 0.0002\) and that probability falls to 0.095. Rare times cheap. The species you miss are the ones whose abundance times your effort is small, and no cleverness downstream recovers them.

Effort enters only through the product \(np_i\). Doubling the effort and halving every abundance leaves the curve exactly where it was, which is the one piece of good news in this section. Richer communities therefore cost more for the same result, because the same individuals get split across more species and every \(p_i\) shrinks accordingly, which pushes the whole curve rightward along the effort axis.

Rarefying a sample to a common effort starts with Sanders in 1968 [4]. His arithmetic was wrong. It inflated the curve, and Hurlbert fixed it three years later [5] with the hypergeometric version everybody now uses, though the 1968 original still picks up citations from people who have not checked.

Figure 1 has the five curves. Circles are simulated means; lines are the formula. At \(n = 50\) the log-normal \(\sigma = 2.0\) community yields 19.95 species out of a hundred, and at \(n = 50{,}000\) it yields 98.25, which is 1.75 species short. A thousand-fold increase in effort bought 78 species and still did not finish the community off.

Six formulas and what each one promises

All six work the same way. Count the species you saw exactly once, count the ones you saw exactly twice, and let the ratio stand in for the species you saw no times at all. Singletons far outnumbering doubletons means a lot of undetected material out there. No singletons left, and you are done, though that state arrives later than anybody hopes.

Chao1 came out of a 1984 paper on estimating the number of classes in a population [9], and Chao framed it there as a lower bound rather than as a point estimate. We use the bias-corrected form [10]:

Chao2 runs the same idea on incidence data, counting the species found in exactly one and in exactly two of the \(T\) sampling units [10], so it never looks at how many individuals of anything you found. ACE takes another route. It estimates your sample's coverage of the community, then corrects for rare-species variability [11][12]. The jackknives are older, out of capture-recapture [6] and adapted to species counts by Heltshe and Forrester in 1983 [7], and the first-order version is \(S_{\text{obs}} + f_1(n-1)/n\). For any reasonable \(n\) that is the singleton count, added on. The second-order form subtracts a doubleton term, and neither version pretends to know anything about how abundance is distributed.

Does the sampler agree with the algebra?

Three checks, and one of them nearly convinced us the code was broken.

The first compares the simulated accumulation curve with the closed form above, paired one community at a time so that community-to-community variation cancels, which is the only way to see a discrepancy of a hundredth of a species through the noise of six thousand redraws. Across 50 cells the largest disagreement is 0.087 species out of 100, the spread of \(z\) scores is 1.01 against the 1.00 a correct sampler should give, and the signed mean is \(-0.09\). No cell sits past three standard errors.

It did not start that way. The first version of that table had a cell at \(z = 12.37\) and two more past 3, and we spent a while looking for a bug that was not there. Every bad cell sat where the simulation found all 100 species every time. The simulated mean was exactly 100.000. The analytic value was 99.99994. The difference was six hundred-thousandths of a species, and the \(z\) score was enormous only because the empirical variance of a column of identical numbers is zero, so we were dividing a tiny number by nothing.

The second check is Hurlbert's within-sample rarefaction [5], the hypergeometric version, which asks how many species a random subsample of \(k\) individuals from your collection would have held. A different calculation from the first, and one that draws without replacement. Against 20,000 direct subsamples of a 2,000-individual collection, at eight subsample sizes, the largest disagreement runs to 0.71 standard errors, which is noise rather than disagreement. At \(k = n\) both agree on 99.0000 exactly, which they must.

The third check is the boring one. Sample four communities with four million individuals each. Every species turns up dozens of times, no singletons survive, and all six estimators collapse to the raw count, returning 100.000 against a true value of 100.000, a difference of zero. A weak test. It proves the formulas are typed correctly, and a test that only catches typing errors is still worth running before you trust a table.

Notes from the club table

What the estimators actually delivered

Figure 2 answers the question we came in with. Every curve starts deep in negative territory. The raw count is 67.6 species short at \(n = 50\). Chao1 recovers most of that. Not all: 31.4 short at the same effort, 4.58 short at \(n = 500\). ACE does slightly better and Chao2 slightly worse, which is the ordering that holds across most of the grid. The two jackknives cross zero and overshoot somewhere between 100 and 500 individuals, which is a different failure mode and, for a conservation report, a more dangerous one.

Now the part we enjoyed more than we should have. Chao1 is the most cited richness estimator in ecology, and in our grid the second-order jackknife beats it on effort-to-target in four configurations out of five, by a factor between 2.5 and 3.3. On the broken-stick community Chao1 needs 406 individuals to reach the 10 percent band. The jackknife needs 124, and that ratio holds roughly across the other four communities.

None of this is news to the people who test estimators for a living. Walther and Moore's review of estimator performance [15] reports the same low bias for the Chao family and the same tendency of the jackknives to overshoot, across dozens of studies and real datasets. A simulation adds one thing: the truth sits on the page beside the estimate, rather than being inferred from the same data that produced it.

Nothing about the gap is mysterious once you check what each formula was built to return. Chao1 is a lower bound [9]. A bound that is working correctly sits below the truth, and asking it to be unbiased is asking it to stop being a bound.

Figure 3 puts the same six estimators at a single effort against evenness. The lines fan out as the community gets less even, and the fan is wide enough that the estimator you picked matters less than the community you picked it for. At \(J = 0.973\) everything sits inside the 10 percent band. At \(J = 0.652\) nothing does. The best of them is still 20.7 species short, which is a fifth of the community gone missing.

The bill, in individuals

The table is the whole study compressed. Relative bias at four efforts, root mean square error at 500 individuals, and the effort at which each estimator's bias first falls inside 10 percent of the truth, interpolated in \(\log n\) between rungs of the ladder. Because the truth is 100 species in every row, relative bias in percent and bias in species are the same number, which is one small convenience of building the community yourself. RMSE is in species. The bolded estimator in each block is the one that reaches the target first.

Read down the last column. On a nearly even community the whole question is cheap: ACE gets there at 50 individuals, everything else by 96, and even the uncorrected count manages it at 342. Push evenness down to \(J = 0.652\) and the same target runs between 1,258 and 7,829 individuals, depending on which formula you trust, which is a sampling programme rather than a morning. An afternoon, against a summer.

Two numbers from that figure are worth holding on to. Across rows, going from the most even community to the least multiplies the effort requirement by 15 for the second-order jackknife and by 73 for ACE. Within a row, the five estimators differ from one another by a factor of 1.9 to 3.6, and all of them beat the uncorrected count by a factor of 6 to 7. So the community you are counting matters more than the formula you count it with. Both matter less than whether you corrected at all, which is the one decision nobody argues about.

The strongest objection we can make

Our quadrats are too well behaved. Not a small point, and arguably fatal to the incidence-based half of this study.

Real organisms are clumped. Plants have seed shadows and insects sit on their host plants. Nothing in a real assemblage spreads itself evenly across a field. A real quadrat is not an independent draw from the whole community; it is a draw from a local patch that looks a lot like the quadrat beside it. Chao2 and the incidence jackknives exist because incidence data resists the problem of abundance being hard to count, and their behaviour under aggregation is a known and studied problem [8]. Our model has no space in it at all, no patches and no gradient from one corner of the plot to the other. Every quadrat is an independent multinomial draw with the same parameter vector, which is the friendliest possible world for an incidence estimator, and Chao2 still finished behind Chao1 in every configuration.

You could read that two ways. The generous reading is that Chao2's poor showing here is a floor, and under aggregation it does relatively better because incidence stops throwing away information that abundance was mismeasuring anyway. The harsh reading is that we have shown Chao2 to be dominated by Chao1 in the one case where we can compute both exactly, and that everything beyond that is speculation. We think the generous reading is probably right. We cannot demonstrate it. Demonstrating it needs a spatial model we did not build, and until somebody builds one the incidence half of this paper should be read as provisional.

A second objection, from the same direction. Every individual in our model is identified correctly. Some fraction of apparent singletons are misidentifications, or juveniles keyed to the wrong genus. Those inflate \(f_1\), and \(f_1\) is what every formula here runs on. Every estimator in this paper is increasing in \(f_1\), so every one of them is inflated by identification error, and Chao1, whose correction term goes as \(f_1^2\), is inflated hardest. A Chao1 estimate from a difficult taxonomic group is biased low by undersampling and biased high by misidentification at the same time, and nothing in our results tells you which one wins.

Third, and smallest. We fixed richness at 100 species in every configuration, and estimator performance depends on richness as well as on evenness, which we did not vary and cannot speak to. A hundred species makes a decent insect trap and a very poor tropical tree plot.

Where another choice would have changed the answer

We chose ten quadrats. That choice is not neutral for Chao2, which reads its singletons off the incidence table rather than the abundance table, and so depends on how the same individuals were bagged. Holding effort at 1,000 individuals and varying only the split, across 1,000 replicates each:

Finer units give a less biased and noisier answer. Coarser units give the reverse. The bias swings by 63 percent across that range on identical data, and the RMSE is worst at the finest split, so there is no setting that wins on both. Everything this paper says about Chao2 is conditional on T = 10. In a real survey the number of quadrats is usually decided by how far a person can walk before lunch, which is not a statistical criterion.

The abundance distribution is the other lever, and a larger one. Log-normal is a defensible default with sixty years of support behind it [2], but a log-series community of the kind Fisher, Corbet and Williams fitted to Malayan butterflies [1] has a fatter tail of rarities at the same richness, and a Zipf-Mandelbrot community fatter still. Every effort figure in this article would rise. The ordering of the estimators might not survive either, since the jackknives' advantage comes from overcorrecting, and a heavier tail means there is more to correct for.

One choice turned out not to matter. Should we draw a fresh community for every replicate, or reuse one? We drew fresh. Our standard errors therefore include community-to-community variation as well as the sampling noise. Repeating the convergence cell as an independent 150,000-replicate run gave means agreeing with the main sweep to within 1.83 standard errors on all six estimators, the largest absolute discrepancy being 0.23 species.

What to do on Monday, and how to reproduce this

Suppose you have one afternoon and a pooter, and you want a number you can defend in a report. Here is what our grid says.

Report the raw count first, always, and report the effort beside it. A richness figure without an effort figure is uninterpretable, and Gotelli and Colwell made this argument twenty-odd years ago in terms nobody has improved on [13]. Then report an estimator, and say which one and why. If your sample looks even, meaning your commonest species is not more than about twenty times your median species, almost anything works and the cheapest correct answer is ACE. If it looks uneven, nothing works well at the effort you can afford, and the honest move is to say so rather than to pick the estimator with the friendliest number.

The modern treatment of all this folds richness into a family of diversity measures and puts rarefaction and extrapolation on a common footing [16], which is a better place to start than any single estimator. It also makes the trade we have been describing explicit rather than implicit. We did not test it. A framework is not a formula, and our grid needed one number out.

Report singletons and doubletons. Two integers, \(f_1\) and \(f_2\). Those two numbers let any reader recompute every estimator in this paper straight off your table, which is more than most published richness figures allow, and it costs you one line. A large \(f_1\) relative to your species list means the survey is not finished. No formula will finish it for you, whatever the estimate on the bottom line says.

And stop correcting once the correction stops helping. When \(f_1\) falls to two or three species, the estimator is adding variance to a count that is already close enough, and the variance you add is not small. Our grid puts that crossover at around \(n = 2{,}000\) to \(5{,}000\) for the middle communities, where the raw count starts beating Chao1 on root mean square error. Strange advice, we know. The numbers are in the table, and you can check the crossover for your own community.

Rerun it from analysis/:

python sampling-effort.py > sampling-effort-output.txt

Python 3.12 and numpy, nothing else. On the club laptop, a four-core machine from 2021, it takes 4 minutes 48 seconds of wall time, of which about two and a half minutes is the main sweep and 82 seconds is the 150,000-replicate convergence run. Everything is seeded from numpy.random.SeedSequence(20241104), spawned per configuration. The output is deterministic and each configuration is independent of the others. If your numbers differ from ours by more than the printed standard errors, one of us has a problem and we would like to know about it. The interactive model runs a smaller version of the same machinery in your browser, and its default settings reproduce the headline numbers in this article, if you would rather poke at it than read it.

Two companion pieces in this issue come at the same problem from other angles. Absence of Evidence asks how many blank surveys it takes before you can call a species gone, which is the single-species version of the question we have been asking about whole assemblages. A Population That Grows on Average takes the other route by which chance removes a species from a list, which is killing it rather than hiding it.

References