Species Abundance Distribution and Species Accumulation Curve: A General Framework and Results

Estimating the diversity of classes in a population is a problem encountered in many fields. We may be interested in the diversity of words a person know from his/her writings [18], the illegal immigrants from the apprehension records [3], the distinct attributes in a database [24, 16], or the distinct responses to a crowdsourcing query [48]. Among different applications, species abundance is the one that receives most attention. For this reason, it is chosen as the theme of this paper with the understanding that the proposed framework and methods are applicable in other applications as well.

Understanding the species abundance in an ecological community has long been an important task for ecologists. Such knowledge is paramount in conservation planning and biodiversity management [35]. An exhaustive species inventory is too labor and resource intensive to be practical. Information about species abundance can thus be acquired mainly through a survey.

Let $N=(N_{0},N_{1},...)$ where $N_{k}$ is the number of species in a community that are represented exactly $k$ times in a survey. We do not observe the whole $N$, but $\tilde{N}=(N_{1},N_{2},...)$ which is the zero-truncated $N$. In other words, we do not know how many species are not seen in the survey. We call the vector $\tilde{N}$, the frequency of frequencies (FoF) [21].

A plethora of species abundance models have been proposed for $\tilde{N}$. Comprehensive review of the field can be found in Bunge and Fitzpatrick [6] and Matthews and Whittaker [34]. A typical assumption in purely statistical models is

where $N_{+}=\sum_{k=1}^{\infty}N_{k}$ is the total number of recorded species, and $p=(p_{1},p_{2},\ldots)$ is a probability vector, such that $p_{k}$ is the probability for a randomly selected recorded species to be observed $k$ times in the survey for $k=1,2,\ldots$. We call this vector $p$, the species abundance distribution (SAD). The great significance of (1.1) relies on three assumptions: (A1) the species names are noninformative; (A2) the observed data for different recorded species are independent and identically distributed (iid), and (A3) for each recorded species, its observed frequency contains all useful information.

Species accumulation curve (SAC) is another popular tool in the analysis of species abundance data. The survey is viewed as a data-collection process in which more and more sampling effort is devoted. The individual-based SAC is the number of recorded species expressed as a function of the amount of sampling effort.

Despite the different emphases of SAD and SAC, the two approaches have an overlapped target: Estimating $D$, the total number of species in the community. For SAD, it means estimating $N_{0}$, the number of unseen species as $D=N_{0}+N_{+}$. For SAC, $D$ is the total number of seen species when the sampling effort is unlimited.

As SAD depends largely on the sampling effort, it is necessary to include sampling effort explicitly in the model in order to make comparison of different SADs possible. This addition establishes a link between SAD and SAC. Sampling effort can be of continuous type, such as the area of land or the volume of water sampled, or the duration of the survey. Discrete type sampling effort can be the sample size. To emphasize the sampling effort considered, the SAC is called species-time curve, species-area curve, or species-sample-size curve when the sampling effort is time, area, or sample size respectively. Species-area curve has been studied extensively in the literature. Review of it can be found in Tjørve [46, 47], Dengler [15] and Williams et al. [50]. In this paper, we use time as the measure of sampling effort. We consider the vector $N$ to be a function of the time $t$, denoted as $N(t)$. Notations $N_{k}$, $\tilde{N}$, $N_{+}$, $p$ and $p_{k}$ are likewise denoted as $N_{k}(t)$, $\tilde{N}(t)$, $N_{+}(t)$, $p(t)$ and $p_{k}(t)$ respectively. The (empirical) SAC is $N_{+}(t)$. Throughout the paper, we assume without loss of generality that the survey starts at time $0$ and ends at time $t_{0}>0$. Because of the great similarity of the species-time relationship and species-area relationship [39, 32], time and area are treated as if they are interchangeable in this paper. When SAC is a species-area curve, we refer to the Type I species-area curve [43], where the areas sampled are nested (smaller areas are included in larger areas), resulting in a curve that is always nondecreasing.

The study of the relation between the (empirical) SAD and the (empirical) SAC started early since the introduction of rarefaction curve ([1, 42]) which describes how we interpolate SAC using the observed FoF. Let $n_{k}(t_{0})$ be the observed $N_{k}(t_{0})$ for $k=1,2,\ldots$, $n_{+}(t_{0})=\sum_{k=1}^{\infty}n_{k}(t_{0})$, and $\tilde{n}(t_{0})=\{n_{k}(t_{0})\}_{k\geq 1}$. The rarefaction curve for species-time relationship is

The rationale of (1.2) is that $1-(1-t/t_{0})^{k}$ is the probability for a species with frequency $k$ in time interval $[0,t_{0}]$ to be observed before time $t$ if the appearance times of it in $[0,t_{0}]$ are iid $U[0,t_{0}]$ distributed. Equivalent formula for species-area curve appears earlier in Arrhenius [1]. Good and Toulmin [22] proposed an extrapolation formula for species-sample-size curve which when expressed as species-time curve is

Equation (1.3) is (1.2) with a change of domain. Define the empirical SAD, $\hat{p}_{k}(t)=n_{k}(t)/n_{+}(t)$. Equation (1.3) becomes

where $\hat{h}_{t}(s)$ is an estimator of $h_{t}(s)$, the probability generating function of the SAD at time $t$. Equation (1.4) delineates a simple and yet convincing formula to find SAC from SAD and vice versa. The aim of this paper is to propose and study a statistical framework under which model in (1.1) and the relation in (1.4) hold when $\hat{N}_{+}(t)$, $n_{+}(t_{0})$ and $\hat{h}_{t_{0}}(s)$ are replaced by $E(N_{+}(t))$, $E(N_{+}(t_{0}))$ and $h_{t_{0}}(s)$ respectively. Our framework assumes (A1), (A2) and that the appearance times of each species follow a Poisson process which is sufficient for (A3) and the validity of (1.2).

Finite $D$ is commonplace in SAD approach although no consensus has been reached. Empirically, log-series distribution, an SAD that assumes infinite $D$, is one of the most successful models. Many SACs do not have asymptote. It is common that the observed number of rare species is large, and shows no sign to decrease. In Bayesian nonparametric approach in genomic diversity study, Poisson-Dirichlet process which assumes infinite $D$ is used (see for example Lijoi et al. [31]). In reality, the existence of transient species (species which are observed erratically and infrequently [38]), and the error in the species identification process (a well-known example is the missequencing in pyrosequencing of DNA (see for example Dickie [17])) are continual sources of rare species making the species number larger than expected. In our framework, $D$ is random and can be finite or infinite with probability one. We use species richness to refer to $E(D)$ when $D$ is random, and $D$ when $D$ is deterministic.

A special feature of our framework is that species can have zero detection probability. We call such species, zero-rate species, and all other species, positive-rate species. Zero-rate species can either be seen only once or unseen in a survey. If there are finite number of zero-rate species, the probability of observing any of them is zero. Therefore, if zero-rate species is observed in a survey, almost surely there are infinite number of them. We interpret observed individuals of positive-rate species as outcomes of a discrete distribution, where individuals of the same species can appear any nonnegative number of times in a survey. On the other hand, individuals of zero-rate species are outcomes of a continuous distribution, and no two such individuals belong to the same species. In our framework, the distribution is allowed to be a mixture of the two. Suppose we want to estimate the population of a town through recording each person we meet on street. Then tourists from distant countries can be viewed as “zero-rate species”. In the first example in Section 9, it is suspected that the sequencing error in pyrosequencing of DNA may be a source for zero-rate species.

Poisson distribution is a main component in our framework. Though Poisson distribution is common in existing SAD models, time $t$ scarcely plays a role. In mixed Poisson model (see for example Fisher et al. [20], Bulmer [5]), the observed frequencies of the species are independent and each follows a Poisson distribution with its own rate, say $\lambda_{i}$ for species $i$. The value $D$ is a fixed unknown finite value and $\{\lambda_{i}\}$ are iid sample from a mixing distribution. Another related model was proposed in Zhou et al. [51]. The paper focuses mainly on finite $D$. Neither time nor zero-rate species are included in the model.

The outline of this paper is as follows. We propose in Section 2 a new model for the sampling process, called the mixed Poisson partition process (MPPP). We emphasize on the parametric approach where additional assumptions are made on top of the framework so that the process depends only on a few parameters. Once the parameters are estimated, we can make inference on different characteristics of the population, say the SAD at any fixed time, the Hill numbers, or the expected future data. We study $\tilde{N}(t)$ in Section 3. In Section 4, we consider the expected species accumulation curve (ESAC). We prove the one-to-one correspondence among (i) an MPPP, (ii) an ESAC which is a Bernstein function that passes through the origin, and (iii) the expected number of recorded species and the SAD at time $t_{0}$ such that the first derivative of its probability generating function is absolutely monotone in $(-\infty,1)$. In Section 5, we introduce $\mathrm{LDR}_{1}$, a parametric family of ESAC. The SAD of $\mathrm{LDR}_{1}$ is the Engen’s extended negative binomial distribution. This family has an attractive property that the ratio of the first and the second derivatives of the ESAC is a linear function of $t$. We extend $\mathrm{LDR}_{1}$ to $\mathrm{LDR}_{2}$ which allows zero-rate species. In Section 6, a D1/D2 plot for $\mathrm{LDR}_{1}$ and some diagnostic plots for specified $\mathrm{LDR}_{1}$ distributions are proposed. Extrapolation of the curve in D1/D2 plot is considered in Section 7. Estimator of species richness is suggested basing on a modified first-order extrapolation. In Section 8, $\mathrm{LDR}_{1}$ is generalized so that the derivative ratio is a rational function instead of a linear function of $t$. Four real data are analyzed in Section 9 to demonstrate the applications of the proposed models and the suggested plots. In Section 10, we propose and study a design where only a few leading appearance times of each species are recorded. In Section 11, we consider the maximum likelihood approach on the empirical SAC. We give a discussion in Section 12.

Poisson process is the backbone of the framework. A Poisson process is a point process characterized by an intensity measure over the $n$-dimensional space $\mathbb{R}^{n}$ (we usually have $n=1$ in this paper). The intensity measure which we denote as $\omega$ delineates how many points are present on average in different parts of $\mathbb{R}^{n}$. More precisely, the number of points in a set $A\subseteq\mathbb{R}^{n}$ follows the $\mathrm{Poisson}(\omega(A))$ distribution. Furthermore, for any finite collection of disjoint subsets $A_{1},\ldots,A_{k}\subseteq\mathbb{R}^{n}$, $S_{1},\ldots,S_{k}$ are mutually independent, where $S_{i}$ is the number of points in $A_{i}$. If $\omega(dx)=f(x)dx$, we call $f(x)$ the intensity function. A simple example is the homogeneous Poisson process where $\omega(dx)=\lambda dx$ for a constant $\lambda$. The parameter $\lambda$ is called the rate of the process.

If $\omega$ is finite (i.e., $\int\omega(dx)<\infty$), simulation of a Poisson process can be performed in two steps: (i) simulate the total number of points $W\sim\mathrm{Poisson}(\int\omega(dx))$, and (ii) simulate $X_{1},\ldots,X_{W}$ iid from the probability measure $\omega/\int\omega(dx)$. Nevertheless, if $\omega$ is infinite, then the number of points is infinite, and it is impossible to simulate all the points in the process. In this case, we can only simulate a selected finite subset of points of the process using the thinning property of Poisson processes [29]. For each point $x$ in the Poisson process, let $\alpha(x)\in[0,1]$ be the probability for the point $x$ to be selected, and $1-\alpha(x)$ be the probability for it to be discarded. Then the selected points form a Poisson process with intensity measure $\alpha\omega$, where $\alpha\omega(A)=\int_{A}\alpha(x)\omega(dx)$. We can simulate the selected points when the measure $\alpha\omega$ is finite using the aforementioned method. For example, if we are interested only in the points that lie in a bounded region $A\subset\mathbb{R}^{n}$ of a homogeneous Poisson process with rate $\lambda$, then $\alpha(x)=\mathbf{1}\{x\in A\}$ where $\mathbf{1}\{.\}$ is the indicator function. The number of selected points follows $\mathrm{Poisson}(\lambda\int_{A}dx)$ distribution, and the selected points are iid uniformly distributed in $A$.

We model the observations in a species abundance survey by a stochastic process where rates of the species follows a Poisson process.

Measures $\nu$ and $\tilde{\nu}$ may be finite or infinite. Let $\Lambda=\int\nu(d\lambda)$. As the expected number of individuals seen in time interval $[0,t]$ is equal to $t\Lambda$ (see (3.4) in Section 3), we can interpret $\Lambda$ as the expected total rate. When $\nu$ is finite (i.e., $\Lambda<\infty$), conditional on the event that there is an individual observed exactly at time $t$, the distribution of the rate of the species of that individual is $\nu/\Lambda$ (see Proposition A in Appendix A for a proof). Therefore, we can regard $\nu$ as the intensity measure of $\lambda$ of the observed individuals. If $\tilde{\nu}$ is finite, the expected number of positive-rate species in the community is finite. From Step 3 in the definition, if $\nu(\{0\})>0$, there are infinite number of zero-rate species and they arrive at a constant rate. With probability one, $D$ is finite if and only if $\nu(\{0\})=0$ and $\tilde{\nu}(\mathbb{R}_{>0})<\infty$. In such case, $D$ follows a $\mathrm{Poisson}(\tilde{\nu}(\mathbb{R}_{>0}))$ distribution. Measure $\tilde{\nu}$ specifies the distribution of the rates of positive-rate species. More precisely, $\tilde{\nu}([\lambda_{0},\lambda_{1}])$ with $\lambda_{0}>0$ is the average number of species with rate in $[\lambda_{0},\lambda_{1}]$. Measure $\nu$ specifies the distribution of the rates of the species of the individuals including those of the zero-rate species. More precisely, $\nu([\lambda_{0},\lambda_{1}])$ is the average number of individuals (per unit time) belonging to species whose rates lie in $[\lambda_{0},\lambda_{1}]$. Condition (2.1) is essential because it is equivalent to the finiteness of the ESAC (i.e., $E(N_{+}(t))<\infty$ for any finite nonnegative $t$). A proof of it is given in Appendix B. Figure 1 illustrates the generation of the MPPP.

When $\tilde{\nu}$ is infinite, we cannot simulate all $\lambda_{i}$’s in Step 1 in practice. We can use the following method to simulate a realization of the process in interval $[0,t_{0}]$. The probability for a species with rate $\lambda$ to be recorded in $[0,t_{0}]$ is $1-\exp(-\lambda t_{0})$. Applying the thinning property, the intensity of the recorded positive-rate species is $(1-\exp(-\lambda t_{0}))\tilde{\nu}$ which is always finite. To include also the recorded zero-rate species, the intensity is $(1-\exp(-\lambda t_{0}))\lambda^{-1}\nu$ (we take $(1-\exp(-\lambda t_{0}))\lambda^{-1}=t_{0}$ when $\lambda=0$). After generating $\lambda_{1},\lambda_{2},\ldots$ according to a Poisson process with this intensity measure, we simulate for each $i$, a realization $\eta_{i}$ (independently across $i$) of a Poisson process with rate $\lambda_{i}$, conditional on the event that $\eta_{i}$ has at least one point in $[0,t_{0}]$ (if $\lambda_{i}=0$, then $\eta_{i}$ contains one uniformly distributed point in $[0,t_{0}]$). An alternative equivalent definition that unifies the generation of individuals of zero-rate species and positive-rate species is given in Appendix C.

A special feature of the framework is that $D$ is random and can be infinite with probability one. This change necessitates modification of biodiversity measures, among which Hill numbers are popular. When $D$ is deterministic, Hill number of order $q$ [25] for $q\geq 0$ and $q\neq 1$ is ${}^{q}\!D=\Big{(}\sum_{i}(p_{\mathrm{sp}}(i))^{q}\Big{)}^{1/(1-q)}$ where $p_{\mathrm{sp}}(i)$ is the relative abundance of species $i$ in an assemblage (the number of individuals of species $i$ divided by the total population). When $q=1$, ${}^{1}\!D$ is defined as $\lim_{q\to 1}\,^{q}\!D=\exp(-\sum_{i}p_{\mathrm{sp}}(i)\log(p_{\mathrm{sp}}(i)))$. Hill numbers are interpreted as the effective number of species. The order $q$ determines the sensitivity of the measure to species relative abundances. When $q=0$, all species regardless its abundance are treated equally. Larger $q$ imposes more weight to abundant species. Hill numbers encompass three important diversity measures as its special cases. They are the species richness ($q=0$), the exponential of Shannon entropy ($q=1$), and the inverse of Simpson index ($q=2$). A notable property of Hill numbers is that they obey the replication principle: The diversity of an assemblage formed by pooling $m$ equally abundant and equally large assemblages with no species in common is $m$ times the diversity of a single assemblage.

Under our framework, each species corresponds to a Poisson process, and the relative abundance of a species is its rate divided by the expected total rate, $\Lambda$. Reasonable modifications to the definition of Hill numbers are to replace $p_{\mathrm{sp}}(i)$ with $\lambda/\Lambda$, where $\lambda$ is the rate of a species, and to replace summation over species with integration with respect to the measure $\lambda^{-1}\nu$ so that the integral of $\lambda/\Lambda$ is one. It works well when $\Lambda$ is finite, but fails when $\Lambda$ is infinite. To fix the problem, we need a meaningful surrogate rate which fulfills two requirements: (i) the expected total rate is finite, and (ii) it approaches the true rate $\lambda$ as a limit. The first appearance time of a species with rate $\lambda$ follows the exponential distribution with density function $\lambda\exp(-\lambda t)$, which can be regarded as the instantaneous rate of its first appearance at time $t$. This rate approaches $\lambda$ when $t$ decreases to zero. The expected total instantaneous rate of first appearances over all species at time $t$ is $\Lambda_{t}=\int\exp(-\lambda t)\nu(d\lambda)$ ($=E(N_{1}(t))/t$ from (3.1)) which is always finite for positive $t$. Replacing $p_{\mathrm{sp}}(i)$ with $\lambda\exp(-\lambda t)/\Lambda_{t}$ and taking $t\to 0$, we define, the Hill numbers, ${}^{q}\!D_{\nu}$ for our framework as

When $\nu(\{0\})>0$ and $0\leq q\leq 1$, ${}^{q}\!D_{\nu}=\infty$. When $\Lambda$ is finite,

Diversity ${}^{q}\!D_{\nu}$ is non-increasing with respect to $q$. Unlike the classical Hill numbers, ${}^{q}\!D_{\nu}$ can be less than one. For instance, from (3.3), ${}^{0}\!D_{\nu}=E(D)$ which can be any positive value. It can be proved that for any positive constant $\gamma$, ${}^{q}\!D_{\nu}=\gamma(^{q}\!D_{\nu/\gamma})$, an analogue of the replication principle for Hill numbers.

A sufficient statistic for a realization $G$ of an MPPP in time interval $[0,t_{0}]$ is $\tilde{N}(t_{0})$. Consider a time interval $[0,t]$. For each $\lambda$ in Step 2 of the definition, it contributes one to one of the counts $N_{k}(t)$, where $k$ follows a Poisson distribution, $\mathrm{Poisson}(\lambda t)$. By the splitting property of Poisson processes [29], their total contribution to $N_{k}(t)$ follows $\mathrm{Poisson}(\int(\lambda t)^{k}(k!)^{-1}e^{-\lambda t}\tilde{\nu}(d\lambda))$ distribution independent across $k$. The zero-rate species only increase $N_{1}(t)$ by a $\mathrm{Poisson}(\nu(\{0\})t)$ random variate. Therefore, for $k\geq 1$,

(we use the convention that $0^{0}=1$). Equations (3.2) and (3.3) follows from $(\ref{eq:Enk})$ for $k\geq 1$. The correctness of (3.1) when $k=0$ follows from (3.2), (3.3) and the fact $E(N_{0}(t))=E(D)-E(N_{+}(t))$.

Again for (3.2), we set $(1-\exp(-\lambda t))/\lambda=t$ when $\lambda=0$. When $\nu(\{0\})>0$, $E(D)$ is infinite. Since all elements in $\{N_{k}(t)\}_{k\geq 1}$ are independent and follow Poisson distribution, variable $N_{+}(t)$ is Poisson distributed, and so do $D$ and $N_{0}(t)$ when their expected values are finite.

Write $S(t)=\sum_{k=1}^{\infty}kN_{k}(t)$ which is the number of individuals observed before time $t$.

From (3.1), for $k\geq 0$,

By the thinning property of Poisson processes,

is the intensity measure of the rate for the species represented $k$ times in $[0,t_{0}]$. We can interpret $(k+1)E(N_{k+1}(t_{0}))/t_{0}$ as the expected total rate for all species represented $k$ times in $[0,t_{0}]$. As pointed out in Section 2 and from (3.4), $\Lambda=E(S(t_{0}))/t_{0}$ is the expected total rate for all species. Conditional on the event that we will observe an individual at a given future time $t_{1}>t_{0}$, the probability that this individual belongs to a species represented $k$ times in $[0,t_{0}]$ is equal to

(formal proof is given in Appendix D) which corresponds to the renowned Good-Turing frequency estimator in Good [21]. It is worth pointing out that the above equation holds for individual observed at a given future time rather than the next observed individual.

From the well-known relation between independent Poisson random variables and multinomial distribution, model (1.1) is valid under the framework with

Formal proof is given in Appendix E. If $E(D)$ is finite, $\lim_{t\to\infty}p_{k}(t)=0$ for any fixed $k$. The joint probability mass function of $\tilde{N}(t_{0})$ is

In terms of the expected FoF, the log-likelihood function is

In terms of $p(t_{0})$ and $E(N_{+}(t_{0}))$, it is

If the unknown vector $p(t_{0})$ and the quantity $E(N_{+}(t_{0}))$ are unrelated, the above log-likelihood function implies that the maximum likelihood estimator (MLE) of $p(t_{0})$ is the conditional maximum likelihood estimator (conditional on the observed $n_{+}(t_{0})$) for the multinomial distribution in (1.1). The MLE of $E(N_{+}(t_{0}))$ is $n_{+}(t_{0})$.

Denote the expected (empirical) SAC (ESAC) as $\psi(t)=E(N_{+}(t))$. Condition (2.1) guarantees that $\psi(t)$ is finite for any finite $t$. For a real-valued function $g(t)$, let $g^{(m)}(t)$ stand for the $m$-order derivative of function $g(t)$. Clearly $\psi(0)=0$. From (3.2),

Note that $\psi^{(1)}(0)=\int\nu(d\lambda)=\Lambda$. From (3.1) and (4.1),

Analogous expression for (4.2) appears in Béguinot [2] as an approximate formula under the multinomial model for fixed total number of observed individuals for the species-sample-size curve with the derivative operator replaced by the difference operator.

Before studying the link between ESAC and SAD, two mathematical terms are needed. A function $g(t)$ is a Bernstein function if it is a nonnegative real-valued function on $[0,\infty)$ such that $(-1)^{k+1}g^{(k)}(t)\geq 0$ for all positive integer $k$ [44]. An infinitely differentiable function $f(s)$ on an interval $A$ is called absolutely monotone in $A$ if $f^{(k)}(s)\geq 0$ for $k=0,1,\ldots$ and $s\in A$. The relation between absolutely monotone function and probability generating function is well known (see for example Strook [45]).

a.

A function $\psi(t)$ is the ESAC of an MPPP if and only if $\psi(t)$ is a Bernstein function that passes through the origin.

b.

A function $h_{t}(s)$ is the probability generating function of $p(t)$ for a fixed positive $t$ of an MPPP if and only if $h_{t}(0)=0,h_{t}(1)=1$, and $h_{t}(s)$ is absolutely monotone in $(-\infty,1)$.

c.

An MPPP is uniquely determined by its $\psi(t)$, or $(h_{t_{0}}(s),\psi(t_{0}))$ for the probability generating function $h_{t_{0}}(s)$ of $p(t_{0})$ for a fixed positive $t_{0}$.

d.

For any MPPP, and $t>0$, we have

$$h_{t}(s)=1-\frac{\psi((1-s)t)}{\psi(t)},~{}~{}~{}(s\in(-\infty,1]).$$

(4.3)

Proof of Proposition 1 is given in Appendix F. Equation (4.3) is the population version of (1.4).

It worths pointing out that Boneh et al. [4] proved that $(-1)^{k+1}\psi^{(k)}(t)\geq 0$ for a different setting ($D$ independent Poisson processes with different rates) and called it “infinite order alternating copositivity”.

From (3.6), the log-likelihood function can be re-expressed as a function of $\psi(t)$.

It can be shown that the Taylor expansion of $\psi(t)$ at $t_{0}$ converges to $\psi(t)$ when $0\leq t<2t_{0}$:

It signifies that Good-Toulmin estimator in (1.3) performs satisfactorily in short-term extrapolation when $t_{0}<t<2t_{0}$. We deduce from (4.5) the following unbiased estimator for different order of derivative of $\psi(t)$

We use (4.6) only for $0\leq t\leq t_{0}$ ¹¹1We include $t=0$ here for completeness. Remember that $\psi^{(1)}(0)=\Lambda$ can be infinite. because outside this interval, $\hat{\psi}^{(j)}(t)$ may not have the correct sign $(-1)^{j+1}$.

Estimator in (4.6) is useful. For example, a concave downward curve when we plot $1/\hat{\psi}^{(1)}(t)$ for $t\in(0,t_{0}]$ ²²2It corresponds to the diagnostic check for the log-series distribution in Table 1 in Section 5. is an indication that $E(D)=\infty$ because if we believe that there is a linear function $b+ct$ with positive $b$ and $c$ such that $b+ct\geq 1/\psi^{(1)}(t)$ for all $t\geq 0$, then

From (4.2), a parallel result of Equation (4.6) is the following relation among expected FoFs

which can be proved using the law of iterated expectations when $0\leq t\leq t_{0}$. Furthermore, for $0\leq t\leq t_{0}$

MPPP is nonparametric in nature. It is defined by an intensity measure, an ESAC or an SAD. When parametric approach is preferred, we restrict our interest to a family of distributions in MPPP, say by putting constraints on the SAD. A way to portray an SAD is to delineate its probability ratio, $p_{j}(t)/p_{j+1}(t)$ for $j=1,2,\ldots$. It is equivalent to examine $-\psi^{(j)}(t)/\psi^{(j+1)}(t)=tp_{j}(t)/[(j+1)p_{j+1}(t)]$, which we call the $j$th derivative ratio. From (4.1) and the Cauchy-Schwarz inequality, for $j=1,2,\ldots$ and $t\geq 0$, $\psi^{(j)}(t)\psi^{(j+2)}(t)\geq(\psi^{(j+1)}(t))^{2}$. It deduces that $j$th derivative ratio is always a nonnegative nondecreasing function of $t$. Among all derivative ratios, the first derivative ratio is most important because $-\psi^{(1)}(t)/\psi^{(2)}(t)=-[d\log(\psi^{(1)}(t))/dt]^{-1}$ which relates to the logarithmic derivative of $\psi^{(1)}(t)=E(N_{1}(t))/t$, the expected total rate for unseen species at time $t$. The following proposition gives a sufficient condition for it.

The simplest nontrivial Bernstein function is the positive linear function on $[0,\infty)$. It suggests the following fundamental family of distributions in MPPP.

Family $\mathrm{LDR}_{1}$ encompasses some common SADs and ESACs. We list the characteristics of all models in $\mathrm{LDR}_{1}$ in Table 1. When $b=0$, $c$ must be larger than 1 (see condition (ii) in Proposition 2), otherwise from the second row of Table 1, $\lim_{b\downarrow 0}\psi(t)=\infty$ for finite $t$. In Table 1, equality of transformed $\psi^{(1)}(t)$ for some models is presented. Such equalities can be used to produce diagnostic check for specified SAD in $\mathrm{LDR}_{1}$ when $\psi^{(1)}(t)$ is replaced by its estimator in (4.6).