A Richness Estimator Based on Integrated Data

1. Introduction

Species richness is the most commonly used quantitative diversity metric and is easily understood. The term "species" can be broadly defined to include biological species, software bugs, words in a book, genes, alleles, or other discrete entities, as reviewed in [1,2,3]. This article focuses on biological applications, specifically the number of detectable species within a given area. Ecological studies use two primary data formats to assess species diversity: individual-based abundance data and sample-based incidence data. Individual-based abundance data involve randomly sampling and identifying individual organisms to species. Sample-based incidence data involve randomly sampling a plot, quadrat, trap, transect, or net from the target region and recording the incidence of species appearing in the sampled unit.

However, producing an inventory species list in a target area is impossible due to resource or sampling limitations. In practice, a random sample from the target community, i.e., a small proportion of area size or community size in the target area, is typically used to assess species diversity; as such, the observed richness in the sampled sample always underestimates true richness. To address the negative bias of observed richness, dozens of estimators have been proposed for these two types of data and have led to considerable progress in various disciplines[1,2,4,5]. In ecological studies, the most widely used estimators are Chao’s lower bound estimators [6,7] and Jackknife-based estimators[8,9]. These estimators were derived without making assumptions about species detection rates. They are classified as nonparametric approaches and, as such, provide more robust estimates. In addition, these nonparametric estimators do not require all the information on observed species; only rare species (singletons and doubletons) are used in the sample to estimate undetected richness. These estimators could show expected robust statistical behavior only when the sampling unit (i.e. a individual in abundance data and a plot in incidence data) is indpendently or randomly sampled. However, due to resource constraints, a random sample is often only feasible in a limited area of a local region and not in a large-scale region.

In recent decades, monitoring species richness to reveal the impact of human activities on a large or global scale is an increasingly urgent task [10,11,12,13]. However, richness estimation of a large-scale region (or multiple communities) is still a statistical challenge for which no reliable estimator has been developed until now. In general, the collected data sets used to estimate the richness of a large-scale region usually consist of many samples that are separately sampled from the subregions by implementing different sampling schemes. Therefore, this integrated data set is composed of different kinds of data formats including individual-based abundance data and sample-based incidence data. However, the widely-used rigorous estimators in the literature have their limitations due to their underlying theoretical assumptions and are not equipped to analyze this type of integrated data.

In this article, I provide a theoretical interpretation of the applicability of Chao's lower bound estimator for estimating species richness based on a pooled sample of integrated data. Additionally, utilizing the Good-Turing frequency formula[14], I address the negative bias inherent in Chao's lower bound estimator and propose a bias-corrected alternative. The variance of the new estimator can be calculated through the asymptotic approach and its 95% confidence interval can be obtained through logarithmic transformation. To evaluate the efficacy of this proposed estimator, three commonly used ecological models and two real datasets are utilized in simulation studies and illustrative examples.

2. Materials and Methods

2.1. Sampling distribution model

Assume there are a total of $S$distinct species in the assemblage of interest. In ecological studies, individual-based abundance data and sample-based incidence data are the most commonly collected data types to assess richness diversity [15]. The sampling unit of individual-based abundance data is an individual randomly sampled and identified to species, and the sampling unit of sample-based incidence data is a plot, quadrat, trap, transect or net randomly sampled from the target subregion and only the incidence (presnece or absence) of species appearing in the selected plot is recorded.

For individual-based abundance data, assume $n$(a small fraction of assemblage size) individuals are randomly sampled and identified to species from the target region by sampling with replacement or sampling without replacement. Let $X_i$ be the number of individuals of species i counted in the sample. The species frequency or species abundance $(X_1,X_2,\dots ,X_S)$ could be assumed to follow a multinomial distribution with size n and probabilities $(p_1,p_2,\dots ,p_S)$ and the species frequency $X_i$ follows a binomial distribution with parameters $n$and $p_i$,

$$X_i~Binomial\left(n,p_i\right), i=1, 2,\dots ,S,$$

where $p_i$ is the relative detection probability of species i. Let $f_k={\sum }_{i=1}^{S}I(X_i=k)$ be the number of species that are observed exactly k times in the sample, $k=1, 2,\dots ,n$. Therefore, $f_0, f_1, f_2, f_3$ are respectively the unseen richness, singleton richness, doubleton richness and tripleton richness in abundance sample.

For sample-based incidence data, assume $t$ sampling units are randomly sampled from the target region and only the incidence (presence or absence) of species in the sampled unit is recorded. Let $Y_i$ be the number of units in which species i is detected in the t sampled units. Then, $Y_i$ could be assumed to follow a binomial distribution with size t and probability ${\pi }_i, i=1, 2, \dots ,S$,

$$Y_i~Binomial\left(t,{\pi }_i\right), i=1, 2,\dots ,S,$$

where ${\pi }_i$is the detection probability of species i, which depends on the abundance, body size and color of species i as well as the investigator’s capability. Let $Q_k={\sum }_{i=1}^{S}I(Y_i=k)$ be the number of species that are detected in exactly k out of t sampling units, $k=1, 2,\dots t.$ Therefore, $Q_0, Q_1, Q_2, Q_3$ are respectively the unseen richness, singleton richness, doubleton richness and tripleton richness in the incidence sample.

2.2. Richness estimation for one assemblage based on integrated data

Assuming that N samples are randomly collected from the assemblage through various sampling schemes, including individual-unit-based and sample-unit-based sampling methods, the integrated data comprises two formats (i.e., abundance data and incidence data), as commonly seen in ecological studies. To determine the richness of the assemblage, Chao1 and Chao2[6,7] derived without model assumption on species detection rates are the most commonly used estimators, which are briefly outlined below.

2.2.1. Chao’s lower bound estimators

Based on the Cauchy-Schwarz inequality and without making any assumptions on species detection rates, Chao proposed lower bound estimators for richness in 1984 and 1987. These estimators were designed for individual-based abundance data and sample-based incidence data, and are referred to as the Chao1 and Chao2 estimators, respectively. The Chao1 and Chao2 estimators are separately expressed as

$$Chao1=S_{obs}+\left\{\begin{array}{cc}\frac{f_1^2}{2f_2}& if\hspace{1em}{f}_2>0\\ \frac{f_1^{}(f_1-1)}{2}& if\hspace{1em}{f}_2=0\end{array}\right.,$$

(1)

$$Chao2=S_{obs}+\left\{\begin{array}{cc}\frac{t-1}{t}\frac{Q_1^2}{2Q_2}& if\hspace{1em}{Q}_2>0\\ \frac{t-1}{t}\frac{Q_1^{}(Q_1-1)}{2}& if\hspace{1em}{Q}_2=0\end{array}\right.,$$

(2)

Chao’s lower bound estimators only use the frequency counts of the two rarest species (i.e. the numbers of singleton and doubleton species) in the sample to estimate undetected richness.

On the basis of the Cauchy-Schwarz inequality theory, Chao’s lower bound estimators are unbiased when the detection rates of species are homogeneous (i.e. $p_i=\frac{1}{S}$ in Eq.1 or ${\pi }_i=c$ in Eq.2, for $i=1, 2, \dots ,S$). In addition, according to the concept of the Good-Turing frequency formula[14], Chao et al.[16,17] shew that Chao’s lower bound estimators are nearly unbiased estimators only when rare species have approximately homogenous detection probabilities (or rates). Therefore, the degree of heterogeneity of abundant species in the assemblage contains no information about the unbiasedness of Chao’s lower bound estimators. When the detection rate of rare species is highly heterogeneous or the sample size is not large enough, in contrast to other parametric estimators, Chao1 or Chao2 can provide a lower bound and robust richness estimate [2,18]. However, Chao1 and Chao2 were separately derived based on different sampling models for abundance data and incidence data. Importantly, it still has no theoretical evidence or proof that the existing estimators can be used to estimate species richness using a pooled sample of integrated data.

2.2.2. Extension of Chao’s lower bound estimators for integrated data

Many richness estimators proposed in the literature, whether they are parametric or non-parametric, are designed for randomly sampled data. This means that the detection rate of a species for each random trial, such as a selected individual or plot, is assumed to be identical. These estimators assume that the underlying assumptions of the binomial distribution are met

However, if N samples are collected using different sampling methods (e.g. sampling schemes, sampling efforts, plot sizes, or investigators) from the target region, the observed species count in the pooled sample no longer follows a binomial distribution. This violates the theoretical assumption of a random sample. This type of integrated data is often encountered in ecological studies where individual-based abundance data and sample-based incidence data are collected from the same target region. Although integrated data is widely used to estimate the richness of a large-scale or global region, no richness estimator has been rigorously designed to handle integrated data. In this section, I demonstrate theoretically that Chao's lower bound estimator can be modified to handle integrated data.

For individual-based abundance data, according to probability theory, when sample size$n$ is sufficiently large and relative abundance (or detection probability) $p_{ }$is sufficiently small, the species frequency ($X_{ }$) follows a binomial distribution that converges to a Poisson distribution. It implies that the frequency ($X$) of rare species (i.e.$p_{ }$is sufficiently small) in the sample could approximate a Poisson distribution with mean $n_{ }{p}_{ }$( i.e. $X_{ }~ poi\left(n{p}_{ }\right)$) for species with low detection rate. This convergence feature also applies to sample-based incidence data. When the number of plots $t$ is large and the detection rate ${\pi }_{ }$ tends to zero, the incidence count ($Y_{ }$) of rare species in the sample could approximate a poison distribution with mean $t_{ }{\pi }_{ }$( i.e. $Y_{ }~ poi\left(t\pi \right)$) for species with low detection rate.

Without loss of generality, two random samples are collected from the target region through different sampling schemes, namely individual-based sampling and plot-based sampling methods. These samples correspond to individual-based abundance data and sample-based incidence data, respectively. When the two sampled samples are pooled, the pooled species frequency $Z_i=X_i+Y_i$ represents the count of species i in the pooled sample. Here, $Z_i$ is no longer a random variable following a binomial distribution.

Based on the convergence principle between the binomial and Poisson distributions discussed earlier, we have, for species with low detection rates in the combined sampling scheme (i.e. small $p_i$ and small ${\pi }_i$), the abundance ($Z_i$) in the pooled sample approximately follows a Poisson distribution with mean parameter ${\lambda }_i=n{p}_i+t{\pi }_i$. For simplicity, let denote this mean parameter as ${{\lambda }_i=md}_i$, where $m$represents the unknown size of the pooled sample and $d_i$ represents the detection rate of species i. Next, let $G_k={\sum }_{i=1}^{S}I\left(Z_i=k\right)$ be the species frequency count, representing the number of species that are present exactly k times in the pooled sample. When k is small (e.g. k = 0, 1, 2 or 3) and the size of the pooled sample is sufficiently large, $G_k$ is primarily contributed by the rare species, which approximately follow a Poisson distribution. Given a specific sampling scheme, all species in the region can be divided into a set of rare species denoted as $\left\{s_{rare}\right\}$ and a set of abundant species denoted as $\left\{s_{abun}\right\}$. Based on the existing convergence theory between the binomial distribution and the Poisson distribution, we have the approximation of the expectation of $G_k$ for small k,

$E\left[G_k\right]=E\left[{\sum }_{i=1}^{S}I\left(Z_i=k\right)\right]={\sum }_{i\in \left\{S_{rare}\right\}}P\left(Z_i=k\right)+{\sum }_{i\in \left\{S_{abun}\right\}}P\left(Z_i=k\right)$.

Since when k is small, the probability that abundant species have a count of k tends to zero. Therefore, ${\sum }_{i\in \left\{S_{abun}\right\}}P\left(Z_i=k\right)$ is roughly equal to 0. We have

$$E\left[G_k\right]\approx {\sum }_{i\in \left\{S_{rare}\right\}}P\left(Z_i=k\right)+0.$$

According to the convergence property between binomial and Poisson distribution for rare species, the following approximation is held:

$$E\left[G_k\right]\approx {\sum }_{i\in \left\{S_{rare}\right\}}P\left(Z_i=k\right)\approx {\sum }_{i\in \left\{S_{rare}\right\}}\frac{{\lambda }_i^k}{k!}{e}^{-{\lambda }_i}\approx {\sum }_{i=1}^S\frac{{\lambda }_i^k}{k!}{e}^{-{\lambda }_i}.$$

Therefore, we can derive the following four approximation equations for the expectation of undetected richness, singleton richness, doubleton richness, and tripleton richness, which represent the number of rare species in the pooled sample.

$$E\left[G_0\right]=E\left[{\sum }_{i=1}^{S}I\left(Z_i=0\right)\right]={\sum }_{i=1}^{S}P\left(Z_i=0\right)\approx {\sum }_{i=1}^S{e}^{-{\lambda }_i}$$

(3.1)

$$E\left[G_1\right]=E\left[{\sum }_{i=1}^{S}I\left(Z_i=1\right)\right]={\sum }_{i=1}^{S}P\left(Z_i=1\right)\approx {\sum }_{i=1}^S{{{\lambda }_i}^{}e}^{-{\lambda }_i}$$

(3.2)

$$E\left[G_2\right]=E\left[{\sum }_{i=1}^{S}I\left(Z_i=2\right)\right]={\sum }_{i=1}^{S}P\left(Z_i=2\right)\approx {\sum }_{i=1}^S\frac{{{\lambda }_i}^2}{2}{e}^{-{\lambda }_i}$$

(3.3)

$$E\left[G_3\right]=E\left[{\sum }_{i=1}^{S}I\left(Z_i=3\right)\right]={\sum }_{i=1}^{S}P\left(Z_i=3\right)\approx {\sum }_{i=1}^S{\frac{{{\lambda }_i}^3}{6}e}^{-{\lambda }_i}$$

(3.4)

It is worth emphasizing once again that these approximations are valid only for lower species frequency counts in the sample, under the condition that sample sizes ($n$and $t$) are sufficiently large. In Appendix A, I provide evidence that the aforementioned approximate equations hold by demonstrating their validity through numerical simulations.

Based on the Cauchy-Schwarz inequality, we have the inequality

$${\sum }_{i=1}^S{e}^{-{\lambda }_i}{\sum }_{i=1}^S{{\lambda }_i}^2{e}^{-{\lambda }_i}\ge {\left({\sum }_{i=1}^S{{\lambda }_i}^{}{e}^{-{\lambda }_i}\right)}^2,$$

(4)

According to Eqs. 3.1-3.2, Eq. 4 is equivalent to $E\left[G_0\right]E\left[2G_2\right]\ge E{\left[G_1\right]}^2$. This inequality is also held when species detected mean abundance ${\lambda }_i$ is assumed to be a random variable with probability density function $f\left(\lambda \right)$, expressed as

$$\left(S\int e^{-{\lambda }_{}}f\left(\lambda \right)d\lambda \right)\left(S\int {\lambda }^2{e}^{-{\lambda }_{}}f\left(\lambda \right)d\lambda \right)\ge {\left(S\int {\lambda }^{}{e}^{-{\lambda }_{}}f\left(\lambda \right)d\lambda \right)}^2.$$

Therefore, we have the lower bound estimator of undetected richness ${\widehat{G}}_0=\frac{G_1^2}{2G_2},$ which uses the number of singletons and doubletons in the pooled sample to estimate undetected richness. Therefore, the proposed richness estimator could be interpreted as an extension of Chao1 or Chao2. It is denoted as Chao3, and the modified formula can be expressed as

$$Chao3=S_{obs}+\left\{\begin{array}{cc}\frac{G_1^2}{2G_2}& if\hspace{1em}{G}_2>0\\ \frac{G_1\left(G_1-1\right)}{2}& if\hspace{1em}{G}_2=0\end{array}\right.$$

(5)

2.2.3. Adjusted Lower bound estimator of species richness for integrated data

According to the Cauchy-Schwarz inequality, we know Chao3 is a lower bound estimator. Based on the concept of the Good-Turing frequency formula, Chao3 will be severely negatively biased when the rare species have a high degree of heterogeneity. In this section, the negative bias of Chao3 can be corrected based on the concept of Good-Turing frequency formula[14]. The bias of Chao3 is approximately equal to

$$E\left[\frac{G_1^2}{2G_2}\right]-E\left[G_0\right]\approx \frac{{\left({\sum }_{i=1}^S{m{d}_i}^{}{e}^{-m{d}_i}\right)}^2}{2{\sum }_{i=1}^S\frac{{\left(m{d}_i\right)}^2}{2}{e}^{-m{d}_i}}-{\sum }_{i=1}^S{e}^{-m{d}_i}.$$

Let $d_{\left(r\right)}={\sum }_{i=1}^S{d}_{i}I(Z_i=r)/G_r$be the mean detection rate of species which are present r times in the pooled sample; it can be estimated using the modified Good-Turing frequency formula as

$${\widehat{d}}_{\left(r\right)}=\frac{\left(r+1\right)G_{r+1}}{{mG}_r}.$$

Using the modified Good-Turing frequency formula, we have the following approximate equations

$$E\left[G_1\right]\approx {\sum }_{i=1}^S{m{d}_i}^{}{e}^{-m{d}_i}={\sum }_{i=1}^S\frac{2}{m{d}_i}E\left[I\left(Z_i=2\right)\right]\approx \frac{2}{{md}_{\left(2\right)}}E\left[G_2\right]$$

(6.1)

$$E\left[G_0\right]\approx {\sum }_{i=1}^S{e^{-m{d}_i}}^{}={\sum }_{i=1}^S\frac{1}{m{d}_i}E\left[I\left(Z_i=1\right)\right]\approx \frac{1}{{md}_{\left(1\right)}}E\left[G_1\right]$$

(6.2)

According to Eqs. 6.1-6.2, the bias of $\frac{G_1^2}{2G_2}$ can be approximately derived as

$$Bia{s}_{Chao3}=E\left[\frac{G_1^2}{2G_2}\right]-E\left[G_0\right]\approx (\frac{1}{{md}_{\left(2\right)}}-\frac{1}{{md}_{\left(1\right)}})E\left[G_1\right]$$

Therefore, the bias of Chao3 can be estimated by replacing $d_{\left(1\right)}$ and $d_{\left(2\right)}$with ${\widehat{d}}_{\left(1\right)}$ and ${\widehat{d}}_{\left(2\right)}$, respectively. It is given as

$${\widehat{Bias}}_{Chao3}=\frac{G_1^2}{2G_2}(\frac{2G_2^2}{3G_1{G}_3}-1)$$

Then, we have the bias-corrected estimator of Chao3, expressed as

$$Chao{3}_{Adj}=S_{obs}+{\frac{G_1^2}{2G_2}\left(2-\frac{2G_2^2}{3G_1{G}_3}\right)}^{-},$$

(7)

where ${\left(A\right)}^{-}$ equals $A$ if $A<1$ and 1 if $A\ge 1$. Here, $G_3$ (or $G_1)$ is replaced by 1 as $G_3=0$ (or $G_1=0$) to make Eq. 7 always well-defined. The mathematic form of the estimator shown in Eq. 7 is identical to the parametric estimator proposed by Chiu[18, 19] which was derived based on the Beta-Binomial mixture model for sample-based incidence data or based on the Gamma-Poisson mixture model for individual-based abundance data. The new estimator can also be proved to be a lower bound of richness under the incidence-based Beta-Binomial mixture model or the abundance-based Gamma-Poisson mixture model [19,20].

Since the Chao3_Adj estimator utilizes the first three rarest species in the sample to estimate undetected richness, it can be applied to integrated data consisting of multiple samples randomly collected from the target region without adhering to a specific sampling model or scheme. In order to aid readers, a table has been created in Appendix B, which outlines the equations and symbols used in the proposed estimators, along with their meanings, origins (with references), as well as their advantages and disadvantages.

2.3. Regional richness estimation based on integrated data

When the target region is divided into N subregions, species sampling data are collected independently and separately from each subregion. The sampling data can be collected by either an individual-unit-based sampling method or a sample-unit-based sampling method. Let $X_{ij}$ represent the number of sampling units (such as the number of individuals in individual-based abundance data or the number of plots in sample-based incidence data) of species $i$ in the sample $j$, which is collected from the $jth$ subregion. Here, i ranges from 1 to $S$ for the $S$ species, and j ranges from 1 to $N$ for the $N$ subregions$.$If the sample size (i.e. the total number of individuals in abundance data or the total number of plots in incidence data) is sufficiently large in each sample, the counts ($X_{ij}$) of species with low detection rates will approximate a Poisson distribution. Then, the total count of species i in the pooled sample, denoetd as $X_{i+}={\sum }_{j=1}^N{X}_{ij}$, will approximate a Poisson distribution when the detection rate of species i in each subregion is uniformly low.

Letting $G_k={\sum }_{i=1}^{S}I\left(X_{i+}=k\right)$ be the number of species with a count of exactly k in the pooled sample. The approximate equations shown in Eqs. 3.1-3.4 are also applicable to the pooled sample of integrated data. Additionally, formulas for Chao3 and Chao3_Adj can be derived to estimate regional species richness based on the Good-Turing frequency formula, without making any specific model assumptions. Similarly, their variance estimators can be obtained using the asymptotic approach, and the 95% confidence interval (CI) of species richness can be derived by referring to the discussion surrounding Eq. 9.

According to the derivation, the proposed richness estimator possesses the following properties: (i) When the samples are individually and randomly collected from each subregion, the sampled samples can be directly combined to estimate undetected richness, regardless of whether the data format in each sample is identical or not. (ii) The estimation of undetected richness is solely based on the frequency counts of the rarest species in the pooled sample. (iii) When the detection rates of rare species are homogeneous (including the homogeneous model as a special case) or the sample size is sufficiently large, the proposed estimators are nearly unbiased.

2.4. Variance estimation

To derive the variance estimators for the discussed richness estimators, an asymptotic approach is employed. By defining $G_{2+}$ as the total frequency count of species with a count of at least 3 in the sample (i.e., $G_{2+}={\sum }_{k\ge 3}^{ }{G}_k$), the estimator ${\widehat{S}}_{Chao3}$ can be expressed as a function of $(G_1,G_2, G_{2+})$. The estimator of ${\widehat{S}}_{Chao3}$’s variance could be obtained by the asymptotic approach in which $(G_0,G_1,G_2,G_{2+})$ approximate a multinomial distribution with parameters $(S,\frac{E\left[G_0\right]}{S},\frac{E\left[G_1\right]}{S},\frac{E\left[G_2\right]}{S},\frac{E\left[G_{2+}\right]}{S})$. Additionally, let$G_{4+}$be the total frequency count of species with a count of at least 4 in the sample (i.e., $G_{4+}={\sum }_{k\ge 4}^{ }{G}_k$), then Chao3_Adj becomes a function of $(G_1,G_2, G_3,G_{4+})$. The estimator of $S_{Adj.Chao3}$'s variance can also be obtained using the asymptotic approach, where$(G_0,G_1,G_2,G_3,G_{4+})$ approximately follow a multinomial distribution with parameters $(S,\frac{E\left[G_0\right]}{S},\frac{E\left[G_1\right]}{S},\frac{E\left[G_2\right]}{S},\frac{E\left[G_3\right]}{S},\frac{E\left[G_{4+}\right]}{S})$.

The variance estimator of the Chao3 or Chao3_Adj can be derived by the delta method and is expressed as

$${\widehat{var}({\widehat{S}}_{}}_{})\approx \sum _i^{}\sum _j^{}\frac{\partial {{\widehat{S}}_{}}_{}}{\partial G_i}\frac{\partial {{{\widehat{S}}_{}}_{}}_{}}{\partial G_j}\widehat{cov}\left(G_i,G_j\right),$$

(8)

where $\widehat{cov}\left(G,G_j\right)=\left\{\begin{array}{cc}{G}_i(1-G_i/{{\widehat{S}}_{}}_{})& ifi=j\\ -G_i{G}_j/{\widehat{S}}_{}& ifi\ne j\end{array}\right.$.

To derive the 95% confidence interval (CI) of species richness and to ensure the lower bound of the 95% CI of species richness is larger than the observed richness, assume ${\widehat{S}}_{ }-S_{obs}$ follows a log-normal distribution [18,21]. Then the two-sided 95% CI of species richness is obtained as

$$[S_{obs}+\frac{{\widehat{S}}_{}-S_{obs}}{R},S_{obs}+(\widehat{S}-S_{obs}\left)R\right],$$

(9)

where$R=exp\left\{1.96{\left[\log \left(1+\frac{Var\left({\widehat{S}}_{}\right)}{{\left({\widehat{S}}_{}-S_{obs}\right)}^2}\right)\right]}^{\frac{1}{2}}\right\}$.

When samples are randomly collected, Chao3 consistently provides a lower bound estimate of species richness. Similarly, the Chao3_Adj also provides a lower bound estimate of species richness under the Gamma-Poisson model or the Beta-Binomial model [19,20]. Therefore, in cases where the community exhibits high heterogeneity or the sample size is small, these two estimators can offer lower bound estimates and more informative one-sided 95% confidence intervals (CIs) of species richness, shown as

$$[S_{obs}+\frac{{\widehat{S}}_{}-S_{obs}}{R},\mathrm{\infty }],$$

where $R=exp\left\{1.65{\left[\log \left(1+\frac{Var\left({\widehat{S}}_{}\right)}{{\left({\widehat{S}}_{}-S_{obs}\right)}^2}\right)\right]}^{\frac{1}{2}}\right\}$ .

3. Results

3.1. Simulation study

A simulation study was conducted to examine the statistical behaviors of the new estimators. The study involved the use of three species abundance models to generate individual-based abundance data and three species detection models to generate sample-based incidence data. The number of species was kept constant at S = 600, and the simulated datasets were generated separately and independently using the following models.

3.1.1. For individual-based abundance sampling models:

The species detection probabilities (or species relative abundance) $\left(p_1,p_2,\dots ,p_S\right)=({ca}_1,{ca}_2,\dots ,{ca}_S)$ in each model are provided below, where c is a normalizing constant such that ${\sum }_{i=1}^S{p}_i=1$. The coefficient of variation (CV) of $(p_1,p_2,\dots ,p_S)$ is also presented to indicate the degree of heterogeneity of $(p_1,p_2,\dots ,p_S)$.

a.

Abundance model 1, random uniform model (CV = 0.53): with $p_i=c{a}_i,i=1,2,\dots ,S$, where $\left(a_1,a_2,\dots ,a_S\right)$is a random sample from a uniform distribution.

b.

Abundance model 2, broken-stick model (CV = 0.97): with $p_i=c{a}_i, i=1, 2,\dots ,S$, where $\left(a_1,a_2,\dots ,a_S\right)$ is a random sample from an exponential distribution with parameter 1. This model is commonly used in the literature and is equivalent to the Dirichlet distribution.

c.

Abundance model 3, log-normal model (CV = 1.56): with $p_i=c{a}_i, i=1, 2,\dots ,S$, where $\left(a_1,a_2,\dots ,a_S\right)$ is a random sample from a log-normal distribution with parameters 0 and 1.

3.1.2. For sample-based incidence sampling models

The species detection probabilities $\left({\pi }_1,{\pi }_2,\dots ,{\pi }_S\right)=({ca}_1,{ca}_2,\dots ,{ca}_S)$ in each model were determined, where c is a rescaling constant such that the maximum detection probability is a fixed at a constant value. The coefficient of variation (CV) of $({\pi }_1,{\pi }_2,\dots ,{\pi }_S)$ is also calculated to indicate the degree of heterogeneity of $\left({\pi }_1,{\pi }_2,\dots ,{\pi }_S\right).$

d

Incidence model 1: the random uniform model (CV = 0.57): where ${\pi }_i=c{a}_i, i= 1, 2,\dots ,S$, and $\left(a_1,a_2,\dots ,a_S\right)$is a random sample from a uniform distribution

with parameters (0, 1), and scale c is used to control the maximum of ${\pi }_i$.

e

Incidence model 2: the broken stick model (CV = 0.99): where ${\pi }_i=c{a}_i, i= 1, 2,\dots ,S$, and $\left(a_1,a_2,\dots ,a_S\right)$ is a random sample from an exponential distribution with parameter 1, and scale c is used to control the maximum of ${\pi }_i$. This model is commonly used in the literature and equivalent to the Dirichlet distribution.

f

Incidence model 3: the log-normal model (CV = 1.23): where ${\pi }_i=c{a}_i, i=\mathrm{1,2},\dots ,S$, and $\left(a_1,a_2,\dots ,a_S\right)$ is a random sample from a log-normal distribution, and scale c is used to control the maximum of ${\pi }_i$.

The coefficient of variation (CV) in these six models ranged from 0 to 1.56, encompassing a wide range of values that encompass most practical cases in real-world applications.

In the simulation study, different sample sizes are considered to represent varying levels of sampling effort. For each simulation scenario, 1000 simulated datasets are generated. The estimates and their corresponding estimated standard errors (SE) are averaged across the 1000 simulated datasets to obtain the mean estimate and mean estimated SE. The sample SE and root-mean-square error (RMSE) are calculated based on the 1000 estimates to determine the sample SE and sample RMSE. The percentage of 95% confidence intervals (CIs) that cover the true value and the average observed richness are also calculated. All the simulation results are presented in Table 1 and Table 2. For simplicity, the average estimates of the discussed estimators are plotted in Figure 1 and Figure 2 to illustrate their statistical behavior as a function of sampling effort.

3.2. Simulation results for richness estimation of one assemblage (a local region)

In this case, the integrated dataset consists of three random samples that are independently collected from the same assemblage. Each sample is simulated separately based on one of the three discussed abundance/incidence models, representing different sampling situations or methods. Different sample sizes are considered to indicate varying levels of sampling efforts, ranging from n = 200 to 600 with an increment of 50 for abundance data and t = 10 to 50 with an increment of 5 for incidence data.

Four different scenarios are examined, including:

• Three abundance models: random uniform, broken-stick, and log-normal.
• Two abundance models: random uniform and broken-stick, along with one incidence model: log-normal.
• One abundance model: random uniform, along with two incidence models: broken-stick and log-normal.
• Three incidence models: random uniform, broken-stick, and log-normal.

The simulation results for these four scenarios are presented separately in Figure 1a–d and Table 1.

3.3. Simulation results for richness estimation of multiple assemblages (a large-scale region)

To evaluate the statistical behavior of the discussed estimators for regional richness estimation based on integrated data, the target region is divided into three subregions. The integrated data consists of three samples that are collected separately from each subregion. It is assumed that the target region comprises S = 600 species, with each subregion having 300 species. There are varying numbers of shared species and unique species in each subregion.

Different sample sizes are considered to represent different sampling efforts, ranging from n = 400 to 4000 with an increment of 450 for abundance data, and t = 10 to 50 with an increment of 5 for incidence data. Four different scenarios are examined, including:

• Three abundance models: random uniform, broken-stick, and log-normal.
• Two abundance models: random uniform and broken-stick, along with one incidence model: log-normal.
• One abundance model: random uniform, along with two incidence models: broken-stick and log-normal.
• Three incidence models: random uniform, broken-stick, and log-normal.

The simulation results for these four scenarios are presented separately in Figure 2a-d and Table 2.

Undoubtedly, for a fixed sample size, a superior species richness estimator should exhibit lower bias and variance (i.e., low RMSE). Additionally, the coverage rate of its associated 95% confidence interval should be close to 0.95. As the sample size increases, the estimator should demonstrate the following essential behaviors: its bias, accuracy (as measured by RMSE), and the coverage rate of its confidence interval should generally improve, eventually converging to the true species richness when the sample size is sufficiently large. Based on these criteria, the following findings can be concluded from the simulation results:

• In all simulation scenarios presented in Table 1 and Table 2 and Figure 1 and Figure 2, both Chao3 and the proposed adjusted Chao3 (Chao3_Adj) consistently provide robust lower bound estimates in all hypothetical models, and they tend to approach the true richness as the sample size increases.
• Both Chao3 and Chao3_Adj exhibit the essential statistical behaviors: their bias and RMSE decrease, resulting in more accurate 95% confidence intervals as the sample size increases (Table 1 and Table 2).
• The estimators of the discussed estimators' variance, derived using the asymptotic approach, perform well across all simulation scenarios (Table 1 and Table 2).
• Compared to Chao3, the adjusted estimator (Chao3_Adj) exhibits lower bias, larger standard errors, and more accurate 95% confidence intervals for the true richness in all simulation scenarios (Table 1 and Table 2).
• When samples are directly collected from the entire region (Table 1), Chao3_Adj has higher RMSEs compared to Chao3. However, when samples are separately collected from subregions within the target region (Table 2), Chao3_Adj demonstrates lower RMSEs.

These findings collectively demonstrate the favorable performance of the proposed adjusted Chao3 estimator (Chao3_Adj) in terms of bias, standard error, and accuracy of the 95% confidence interval, particularly when samples are collected from subregions within the target region.

3.4. Using data sets as true assemblages

I utilized two biological survey datasets representing true assemblages and generated separate datasets from these two assemblages. In each dataset, the observed species relative abundance was considered as the true species relative abundance or detection probability. A sample of size n (or t) was then generated through sampling with replacement to create the sampling dataset. Different sample sizes were considered to indicate varying levels of sampling efforts.

The average estimate and other relevant statistics obtained using the 1,000 generated datasets, as a function of sample size, are depicted in Figure 3 and Figure 4 and supplemental Table C1 and Table C2 (refer to Appendix C for detailed information). These evaluations aimed to assess the statistical behaviors of the discussed richness estimators across four different sampling scenarios: three abundance data, two abundance data and one incidence data, one abundance data and two incidence data, and three incidence data.

3.4.1. Moth species data

Moth species data were collected in the Golfo Dulce region of the Costa Rican rainforest from July to October 2014 [22]. The target region was divided into three types of forest: creek forest, slope forest, and ridge forest. Light traps were set up at 18 sites, with six replicates within each forest type. Further details can be found in [22].

Table 3 presents a summary of the data, including the sample size, observed richness, and the first five species frequency counts for each forest type. In the pooled sample, a total of 421 species were recorded, with 115, 285, and 356 species observed in the creek, slope, and ridge forests, respectively. In this case, the survey data sets are considered as the true assemblages. The proportion of species in the sample is assumed to represent the species' relative abundance for generating individual-based abundance data, while the ratio between species abundance and the maximum abundance is considered as the species' detection probability for generating sample-based incidence data. Consequently, each type of forest has its corresponding abundance model and incidence model. Four different scenarios are examined, including:

• Three abundance models: creek, slope, and ridge.
• Two abundance models: creek and slope, along with one incidence model: ridge.
• One abundance model: creek, along with two incidence models: slope and ridge.
• Three incidence models: creek, slope, and ridge.

The simulation results for each scenario are depicted separately in Figure 3a-3d and Supplementary Table C1 (Appendix C).

3.4.2. xylobiont beetle species data

The second real dataset comprises xylobiont beetle species data collected from the Leipzig floodplain forest in 2016 [23]. The beetle species data were collected separately from three dominant tree species: Quercus robur (QR), Tilia cordata (TC), and Fraxinus excelsior (FE) in the Leipzig floodplain forest area. Table 3 provides information on the sample size, observed richness, and the first three species frequency counts for each tree species. In total, 307 beetle species were observed, with 174, 198, and 184 species recorded in QR, TC, and FE tree species, respectively.

In this case, the survey data sets are treated as the true assemblages, and the species abundance/incidence model is constructed using the same method discussed earlier for each tree species. Four different scenarios are considered, including:

• Three abundance models: QR, TC, and FE.
• Two abundance models: QR and TC, along with one incidence model: FE.
• One abundance model: QR, along with two incidence models: TC and FE.
• Three incidence models: QR, TC, and FE.

The simulation results for each scenario are presented separately in Figure 4a-d and Supplementary Table C2 (Appendix C).

The results from the analysis, as depicted in Figure 3 and Figure 4 and Supplemental Tables in Appendix C, demonstrate that both Chao3 and the adjusted Chao3 (Chao3_Adj) effectively reduce the underestimation of observed richness. From a theoretical standpoint and considering the results of the simulation study, it is expected that Chao3_Adj exhibits lower bias compared to Chao3, particularly when there is high heterogeneity as indicated by a high coefficient of variation (CV). Furthermore, the supplemental Table C1 and Table C2 in Appendix C confirm that Chao3_Adj has lower bias, higher standard error, and lower root-mean-square error (RMSE). The higher estimated standard error of Chao3_Adj compared to Chao3 suggests that the former estimator may provide a more accurate 95% confidence interval for true richness. This observation aligns with the findings of the simulation study presented in Table 1 and Table 2.

Overall, the results support the notion that Chao3_Adj has less bias and performs better in terms of standard error and RMSE, indicating its potential to provide more accurate estimates and confidence intervals for true richness, particularly in scenarios with higher heterogeneity.

4. Discussion and Conclusion

Estimating species richness for a large-scale region or multiple assemblages poses a statistical challenge due to the difficulty of obtaining a random sample from the entire region. Typically, data collected for assessing species richness in such cases consist of multiple samples that are individually sampled from each assemblage or subregion. Additionally, these samples may employ different sampling schemes or strategies. As a result, the pooled sample cannot be considered a random sample from the entire region, even though each individual sample is randomly collected from its respective subregion or assemblage. Consequently, the pooled sample cannot be modeled using a traditional sampling distribution model.

However, certain richness estimators that rely only on the frequency counts of rare species in the sample have been theoretically demonstrated to be applicable to the pooled sample, as long as the samples are randomly collected and the sample size is not excessively small. In this context, the detection probability of a species may vary across the samples, and the data format (individual-based abundance data or sample-based incidence data) can differ among the samples.

In this research, Chao's lower bound estimator (Chao3) and its bias-corrected estimator (Chao3_Adj) are theoretically proven to be suitable for estimating richness in multiple assemblages. Chao3 is derived using the Cauchy-Schwarz inequality, while the adjusted estimator (Chao3_Adj) corrects the bias of Chao3 based on the Good-Turing frequency formula[14]. Simulation results demonstrate that both Chao3 and Chao3_Adj can be used to estimate regional richness based on the pooled sample from integrated data, aligning with the theoretical findings. These estimators provide lower bound estimates in all hypothetical models and tend to converge to the true richness as the sample size increases. Notably, when the sample size is small or the community exhibits high heterogeneity, Chao3_Adj outperforms Chao3 with lower bias, lower root-mean-square error (RMSE), higher standard error (s.e.), and a more accurate 95% confidence interval (CI) for true richness.