Introduction
The incompleteness of the fossil record has been debated heavily ever since a groundbreaking study (Raup Reference Raup1972) brought analytical methods into the arena. However, there is still no real agreement on the very basic question of which method, or even which family of methods, should be used to compensate for incompleteness when counting taxa. A survey of all options, of which there are many (Colwell and Coddington Reference Colwell and Coddington1994; Bunge et al. Reference Bunge, Willis and Walsh2014), would be impractically complex and lengthy. Instead, this paper focuses on just four methods: shareholder quorum subsampling (Alroy Reference Alroy2009, Reference Alroy2010a,Reference Alroyb), which is now widely used in paleobiology and ecology; a new extrapolation method called squares (Alroy Reference Alroy2018); another called the corrected first-order jackknife (cJ1), which belongs to a family of relatively accurate and precise equations that are backed by solid and simple proofs; and another called true richness estimated using a Poisson sampling model (TRiPS) (Starrfelt and Liow Reference Starrfelt and Liow2016). The analysis will elaborate on a recent comparison of quorum subsampling and TRiPS (Close et al. Reference Close, Evers, Alroy and Butler2018) by using different tests and different data, and also by considering squares and cJ1.
The four approaches are very different. Quorum subsampling, like multiton subsampling (Alroy Reference Alroy2017) and the much older approach of rarefaction (Sanders Reference Sanders1968; Hurlbert Reference Hurlbert1971), seeks only to provide unbiased relative estimates of species richness, as opposed to estimating the full number of taxa in each sampling universe. Squares, cJ1, and TRiPS do seek to estimate true richness, but they otherwise have nothing in common, except that the latter two methods assume that the sampling process is Poisson.
The three extrapolation methods also all assume that the number of species is finite. So, for example, it makes sense to apply extrapolators when real-world distributions follow the log-normal or double geometric distributions, which they often do (Alroy Reference Alroy2015). Richness would not be defined if real-world distributions followed the log series or Zipf distributions, in which case no richness interpolator or extrapolator would be of any use. However, global richness is finite by definition, the log series highly resembles the log normal when sampling is poor, and the Zipf fits hardly any real-world data well (Alroy Reference Alroy2015). Thus, extrapolation is a reasonable thing to do.
Putting all that aside, differences among the methods are worth emphasizing. Squares and cJ1 are simple nonparametric indices based on very different reasoning, neither of which bears any relationship to TRiPS. Squares basically assumes that richness estimates should be increased when the most common species are unusually common and many species are each represented by a single observation (i.e., are singletons). Although cJ1 also assumes that estimates should be increased when there are many singletons, it pays attention to the overall sample size instead of computing some function of all the individual abundances. Both methods are designed for general use. By contrast, TRiPS is a parametric method that back-estimates a single sampling probability per time interval from the distribution of taxonomic occurrences and then scales up a time series of taxon counts by dividing observed values by the sampling probabilities.
It is important to stress that TRiPS rests strongly on the assumption that every species within a given time interval has the same sampling probability. Assuming uniformity in this way seems to be responsible for many of the method's properties. However, back-estimating sampling probabilities is not problematic by itself. As I discuss later, this idea goes back to Buzas et al. (Reference Buzas, Koch, Culver and Sohl1982) and was given a much more realistic treatment by them and by Wagner and Marcot (Reference Wagner and Marcot2013), both of whom relaxed the uniformity assumption. Meanwhile, the additional step of rescaling diversity curves that figures in TRiPS was employed by Alroy (Reference Alroy2014).
The results of the analyses presented here are straightforward. Based on their very nature, most extrapolation methods do more or less obey Hill's replication principle (Hill Reference Hill1973) in the same way that quorum subsampling does (Alroy Reference Alroy2010b): doubling the number of species in the system, known or unknown, is mathematically required to double the extrapolated richness estimates, even though this does not always occur in practice. However, unlike the other three methods discussed here, TRiPS does not meet other basic criteria for a richness estimation protocol. In particular, I demonstrate that TRiPS fails on three counts. First, its estimates seem to depend solely on the numbers of individuals and species to be found in particular samples, not on the shapes of abundance distributions. Second, it returns very different estimates when the raw data are first subjected to subsampling (in this case rarefaction: Sanders Reference Sanders1968; Hurlbert Reference Hurlbert1971; Heck et al. Reference Heck, van Belle and Simberloff1975). Third, it produces unrealistic values when subsets of inventories are summed and then compared with estimates generated on the basis of the raw, unsplit inventories. A fourth and final test shows that unlike the other methods, quorum subsampling is influenced substantially by differences in evenness, which is not a good thing. The possibility that both rarefaction and quorum subsampling have such a flaw has been suggested previously (Mitchell Reference Mitchell2015). It is noteworthy that this result rests on a new metric of evenness having low bias, although the very widely used and less reliable Pielou's J (Pielou Reference Pielou1966) shows the same thing.
This paper is far from being the last word on the topic of richness measurement. However, the stakes here are still high. If it is really true that new methods need to be used, then large bodies of literature in paleobiology and ecology could be called into question. Not only have scores of paleontological and ecological studies used quorum subsampling, but thousands of ecological studies have relied upon on species richness estimators in general. Thus, it is imperative that we start to evaluate such methods in a consistent and clear way using empirical data.
Data
As I will show, the differences among methods are sufficiently large that the choice of a dataset is not really important. However, I favor using ecological data for present purposes, because they escape the vagaries of the fossil record and of literature databases. Therefore, all the analyses involve tree and bat inventories that were drawn from the Ecological Register, a relational online database that has been described previously (Alroy Reference Alroy2015, Reference Alroy2017, Reference Alroy2018). Ten tree samples from Queensland derived from Bradford et al. (Reference Bradford, Murphy, Ford, Hogan and Metcalfe2014). The download was carried out on 25 September 2019, and the files are reposited on the Ecological Register site (http://ecoregister.org/?page=data). The tree and bat datasets were selected for use because they are large; the underlying surveys were conducted using consistent field methods; and, thanks to being global, the datasets capture a substantial amount of variation in raw species richness. For trees, 201 surveys including 12,326 register records were based on counts of individuals with diameters at breast height of about 10 cm. For bats, there were 220 surveys and 3386 records, and the counts were of mist-netted individuals. Default download criteria were otherwise used, except that replicate samples from the same 111 × 111 km grid cells were lumped instead of being discarded.
Note that in paleobiological taxonomic diversity studies, individuals are often replaced with counts of taxonomic occurrences (each occurrence is a presence in a sample) and inventories are then replaced with frequencies of occurrences over all samples within a geological time interval. In other words, paleobiological diversity curve studies concern distributions of taxon occurrences within intervals, whereas ecological studies concern distributions of individuals within samples. A few examples of the paleobiological approach, out of many, include Alroy (Reference Alroy2010a,b) and Starrfelt and Liow (Reference Starrfelt and Liow2016). This difference is scientifically substantive, but really should not be methodologically important: if a method can estimate richness from any frequency distribution, it should not matter whether the data are resolved to the level of individuals.
Methods
Shareholder Quorum Subsampling
The basic idea of quorum subsampling is to estimate the number of species that would be found in a subsample of an inventory if that subsample had a fixed amount of frequency distribution coverage, that is, met a quorum (Alroy Reference Alroy2010a,Reference Alroyb). Coverage is estimated as one minus the number of species sampled once (= singletons) divided by the number of individuals (Good Reference Good1953). A sample size correction to this equation has been proposed (Chao and Shen Reference Chao and Shen2010). Randomization must be used to handle complexities in the data when treating literature compilations of taxon presences (Alroy Reference Alroy2014). Otherwise, subsampled richness values can be computed analytically (Chao and Jost Reference Chao and Jost2012) by walking through all possible sample sizes in terms of the number of individuals and at each step computing expected richness using the method of Hurlbert (Reference Hurlbert1971: eq. 13) and expected coverage using an obvious and exact combinatorial expression that involves estimating the number of singletons, which is given by Chao and Jost (Reference Chao and Jost2012: eq. 4b). The returned richness value is the mean of the estimates found at the steps immediately below and above the requested quorum.
Note that Chao and Jost (Reference Chao and Jost2012) renamed the method “coverage-based rarefaction” on the grounds that equations and algorithms are different even if they seek the same goal. By contrast, when the original rarefaction algorithm of Sanders (Reference Sanders1968) was replaced by the equation of Hurlbert (Reference Hurlbert1971), there was universal agreement that the name would continue to be used.
I recently proposed a substantively different but related method called multiton subsampling (Alroy Reference Alroy2017). Because it produces estimates that are much more strongly correlated with quorum subsampling estimates than with extrapolations (Alroy Reference Alroy2017), I put it aside here in order to save space.
Squares
Squares (Alroy Reference Alroy2018) is structurally related to the Chao 1 equation (Chao Reference Chao1984) and Simpson's D (Simpson Reference Simpson1949), but aims to do something very different: provide reasonable richness estimates even when distributions are highly uneven. Let R′ = estimated true richness, n = the vector of species counts, N = the sum of n, S = the number of species in a subsample, and s 1 = the number of singletons. The extrapolator is:
$$R^{\prime} = S + s_1^2 \sum {{\bf n}^2} /\lpar N^2 + s_1S\rpar $$As mentioned, squares is expected to obey the replication principle (Hill Reference Hill1973). The reason is that, for example, doubling the counts doubles all the terms and sums of squared terms. Just as importantly, simple computations can show that putting aside random error, squares always yields lower bound estimates on the true value of R when underlying distributions are log normal, sampling is Poisson, and overall species pools are not extremely large (Alroy Reference Alroy2018). This is not true of some other common estimators such as the first- and second-order jackknifes (Burnham and Overton Reference Burnham and Overton1978; see also Colwell and Coddington Reference Colwell and Coddington1994). It has also already been shown that squares yields more accurate estimates than Chao 1 when distributions are at least somewhat uneven, with the gap between the two widening as evenness decreases (Alroy Reference Alroy2018). Squares does have a precision problem, but I will put Chao 1 aside in this paper for reasons of brevity and because the estimator introduced in the next section is just as accurate as Chao 1 but much more precise. The λ5 estimator (Alroy Reference Alroy2017) is also put aside, because it does not always yield lower-bound estimates and is demonstrably not as accurate as squares.
Corrected First-Order Jackknife
Chao 1 can be interpreted (Alroy Reference Alroy2017) as assuming that sampling is Poisson and that the rare species, including the singletons, doubletons, and those not sampled, have frequencies all governed by the same Poisson sampling rate λ. As suggested, this method has a major shortcoming: it is often highly imprecise because it relies heavily on two counts that may be small, that is, s 1 and the doubleton tally s 2. Chao 1 also handles uneven distributions poorly. To be specific, whenever a distribution of logged frequencies is normal, double geometric (Alroy Reference Alroy2015), or otherwise symmetrical, we expect the presence of very common species to imply the presence of very rare species. Chao 1 makes no use of this fact. Thus, its naive estimate of λ based on just s 1 and s 2 will be too high. Because of these issues, the method's extrapolations fall quite short when evenness is low (Alroy Reference Alroy2018).
The evenness problem is extremely difficult, and there may be no pure lower-bound estimator that is not only precise and has S as a minimum but outperforms Chao 1 by a large margin when it comes to this matter. However, the precision problem can be addressed by ignoring the s 2 count and focusing only on the s 1 count. A much older extrapolator called the first-order jackknife (Burnham and Overton Reference Burnham and Overton1978) does just that. The basic version of this equation is simply S + s 1. Alroy (Reference Alroy2017) noted that the counts s 1, s 2, etc. can be predicted by assuming that sampling follows a Poisson sampling process, in which case s 1 is easily estimated:
$$s_1 = R \; {\rm \lambda }\;{\rm e}^{-{\rm \lambda }} $$where, again, R = the true number of species in the sampling pool. Observed richness S is implied by a similarly straightforward equation:
$$S = R\lpar 1 - {\rm e}^{-{\rm \lambda }}\rpar.$$It follows that the jackknife is a simple function of R and λ:
$$S + s_1 = R\lpar 1 - {\rm e}^{-{\rm \lambda }} + {\rm \lambda } \,{\rm e}^{-{\rm \lambda }} \rpar.$$We now have a simple estimator of richness that generalizes the jackknife by allowing for variation in the sampling rate:
$$R' = \lpar {S + s_1} \rpar /\lpar 1 - {\rm e}^{-{\rm \lambda }} + {\rm \lambda }\,{\rm e}^{-{\rm \lambda }} \rpar.$$This equation has an interesting corollary: insofar as sampling is Poisson, Burnham and Overton's estimator must assume that λ = 1, because the terms in the denominator cancel out to 1 only when this is true. This assumption is implausible, because it literally means that the average species in the sampling pool is represented by one individual. The fact that their jackknife makes this flawed assumption explains why S + s 1 falls short when sampling is poor and goes over R when sampling is intermediate, which is easily shown by the kind of simulation illustrated by Alroy (Reference Alroy2018: Fig. S7).
All of this begs the question of how we can estimate λ if we are not to assume it has a fixed value like 1. Poisson theory implies several simple estimators of λ that involve combinations of S, s 1, and N, the total number of observations (Alroy Reference Alroy2017). A particularly robust version is ln(N/s 1), which has been called λ3 (Alroy Reference Alroy2017). It is especially noteworthy because it is a transform of Good's u, which as mentioned is 1 – s 1/N. Thus, λ3 equals ln[1/(1 – u)] and 1 – e−λ3 = u. Note the appearance of 1 – e−λ in the important equation (3), which suggests that u by itself can be used to extrapolate richness. In any event, I use the λ3 equation in this paper's analyses because it has good precision and yields realistically low values, except when sampling is very strong, in which case equation (5), which I call the corrected first-order jackknife or cJ1, will limit on R one way or another.
Now, equation (3) also can be inverted to produce a simple estimator of R, but one that puts aside the singleton count. This approach fails badly whenever distributions are uneven, because singletons are often still present even when sampling is heavy. Including the singleton count prevents this problem and yields a lower-bound estimator that matches Chao 1 in accuracy. Note, however, that when λ is high, the λ-based terms in the denominator of equation (4) limit on zero and its overall sum limits on 1, so we end back up with the original first-order jackknife equation.
The same simple theory behind cJ1 can be used to develop a large family of estimators, all of which meet the same general description: estimated R is some function of a numerator comprised of S, s 1, and perhaps other counts like s 2 divided by expressions involving λ that predict this numerator when multiplied by R. To save space, I will not elaborate on this family in more detail here, but I will note that some of these equations are less vulnerable to the evenness bias than Chao 1 or cJ1. The following is a good example of a precise estimator that indeed yields higher estimates when distributions are uneven:
$$\eqalign{R' & = [ {S^2/\lpar {S - s_1} \rpar } \rsqb /[ \lpar 1 - {\rm e}^{-\rm \lambda } \rpar ^{2} /\lpar 1 - {\rm e}^{-{\rm \lambda }} \cr & - {\rm \lambda }\,{\rm e}^{-{\rm \lambda }} \rpar \rsqb}.$$Because S – s 1 is the multiton count m (Alroy Reference Alroy2017), a reasonable shorthand for the equation is S 2/m. Note that this estimator and cJ1 have no particular connection to Chao 1 other than that all such methods can be derived using very basic Poisson theory linked with the assumption that underlying distributions are perfectly uniform, in which case they are always entirely accurate.
For the record, I know of only two other equations that consistently produce richness estimates that have little or no systematic bias when distributions are uniform. The first is 1 divided by the sample size–corrected version of Simpson's concentration index D (Simpson Reference Simpson1949; Hurlbert Reference Hurlbert1971), which has a very strong evenness signal and is therefore not usable as a richness estimator. The second is S/u, which as implied by the previous discussion equals the estimator R = S/(1 – e−λ) (a rearrangement of eq. 3) when λ3 is used to fix λ. This expression also suffers from variation in evenness, which explains why I have developed the S 2/m equation. It is of interest to note that S/u can be used to estimate overall richness from any subsampled estimate derived by quorum subsampling. The reason is that the method starts with a user-specified target value of u and then produces a subsampled richness estimate that can be substituted for S. Again, however, the evenness bias renders such rescalings untrustworthy.
cJ1 is easy to implement in R:
cJ1 and S 2/m have two advantages over squares. First, they have very good precision, as will be illustrated later for cJ1. Second, because they use N only indirectly (to obtain the sampling rate λ), they are less likely to be fooled by an outlying high count. To be specific, they conservatively assume that the only indicator of missing richness is the s 1 count when samples are very uneven, instead of allowing estimates to explode in such circumstances.
TRiPS
TRiPS was implemented using R functions written by J. Starrfelt and reposited on Dryad (https://doi.org/10.5061/dryad.86922). The two key functions are estimatePoiss and estimatetrue (equivalent to eqs. 2.3 and 2.5 in Starrfelt and Liow [Reference Starrfelt and Liow2016]). The second one is the actual richness estimator, and it requires output from the first one, which computes a Poisson sampling rate called λt that applies to an entire geological interval (= sample) labeled t. It is worth stressing that the sampling and richness estimates are not calculated simultaneously: λt is computed first, and richness is then estimated by scaling up on the basis of this parameter.
In any case, λt is related to a per-interval sampling probability estimated on the basis of taxonomic duration frequency distributions (Foote and Raup Reference Foote and Raup1996) or by examining sampling patterns over three intervals (Alroy Reference Alroy2008). For example, if λt = 0.1 and there are 10 time steps within an interval, then the chance of not drawing a taxon within the interval is exp(−10 λt) = 0.3679 and the chance of drawing it is 1 − 0.3679.
Sampling probabilities such as this one could be employed to produce richness estimates just by dividing observed values through by the per-taxon sampling probabilities. Alroy (Reference Alroy2014: Fig. 8) employed such a correction and suggested that it smoothed out some residual sampling variation even after quorum subsampling. It is easily shown that equation (2.5) in Starrfelt and Liow (Reference Starrfelt and Liow2016) reduces to the same calculation: the maximum-likelihood estimate of true richness given any sampling probability p (such as pbinom,t of Starrfelt and Liow [Reference Starrfelt and Liow2016]) and given any observed richness value N is simply N/p. Thus, the only new feature of TRiPS is the maximum-likelihood calculation used to estimate the Poisson sampling rate λt.
To make that calculation possible, Starrfelt and Liow (Reference Starrfelt and Liow2016) assumed the distribution of occurrences to be generated at random from a uniform underlying distribution. However, this assumption would be violated by any data, because all taxa are not created equal: many are legitimately rare and a few are legitimately common (Buzas et al. Reference Buzas, Koch, Culver and Sohl1982). Indeed, one would expect on first principles that occurrence counts would follow the log-normal distribution, and there is good evidence that this intuition is valid (Wagner and Marcot Reference Wagner and Marcot2013). Furthermore, no ecologist would assume a uniform underlying abundance distribution: uneven options such as the log normal and log series are favored instead (May Reference May, Cody and Diamond1975), and the double geometric distribution fits most ecological data very well (Alroy Reference Alroy2015). As I will show, although methods like Chao 1 and cJ1 also assume uniformity and also involve Poisson sampling models, their focus on singleton counts prevents them from falling into the same traps as TRiPS.
Regardless, the main point of this paper is not to challenge novelty or assumptions. Instead, my focus here is on performance.
First Test: Ability to Extrapolate
On a very basic level, an extrapolation method should always extrapolate unless we have strong a priori reasons to believe that sampling is complete. Usually, that means having found no singletons in a sample (Good Reference Good1953; Chao Reference Chao1984). A large majority of ecological samples do include singletons, so any extrapolative analysis of ecological data really should extrapolate.
TRiPS has very unexpected properties in this respect. When applied to the tree and bat data, TRiPS estimates richness values that are almost exactly the same as the raw, observed species counts (Fig. 1A). They are higher only in a handful of cases and only when richness is very high. The only apparent reason is the fact that ecological data are truly not close to being uniform in almost any real-world case, which seemingly causes TRiPS to assume that sampling is thorough.
Comparisons between observed counts of species and extrapolated richness estimates generated by the TRiPS method of Starrfelt and Liow (Reference Starrfelt and Liow2016). TRiPS adds little information to the counts. Points are ecological samples of trees and bats. Bat data points are nearly invisible because they are almost entirely overlain by tree data points. Lines of unity are shown. A, Based on raw data with no subsampling. B, Based on data that have first been randomly subsampled without replacement (i.e., rarefied) down to a level of 100 individuals.
A completely different and equally surprising pattern is found when we first standardize the data by rarefying them down to a level of 100 individuals (i.e., drawing 100 individuals from the abundance distribution without replacement). This cutoff seems challenging but fair: the median sample meeting it includes either 562 individuals (trees) or 297 individuals (bats), meaning that well over half the individuals are thrown out during each rarefaction trial. Note that there is one trial per sample.
In this case, TRiPS does consistently extrapolate beyond the data, but only when richness is relatively high and always in exact proportion to raw richness (Fig. 1B). In other words, the data seem to indicate a remarkable fact: if you know the number of individuals and you know the number of species in the raw dataset, then you know exactly what TRiPS will do. For example, if you start with 100 individuals and 50 species, TRiPS will give you back 62 species every single time. The higher the count of individuals, the lower the extrapolation. So with 200 individuals you will get back 51 species, and with 300 or more you will get 50. Indeed, 50 is the answer regardless of whether the 50 species each include six individuals or 49 species are singletons and the other species includes 251 individuals. Therefore, the reason the relationship is linear in the raw data instead of curvilinear (Fig. 1A) is that most of these samples are so large TRiPS is unwilling to extrapolate.
Thus, although TRiPS considers sample sizes and raw richness values, it is entirely insensitive to the shapes of abundance distributions. Other methods also assume that distributions are uniform, at least when it comes to the rarer species: Chao's famous extrapolation indices (Chao Reference Chao1984; Colwell and Coddington Reference Colwell and Coddington1994) are a good example (Alroy Reference Alroy2017). However, these methods still exploit the fact that there is powerful information in the frequency distribution, especially with regard to the number of singletons and how that compares to the size of other frequency classes.
Second Test: Consistency of Estimates
The number of taxonomic occurrences in an interval, or of individuals in a census, should bear no relation to reconstructed richness: we should obtain the same extrapolated estimate on average regardless of how much data we happen to have on hand, although obviously the estimate will be less precise if the sample is small. Only a few ecological extrapolators such as Chao 1 (Chao Reference Chao1984) even come close to behaving in this way (Colwell and Coddington Reference Colwell and Coddington1994), and even Chao 1 underestimates when sampling is truly poor (Alroy Reference Alroy2018). Any estimator that is highly robust should also be at most weakly sensitive to sample size.
One very simple way to test for such consistency is to compare the results obtained using full datasets to those obtained using the rarefied data also considered in the last analysis. In addition, to make the test more stringent, the data are passed through a “taphonomic” filtering step in another set of trials. Here it is assumed that each species has its own intrinsic per-individual sampling probability, meaning it is systematically either easy or hard to recover. The probabilities are drawn from a uniform distribution varying between 0 and 1. Imposing such a filter necessarily lowers subsampled richness, because some species effectively disappear from the sampling pool. It also necessarily increases the random error in the estimates. However, it does not necessarily create a nonlinear relationship between raw and subsampled richness.
No matter how the data are treated, quorum subsampling passes the rarefaction test and TRiPS fails it (Fig. 2). More specifically, the quorum subsampling estimates are somewhat variable but mostly fall on or very close to the line of unity. Without the taphonomic filter, there is no large difference between pairs of estimates for trees, albeit a small one for bats (Fig. 2A). The taphonomic filter lowers the values, as expected, so the averages are now different (Fig. 2B). However, the relationships are still approximately linear.
Performance of quorum subsampling (Alroy Reference Alroy2009) and TRiPS (Starrfelt and Liow Reference Starrfelt and Liow2016) based on full samples (i.e., raw data) and subsamples of 100 individuals randomly drawn without replacement (i.e., rarefied data). Quorum subsampling estimates are less sensitive to sample size. Data obscured by a “taphonomic” filter (see text) are illustrated on the right-hand side. Lines of unity are shown. A, B, Quorum subsampling (SQS) estimates (quorum = 0.95). C, D, TRiPS estimates.
By contrast, the squares, cJ1, and TRiPS estimates based on rarefied data are consistently and highly significantly lower no matter what the data treatment (Figs. 2C,D, and 3). This problem is predictable, because the rarefied samples in this case are extremely small, and even Chao 1 would fail with such data. Furthermore, a method such as TRiPS that generally fails to extrapolate (Fig. 1) will always be biased downward after rarefaction. More surprisingly, the relationships are highly nonlinear. Thus, squares, TRiPS, and to a lesser extent cJ1 are particularly unreliable when richness is very high relative to the number of individuals. The extrapolation estimates also show greater variability in general than the quorum subsampling estimates (Fig. 2). For example, the Spearman's rank-order correlation coefficient ρ is 0.981 for the quorum subsampling estimates based on the tree data (Fig. 2A) but 0.943 for the TRiPS estimates (Fig. 2E).
Performance of squares (Alroy Reference Alroy2018) and cJ1 (see text) based on raw, rarefied, and filtered data. TRiPS (Fig. 2), squares, and cJ1 are all lower-bound estimators, so their estimates are lower for rarefied data. However, cJ1 rarely overshoots. See caption of Fig. 2 for details. Lines of unity are shown. A, B, Squares estimates. C, D, cJ1 estimates.
Note, however, that the cJ1 data are cleaner than the squares data (Fig. 3): there are fewer low outliers, and there are far fewer cases in which the estimates based on rarefied data are actually higher than those based on the raw data. The absence of these overshoots emphasizes the point that cJ1 is not merely a lower-bound estimator, but a very reliable lower-bound estimator. This property facilitates easy, unambiguous interpretation of the estimates.
An even more important point is that squares and cJ1 extrapolate significantly beyond the raw counts. This explains why there is more variability along both axes in the squares and cJ1 plots (Fig. 3) than in the quorum subsampling plots (Fig. 2A,B). To be specific, the squares and cJ1 estimates for tree richness are spread over two orders of magnitude, but the quorum subsampling estimates range only over one order of magnitude (compare the axis scalings in Fig. 2A,B vs. Fig. 3). The idea that tree diversity varies from the cold temperate zone to the tropics by merely one order of magnitude would come as a surprise to many field ecologists.
To cut to the chase, then, subsampling succeeds in yielding internally consistent estimates but fails because it compresses those estimates. Extrapolation methods are dangerous when sampling is very poor, but it is better to sometimes fall short than always mislead.
Note that the sample sizes for quorum subsampling are lower because the quorum (0.95) is not met by many samples, whereas TRiPS will generate a value regardless of coverage. Thus, TRiPS is susceptible to error introduced by the quirks of sampling when datasets are small.
As for the seemingly high value of 0.95, it was selected because pedestrian ecological samples do meet this hurdle thanks to the fact that counts of individuals are much more uneven than counts of occurrences. This difference explains why quorum levels in paleobiological studies are normally much lower. The exact choice of a quorum is not really important because relative estimates of richness across samples hold up even when very low values are selected (Alroy Reference Alroy2010a). This fact explains why rarefaction of the data beforehand has so little effect on the subsampling results (Fig. 2A,B).
Third Test: Combination
The next challenge involves a much more stringent protocol called a split-analyze-and-sum test. Alroy (Reference Alroy2010a) illustrated analyses of this kind in which a Phanerozoic marine animal dataset was split between brachiopods and all other taxa (Alroy Reference Alroy2010a: Fig. 4) or alternatively split between latitudinal belts (Alroy Reference Alroy2010a: Figs. 7, 8). Alroy (Reference Alroy2010b) presented a more ambitious analysis in which 15 richness curves were generated for separate taxonomic groups and then summed. In all cases, the summed curves closely matched a single curve based on the original data.
The point of such a test is to see whether a method consistently produces accurate estimates given its assumptions. If the summed and original data match, then the partitions presumably have each been fairly sampled. If not, then at least one of the datasets on hand has certainly not been fairly sampled. An analogy involving raw samples of fair data might help. Suppose two people shop at a small grocery store, each with one basket, and they then figure out the total number of different food items they have collectively purchased. Now suppose a third person shops with two baskets. If that person comes up with more or fewer kinds of items, then that person's shopping has either been more or less efficient (and sampling by one team or the other or both must have been unfair).
The seriousness of the problem is even easier to understand in the case of extrapolating from unfair data. Suppose one person shops consistently for commonplace staples like milk and bread, while the other, whose palate is more refined, shops for many different things regardless of whether they are common. Extrapolating from the two baskets and summing should yield the total number of different food items in the store. Extrapolating from a single large basket collected from across the store should yield the same number.
Here I have used a particularly harsh version of the split-analyze-and-sum approach, analogous to the second shopping analogy: the species in each ecological sample are ordered from most common to least, the most common half are put in one dataset, the remainder are put in another, each partition is analyzed, and the richness estimates for the “common” and “rare” partitions are summed. A lower subsampling level (quorum = 0.8) is used in this analysis, because the higher cutoff of the preceding one often cannot be met by the rare partitions.
Quorum subsampling produces mixed results with this setup, because it oversamples the rare partitions. As a result, the summed values are higher than the original dataset values (Fig. 4A). However, the offset is moderate (median differences on a log scale = 0.534 for trees and 0.454 for bats). Furthermore, the relationships between summed richness and original dataset richness are consistent in the sense of being nearly log-linear (the least-squares slope on a log scale = 0.970 for trees and 0.867 for bats). The rank-order correlations between the values are also strong, although they are weaker for bats, because the ranges of values are compressed (Spearman's ρ = 0.969 for trees and 0.869 for bats).
Performance of quorum subsampling, squares, and cJ1 when applied to complete samples (raw data) and to split-analyzed-and-summed samples (split data). Squares and cJ1 are less sensitive to splitting, although there are outliers related to extreme undersampling. In the split data, each sample is split between common and rare species, richness values are estimated separately for each partition, and the estimates are summed. Lines of unity are shown. A, Quorum subsampling (SQS) estimates (quorum = 0.8). B, Squares estimates. C, cJ1 estimates.
Squares and especially cJ1 perform much better (Fig. 4B,C). There is no large, consistent offset in the values, which cluster tightly around the lines of unity despite the presence of some obvious outliers: the residual standard errors around a regression line are small (squares: 0.269 for trees and 0.103 for bats; cJ1: 0.348 for trees and 0.133 for bats). Correlations are also stronger (squares: ρ = 0.987 for trees and 0.988 for bats; cJ1: ρ = 0.971 for trees and 0.981 for bats).
With respect to cJ1, the outliers generally involve cases in which the rare partition only includes singletons and doubletons. Whenever this is true, the method is mathematically required to generate a much higher summed estimate. Putting aside these few cases would result in a better fit. Indeed, when samples having only singletons and doubletons in the rare partition are removed, the correlations for squares and cJ1 rise (quorum subsampling: ρ = 0.962 for trees and 0.865 for bats; squares: ρ = 0.992 for trees and 0.991 for bats; cJ1: ρ = 0.980 for trees and 0.992 for bats).
TRiPS has a much bigger problem: it produces unusual trends that are mostly an artifact of its inability to extrapolate when sample sizes are large (Fig. 5). At first glance, some of the TRiPS data seem to be similar to the quorum subsampling and squares data, also evidencing modest amounts of scatter (Fig. 5A). However, TRiPS places many points right on the line (Fig. 5A)—and these points represent large samples with high coverage of frequency distributions. All of the TRiPS points above the line (Fig. 5B) are undersampled, having been thrown out in the quorum subsampling plot because they failed to meet the sampling target. Furthermore, visible arcs in the undersampled data (Fig. 5B) hint at TRiPS’ unusual property (Fig. 1) of returning a fixed value insensitive to the shape of a frequency distribution whenever it is not uniform. Thus, the entire pattern says nothing about the split-analyze-and-sum test by itself. All it says is that TRiPS fails to extrapolate when sampling is reasonable (Figs. 1A and 5A) and extrapolates in a predictable way when sampling is poor (Figs. 1B and 5B).
Performance of TRiPS when applied to raw data and split-analyzed-and-summed samples. TRiPS does nothing when sampling is strong and is biased when sampling is weak. See caption of Fig. 4 for details. A, TRiPS estimates pertaining to well-sampled data. B, TRiPS estimates pertaining to undersampled data, which are defined as those points omitted in Fig. 4 because of failure to meet the 0.8 quorum.
In sum, what you get with a TRiPS analysis can depend on how you choose to segregate your data, such as by taxon or by geographic region. If you split so far as to render certain subsets badly undersampled, TRiPS will extrapolate in an uninformative way. If you do not, TRiPS will not extrapolate. On this point, Starrfelt and Liow (Reference Starrfelt and Liow2016) presented three curves for major dinosaur groups in addition to an overall curve without showing whether summing the former reproduced the latter. The question might have been interesting to answer.
Fourth Test: Evenness
One final matter remains to be addressed: do any of these methods capture a signal of abundance distribution evenness, and if so, is that a good thing? Here I take the simplest possible approach: using multiple regression to see whether the various estimates are functions of observed richness, observed evenness, or both. My reasoning is that all of these methods are supposed to reflect “true” richness, not evenness, so they had better perform as billed. An additional relationship with evenness would not be justifiable given any scenario.
For example, suppose that in the full sampling universe, there actually was a correlation between true evenness and true richness. After all, it is entirely plausible that rich communities tend to be even. That would create an indirect correlation between evenness and sampled richness, because actually poor, uneven communities would appear to be even more poor, being harder to sample, and vice versa. But a really good method's estimates would not correlate directly with evenness, because they would correlate perfectly with true richness, causing the evenness signal to wash out.
Unfortunately, a host of different evenness estimators exist, and there is no consensus about which one to use, with theorists arguing endlessly about technicalities and ignoring empirical properties (e.g., Smith and Wilson Reference Smith and Wilson1996; Kvålseth Reference Kvålseth2015). Here I employ two indices that are shown by simulation to have good empirical properties, regardless of hypothetical niceties: the classic J = H/ln S, where H is the Shannon index of diversity (Pielou Reference Pielou1966), and a new, very simple summation I call the ratio of means estimator:
$$[ \lpar \ln {\bf n}\rpar ^2/S\rsqb /\lcub [ \lpar \ln {\bf n}\rpar /S\rsqb ^2 + \lpar \ln {\bf n}\rpar /S\rcub $$This index has a lower bound of about 0.409, reached when a distribution consists of an infinite number of doubletons. If the distribution consists of an infinite number of singletons followed by a single doubleton, it has a lower bound of ln 2 (as can be shown algebraically). If the distribution is completely flat and the sample size is infinitely large, the expected value of 1 is actually higher. So, the ratio indicates that when samples are small but still even, they are even more even than expected by chance. Unlike most other estimators, the ratio of means climbs without an upper bound as evenness decreases, albeit very slowly. If a distribution includes one very common species and a series of singletons, it increases at a constant of very nearly ln 10 every time the common species increases in abundance by a factor of 10. The ratio is designed to be correlated only weakly with sampled richness, which is not something that could be said of J.
A simple simulation shows that estimates for J at its low end are so noisy as to be useless (Fig. 6A), but the ratio of means is reasonably precise and has a reasonable signal of variation in underlying abundances over its entire range (Fig. 6B). J does show less scatter at high evenness levels (SDs < 1), but it is particularly imprecise at low evenness values (SDs > 2) that appear to be characteristic of many ecological communities. Both indices have a less than perfect evenness signal (Fig. 6). But as suggested, the vulnerability of J to the richness bias is substantial and that of the new ratio is not (Spearman's ρ = 0.366 for richness vs. J; ρ = −0.154 for richness vs. the ratio of means).
Monte Carlo simulation analysis contrasting two evenness measures. Pielou's J (Pielou Reference Pielou1966) is less precise than the ratio of means estimator (see text). Each point is an estimate for a sample. Underlying abundances in each sample are log-normally distributed. The standard deviation of the abundances ranges randomly between 0 and 3; richness varies randomly on a log scale with a geometric mean of 20 and standard deviation of 1, rounded up to the nearest integer; and the number of drawn individuals ranges randomly on a log scale with a mean of 100 and a standard deviation of 1, rounded up. There are 1000 trials. The x-axis is the underlying standard deviation. A, Pielou's J. B, Ratio of means. The y-axis is logged.
Regardless, the actual richness estimates yielded by the four methods relate to the two indices in much the same way (Figs. 7, 8). Quorum subsampling is quite substantially related to evenness (Fig. 7A,B), which can be quantified by regressing its logged values against each of the indices combined with logged raw richness S (Table 1): the beta coefficients representing the influence of evenness are consistently high when using quorum subsampling, but close to zero when using the three extrapolation methods. The fact that very nearly the same coefficients are seen regardless of the evenness measure bolsters the argument that quorum subsampling, like Shannon's H, is a joint evenness + richness measure. Meanwhile, the extrapolation indices work more like pure richness measures when applied to real data. In other words, evenness biases seen in simulations like those of Alroy (Reference Alroy2018: Fig. S7) may not be a real problem in practice.
Relationships between quorum subsampling, TRiPS, and evenness metrics. Quorum subsampling (SQS) has an evenness signal and TRiPS does not. See also Table 1. A, SQS vs. Pielou's J. B, SQS vs. the ratio of means (see text). C, TRiPS vs. Pielou's J. D, TRiPS vs. the ratio of means.
Relationships between squares, cJ1, and evenness metrics. Squares and cJ1 also lack evenness signals. See also Table 1. A, Squares vs. Pielou's J. B, Squares vs. the ratio of means. C, cJ1 vs. Pielou's J. D, cJ1 vs. the ratio of means.
Beta coefficients produced by multiple regressions of logged richness estimates on evenness plus logged raw richness. cJ1 = corrected first-order jackknife; J = Pielou's J (Pielou Reference Pielou1966); RM = ratio of means; S = log of raw species richness; SQS = shareholder quorum subsampling; TRiPS = true richness estimated using a Poisson sampling model.
Conclusions
No method is perfect, and the fact that all four methods discussed in this paper have shortcomings should not be misinterpreted to mean either that I have overlooked a vastly better estimator—I do not believe one exists—or that the situation is hopeless. Quite to the contrary, I believe that because of its elegance and good performance, cJ1 is preferable to the current alternatives when datasets are anything other than badly undersampled (like the poor samples illustrated by Figs. 2 and 3). It rests on an elegant proof of its perfect accuracy when distributions are uniform, it extrapolates aggressively (Fig. 3C,D), its performance in the split-analyze-and-sum test is quite good (Fig. 4C), it has no large empirical evenness bias (Fig. 8C,D, Table 1), and it is a high lower bound on richness that yields expected values matching those of the less precise Chao 1 whenever distributions are log normal and uneven, which is a fair premise (Alroy Reference Alroy2015). cJ1 belongs to a larger family of equations, and some of them are even more accurate if not necessarily more precise, including S 2/m (eq. 6). I have dwelt on cJ1 only because it modifies a popular index in the older literature. And as for undersampling, when the data are truly poor and the abundances are uneven, it is arguable that species richness is unknowable anyway.
cJ1 and squares both have another advantage over subsampling: because they yield absolute richness estimates while working with raw counts, they also can produce per-interval sampling probabilities that are simply the ratio of the latter to the former. Although assuming one probability for all taxa found in each interval is an oversimplification (Wagner and Marcot Reference Wagner and Marcot2013), it is better than nothing, and such probabilities have good uses (such as cleaning up turnover rates and scaling diversity curves: Alroy Reference Alroy2008, Reference Alroy2014). One potential application would be to aid the reconstruction of phylogenies. The idea of minimizing stratigraphic gaps implied by phylogenies has a long history (e.g., Fisher Reference Fisher1988; Wagner Reference Wagner1998), and the explicit use of sampling probabilities figures in modern likelihood and Bayesian models (e.g., Stadler Reference Stadler2010; Bapst Reference Bapst2013). Interest in such methods has exploded in recent years (Hunt and Slater Reference Hunt and Slater2016). However, even the latest theoretical work in this area, such as Stadler et al. (Reference Stadler, Gavryushkina, Warnock, Drummond and Heath2018), continues to assume that the raw data are stratigraphic ranges and that per-interval sampling probabilities do not need to be informed directly by data. Methods like cJ1 may help to broaden these approaches.
Given the emphasis on cJ1 and squares, do I mean to say that we should give up on interpolation methods such as quorum subsampling? Not really. This method of standardization has some good properties, such as yielding consistent estimates regardless of how the raw data are degraded (Fig. 2A,B), more or less adhering to the replication principle (Alroy Reference Alroy2010a,Reference Alroyb), and not usually failing a split-analyze-and-sum test unless the data are truly ill behaved (as in Fig. 4A). The problem is that the meaning of subsampled values is not completely clear given that they are still compressed (Figs. 2A,B, and 4A) and still have an evenness signal (Fig. 7A,B). This evenness signal can become extreme when quorums are very low. Rescaling subsamples by u, as discussed in the “Methods” section, does not solve this problem at all. Furthermore, the stated goal of subsampling is only to obtain relative richness estimates that are accurate, but accurate absolute estimates are now often attainable, and they are more intuitive. Nonetheless, I do believe that subsampling methods have a place at the table and that better ones can be developed.
TRiPS, meanwhile, is said to be “pregnant with possibilities,” because it sequentially estimates “both sampling rates and true richness directly” (Starrfelt and Liow Reference Starrfelt and Liow2016, emphasis in original). An editorial went further by asserting that it would “revolutionize the ways in which palaeontologists and evolutionists in general treat fossil record data” (Ezard et al. Reference Ezard, Quental and Benton2016). Given those high stakes, paleobiologists should only employ TRiPS—or any other method—if it can be validated firmly. Squares and cJ1 have largely met this challenge; TRiPS has not. And TRiPS substantively disagrees with the other methods (Fig. 9), yielding more compressed estimates regardless of the fact that those estimates are highly correlated with squares and cJ1 values (Fig. 9C). For example, a simple log-log regression (dashed line in Fig. 9C) predicts that, on average, a cJ1 estimate of 200 species would correspond to a TRiPS estimate of only 163 species. The major reason for this failure is that TRiPS simply declines to extrapolate whenever sampling is at least fair (e.g., Fig. 1). Thus, diversity curves produced by TRiPS might often be substantively different, but only for methodological reasons.
Relationship between quorum subsampling (SQS) (quorum = 0.95), squares, and cJ1 estimates on the one hand and TRiPS estimates on the other, all based on raw data. cJ1 is particularly close to TRiPS, but it decompresses the raw data, so its high estimates are higher. Solid lines are lines of unity; dashed lines show least-squares regression fits. A, Quorum subsampling vs. TRiPS. B, Squares vs. TRiPS. C, cJ1 vs. TRiPS.
Regardless of the methods under consideration, however, it is hoped that the empirical validation criteria laid out in this paper will prove to be helpful the next time researchers propose a major innovation in the area of richness inference. Meeting most or all of these criteria is necessary for validation, if not necessarily sufficient. It would not be circumspect to instead go forward only on the basis of simulations, much less without any particular justification. One reason is that any method can be shown to work by simulation merely by making sure its assumptions are met by the model. For example, Starrfelt and Liow (Reference Starrfelt and Liow2016) varied per-species sampling rates by less than a factor of two, consistent with their method's premise of virtually uniform sampling regimes. A more stringent test would have been provided by instead assuming that sampling rates track something like a log-normal distribution and varying the parameters of this distribution over a realistic range, as in Alroy (Reference Alroy2018).
Acknowledgments
I thank R. Benson, M. Clapham, and D. Nicholson for alerting me to TRiPS and some of its properties. R. Benson, G. Lloyd, D. Nicholson, J. Tennant, A. Tóth, P. Wagner, S. Wang, and an anonymous reviewer provided helpful comments. This is Ecological Register publication number 7.
