A Bayesian theory for estimation of biodiversity

In this section, we provide an overview of Gibbs-type priors (Gnedin and Pitman, 2005), focusing on their relevance for ecological applications and quantification of biodiversity. For a mathematical exposition, we refer to De Blasi et al. (2015), while key theoretical developments are described in Lijoi et al. (2007a, b, 2008a, 2008b); De Blasi et al. (2013). Our data set comprises taxon labels $X_{1},\dots,X_{n}$, representing, for example, species or families, and taking values in a set $\mathds{X}$. We assume that $(X_{n})_{n\geq 1}$ is an (ideally) infinite sequence of exchangeable observations which, according to the de Finetti theorem (de Finetti, 1937), is equivalent to the following hierarchical specification

where $\tilde{p}$ is a (random) probability measure whose law $\mathcal{Q}$ is the prior distribution in Bayesian statistics. In our setting, we specify a multinomial model for $\tilde{p}$, namely

where $\delta_{x}$ represents the Dirac delta measure at $x$, whereas $\pi_{h}$ are random probability weights such that $\sum_{h=1}^{\infty}\pi_{h}=1$ almost surely, and $(Z_{h})_{h\geq 1}$ are random $\mathds{X}$-valued locations, denoting distinct taxon labels. This approach is referred to as nonparametric because there are potentially infinitely many probabilities $(\pi_{h})_{h\geq 1}$, for which a suitable prior will be specified. We assume that the $Z_{h}$’s are iid samples from a diffuse probability distribution $P$ on $\mathds{X}$ and that $(Z_{h})_{h\geq 1}$ and $(\pi_{h})_{h\geq 1}$ are independent. The diffuseness of $P$ implies that the labels $Z_{h}$’s are all different. Thus, $\tilde{p}$ belongs to the class of proper species sampling models (ssms, Pitman, 1996). The baseline measure $P$ serves as a technical tool that simplifies the mathematical exposition, but we are not interested in “learning” it. Indeed, without loss of generality, in species sampling problems, we may let $\mathds{X}=\mathds{R}$ and let $P$ be a uniform distribution, meaning that each value $Z_{h}$ is the numerical encoding of the associated taxon.

The primary focus of inference typically lies not in the probabilities $(\pi_{h})_{h\geq 1}$ but rather in the combinatorial structure induced by the multinomial model of equation (2). Due to the discrete nature of $\tilde{p}$, there will be identical values among $X_{1},\dots,X_{n}$ with positive probability, comprising a total of $K_{n}=k$ distinct values with labels $X_{1}^{*},\dots,X_{k}^{*}$ and frequencies $n_{1},\dots,n_{k}$ such that $\sum_{j=1}^{k}n_{j}=n$. The random variable $K_{n}$ represents the observed taxon richness, whereas $n_{1},\dots,n_{k}$ denote the associated abundances. The ties among the observations $X_{1}\dots,X_{n}$ induce a random partition $\Psi_{n}=\{C_{1},\dots,C_{k}\}$ of the indices $\{1,\dots,n\}$, where $C_{j}=\{i:X_{i}=X_{j}^{*}\}$ and $n_{j}=|C_{j}|$. A species sampling model $\tilde{p}$ is a Gibbs-type process if the law of the random partition, called Gibbs partition, is such that

where $(a)_{n}=a(a+1)\cdots(a+n-1)$ denotes a rising factorial, $\sigma<1$ is a discount parameter, and $V_{n,k}$ are non-negative weights satisfying the forward recursion $V_{n,k}=\left(n-\sigma k\right)V_{n+1,k}+V_{n+1,k+1}$ for any $n\geq 1$ and $1\leq k\leq n$. $\Pi_{n}\left(n_{1},\ldots,n_{k}\right)$ is called the exchangeable partition probability function (eppf) and, in conjunction with $P$, characterizes the distribution of random probability measure $\tilde{p}$ (Pitman, 1996). Of considerable importance is the predictive distribution induced by $\Pi_{n}$ (Lijoi et al., 2007a):

The above predictive rule provides a sampling mechanism for the data $X_{1},\dots,X_{n}$. The term $\mathds{P}(X_{n+1}=\text{``new''}\mid X_{1},\dots,X_{n})=V_{n+1,k+1}/V_{n,k}$ denotes the probability of discovering a new taxon, i.e., the probability of $X_{n+1}$ not being observed among the previous $X_{1},\dots,X_{n}$. Conversely, the probability $\mathds{P}(X_{n+1}=\text{``old''}\mid X_{1},\dots,X_{n})=1-V_{n+1,k+1}/V_{n,k}$ corresponds to a Bayesian estimate of the sample coverage, namely the fraction of already observed taxa. The probability of re-observing the $j$th taxon $X^{*}_{j}$ is proportional to $n_{j}-\sigma$, explaining why $\sigma<1$ is sometimes referred to as the “discount parameter”, as it diminishes (or augments) the observed frequencies. Notable Gibbs-type priors are discussed below.

Remarkably, the aforementioned set of weights $V_{n,k}(\sigma,H)$, $V_{n,k}(\alpha)$ and $V_{n,k}(\sigma,\gamma)$ form the foundation of any Gibbs-type prior. In fact, Gibbs partitions can always be expressed as a mixture with respect to the parameters $H$, $\alpha$, and $\gamma$. This result was established by Gnedin and Pitman (2005).

Any Gibbs-type process can be represented hierarchically, involving a suitable prior distribution for the key parameters $H$, $\alpha$, and $\gamma$. Prior distributions for $H$ in the $\sigma<0$ case are discussed in Gnedin (2010) and De Blasi et al. (2013); see also Miller and Harrison (2018) for applications to mixture models and Legramanti et al. (2022) for employment in stochastic block models. In the $\sigma=0$ case, a popular choice is the semi-conjugate Gamma prior for $\alpha$ as in Escobar and West (1995), with the Stirling-gamma prior of Zito et al. (2024) a recent alternative. Finally, in the $\sigma\in(0,1)$ regime, a polynomially tilted stable distribution prior for $\gamma^{-1/\sigma}$ leads to the Pitman-Yor process (Pitman and Yor, 1997), whereas an exponentially tilted stable density leads to the normalized generalized gamma process (Lijoi et al., 2007b); refer to Lijoi et al. (2008b); Favaro et al. (2009) for a thorough investigation and to Favaro et al. (2015) for a general class of priors for $\gamma$.

The simplest measure of biodiversity is arguably the taxon richness. In our notation, we observe $K_{n}=k$ different taxa among $n$ data points $X_{1},\dots,X_{n}$; we refer to this value as the in-sample richness. A priori, the distribution of $K_{n}$ induced by a Gibbs-type prior has a simple form

where $\mathscr{C}(n,k;\sigma)$ denotes a generalized factorial coefficient (Charalambides, 2002); the case $\sigma=0$ is recovered as the limit of the above formula, since $\lim_{\sigma\rightarrow 0}\mathscr{C}(n,k;\sigma)\sigma^{-k}=|s(n,k)|$, which is the signless Stirling number of the first kind. The expected values $\mathds{E}(K_{1}),\dots,\mathds{E}(K_{n})$ can be interpreted as a model-based rarefaction curve (Zito et al., 2023). Moreover, we may want to predict the number of new taxa $K_{n+m}\geq k$ within a future sample $X_{n+1},\dots,X_{n+m}$; we call this value the out-of-sample richness. The posterior distribution of $K_{n+m}$ given the data $K_{n}=k$, or equivalently of the number of previously unobserved taxa $K_{m}^{(n)}=K_{n+m}-K_{n}$, has been derived by Lijoi et al. (2007a) and is

where the term $\mathscr{C}(m,j;\sigma,-n+k\sigma)$ is the noncentral generalized factorial coefficient; see Lijoi et al. (2007a). The collection $\mathds{E}(K_{n+1}\mid K_{n}=k),\dots,\mathds{E}(K_{n+m}\mid K_{n}=k)$ represents a model-based extrapolation of the accumulation curve. Note that $K_{n}=k$ is a sufficient statistic for predictions.

While rarefaction and extrapolation curves are valuable tools (Gotelli and Colwell, 2001), it is useful to summarize biodiversity with a single number. The concept of richness, defined as $\lim_{n\rightarrow\infty}K_{n}$, is problematic since $\lim_{n\to\infty}K_{n}=\infty$ for $\sigma\geq 0$ regardless of the observed data, whereas $K_{n}$ remains finite for $\sigma<0$. Avoiding $\sigma\geq 0$ is tempting, but leads to poor fit and predictions for some datasets. Gibbs-type priors with positive $\sigma\geq 0$ often excel in predicting future values $K_{m}^{(n)}$ for highly diverse taxa, compared to models with $\sigma\leq 0$ (Lijoi et al., 2007a; Favaro et al., 2009). Additionally, the total number of taxa present in a given area can be estimated even in models where $\sigma\geq 0$, e.g. following the strategy of ter Steege et al. (2013). This discussion motivates embracing an alternative concept of diversity, called $\sigma$-diversity, introduced by Pitman (2003), which encompasses Hubbell (2001)’s fundamental biodiversity number when $\sigma=0$.

The $\sigma$-diversity can be seen as a richness measure that has been appropriately rescaled. Proposition 2 highlights the central role of the Dirichlet-multinomial, Dirichlet, and $\sigma$-stable pk processes among Gibbs-type priors. It shows that their parameters $H$, $\alpha$, and $\gamma$ are biodiversity indices. In these three processes, $\sigma$ diversity is deterministic and assumed to be known. However, $\sigma$-diversities $H$, $\alpha$, or $\gamma$ are typically unknown, and can be estimated employing a prior distribution, leading to a Gibbs-type process for $\tilde{p}$ thanks to Proposition 1. In light of this, the posterior law of $\sigma$-diversity is a key quantity for measuring biodiversity. In addition, the posterior distribution of $S_{\sigma}$ has an elegant connection with accumulation curves, as shown in the following theorem.

Theorem 1 states that the posterior law of $S_{\sigma}$ coincides with the $\sigma$-diversity associated with extrapolation of the accumulation curve. Related results for the normalized generalized gamma model are in Favaro et al. (2012). In practice, deriving the posterior distribution of $\sigma$-diversities $H$, $\alpha$ and $\gamma$ is based on the Bayes theorem, where the eppf in (3) acts as likelihood function. Notably, in Gibbs-type priors, the abundances $n_{1},\dots,n_{k}$ appearing in (3) do not provide insights about the diversity, and they do not appear in the posterior for $S_{\sigma}$ because the number of observed species $k$ is a sufficient statistic. This would not be the case for general species sampling models beyond the Gibbs-type.

We show that many classical biodiversity indices can be expressed in terms of $H$, $\alpha$, or $\gamma$, based on the connection with the accumulation curves discussed earlier. Consider the popular Simpson similarity index $\mathcal{S}$ (e.g., Colwell, 2009), which is the probability that two randomly chosen individuals belong to the same taxa. Within the framework of species sampling models, $\mathcal{S}:=\sum_{h=1}^{\infty}\pi_{h}^{2},$ and due to exchangeability $\mathds{E}(\mathcal{S})=\mathds{P}(X_{1}=X_{2}).$ In Gibbs-type priors, it follows that $\mathds{E}(\mathcal{S})=V_{2,1}$. We can specialize this formula in the Dirichlet multinomial and Dirichlet process case, yielding respectively

Thus, there is an inverse relationship between the expected Simpson similarity and the $\sigma$-diversity. The following section shows a similar relationship in the polynomial regime, when $\sigma=1/2$. More generally, most biodiversity indices can be expressed as $\mathcal{I}_{g}:=\sum_{h=1}^{\infty}g(\pi_{h}),$ for some function $g:(0,1)\rightarrow\mathds{R}^{+}$. For instance, $g(\pi)=\pi^{2}$ corresponds to the Simpson index, $g(\pi)=-\pi\log{\pi}$ is the Shannon diversity, and $g(\pi)=\pi^{\tau}$ for some $\tau>0$ is the (unnormalized) Tsallis or Hill diversity (Colwell, 2009; Magurran and McGill, 2011). The expectations $\mathds{E}(\mathcal{I}_{g})$ are a function of $H,\alpha$ and $\gamma$. Additionally, their relationship with $\mathcal{I}_{g}$ can often be explicitly determined: if the random probabilities $(\pi_{h})_{h\geq 1}$ are in size-biased order, then we have the simplification $\mathds{E}(\mathcal{I}_{g})=\mathds{E}\left\{\sum_{h=1}^{\infty}g(\pi_{h})\right\}=\mathds{E}\{g(\pi_{1})/\pi_{1}\}$. For example, in the dp case, the stick-breaking weights are in size-biased order, and we have $\pi_{1}\sim\text{Beta}(1,\alpha)$, which makes computations straightforward.

We propose two approaches to assess the fit of the model. The first approach compares the observed distinct values $K_{1},\dots,K_{n}$ with the model-based rarefaction curve $\mathds{E}(K_{1}),\dots,\mathds{E}(K_{n})$. As shown above, expectations $\mathds{E}(K_{i})$ can often be computed explicitly, and are typically evaluated given an estimate for the diversity parameter, e.g. $\mathds{E}(K_{i}\mid\hat{\alpha}_{\textsc{ml}})$ in the $\sigma=0$ case. Due to exchangeability, expectations $\mathds{E}(K_{i})$ do not depend on the order of the data. However, in most cases, the data $X_{1},\dots,X_{n}$ are not observed in a specific order. Thus, a common practice is to reshuffle the order of $X_{1},\dots,X_{n}$ and then consider the averages $\bar{K}_{1},\dots,\bar{K}_{n}$, over all possible permutations of the data. Conveniently, for this combinatorial problem, there exists an explicit solution, which is:

where $K_{n}=k$ and $n_{1},\dots,n_{k}$ are the abundances; see Smith and Grassle (1977); Colwell et al. (2012). The values $\bar{K}_{i}$ define the “classical rarefaction”, which can be regarded as a frequentist nonparametric estimator arising from a multinomial model. Our Bayesian nonparametric approach instead depends on a few parameters, imposing some rigidity on the functional shape of the rarefaction, which is critical for extrapolation. Graphically comparing $\bar{K}_{1},\dots,\bar{K}_{n}$ with $\mathds{E}(K_{1}),\dots,\mathds{E}(K_{n})$ can give a sense of the suitability of the chosen model. We show a concrete example in Section 5.1.

The second model-checking approach is even simpler and has been used, e.g. in Favaro et al. (2021); see also Thisted and Efron (1987) for early ideas. The abundances $n_{1},\dots,n_{k}$ can be reformulated in terms of frequency counts $m_{1},\dots,m_{n}$, where $m_{r}$ is the number of taxa appearing with frequency $r$ in the sample. Hence, $m_{1}$ represents the number of singletons, $m_{2}$ is the number of doubletons, etc. We denote by $M_{1,n},\dots,M_{n,n}$ the associated random variables. To assess goodness of fit, we compare empirical counts $m_{1},\dots,m_{n}$ with their model-based expectations $\mathds{E}(M_{1,n}),\dots,\mathds{E}(M_{n,n})$. In the Dirichlet process case, these expectations have a simple analytical formula

and therefore when $r=1$, corresponding to the expected singletons implied by Hubbell (2001) theory, we get $\mathds{E}(M_{1,n}\mid\alpha)=n\alpha/(n+\alpha)=n\>\text{pr}(X_{n+1}=\text{``new''}\mid\alpha,X_{1},\dots,X_{n})$. The general analytical formulas for the expectations $\mathds{E}(M_{r,n})$ for any prior of Gibbs type are given in Favaro et al. (2013). These expectations can be approximated via Monte Carlo, by sampling $X_{1},\dots,X_{n}$ from the urn scheme in equation (4).