A Model for Species Accumulation based on the Multinomial Distribution

Miguel Nakamura Centro de Investigación en Matemáticas Apartado Postal 402 Guanajuato, Gto. 36000, Mexico

Felipe Peraza Escuela de Ciencias Físico-Matemáticas Universidad Autónoma de Sinaloa Apartado Postal 1872 Ciudad Universitaria, Culiacán, Sin. 80030, Mexico

Abstract

A species accumulation curve is a plot of the number of species found in a given region as a function of some measure of effort used to find them. In ecological studies, there are important goals achieved when the relationship between these two variables is understood. In this paper, we construct a model—which turns out to involve a multinomial distribution—based on the premise that the form of the accumulation curve depends on the joint behavior of species’ capture probabilities. Maximum likelihood and profile likelihood methods are explored for the estimation of parameters. The expectation-maximization algorithm is described for dealing with incomplete collecting lists. Two examples which illustrate the techniques developed are considered, and comparisons are made with estimates based on other approaches.

Key words: Number of species; EM algorithm; profile likelihood.

Running head: Species accumulation curves.

1 Introduction

Although the problem addressed in this paper finds its motivation in biodiversity studies, it is possible that the setting may be applicable in other sciences as well. In a general context, suppose a population of objects is classified into an unknown number of classes and that one can observe and classify a sequence of individual objects drawn at random. Interest lies in modelling the relationship between the total number of classes identified and the effort applied. In ecological studies, there are important goals achieved when a theoretical relationship is found between these two variables. Examples are found in the planning of future expeditions, quantitative comparisons between two or more collection efforts, and in biodiversity and conservation studies when interest lies in total number of species present.

In this ecological setting, a species accumulation curve or “collector’s curve” is a plot of the number of species found in a given region, S(n), as a function of some measure n of the effort used to find them. The simplest measure of effort is the total number of individuals examined, but in ecological studies typical units of measure are quadrant areas, sampled biomass, hours spent on observation, number of traps per day, and others (see Coldwell & Coddington, 1994).

Numerous references exist in the literature regarding estimation of the number of different types of objects in a population (see Bunge & Fitzpatrick, 1993). For example, biologists and ecologists are interested in estimating the total number of species in a population of plants or animals (Engen, 1976; Soberón & Llorente, 1993; Coldwell & Coddington, 1994). A linguist may be interested in estimating the size of the vocabulary of an author (Efron & Thisted, 1976), or a programmer in the number of bugs in new software, among other examples. But most of these works are either based on a capture-recapture sampling scheme or on n_{x}, the number of objects which appear exactly x times, and x_{s}, the number of individuals which belong to class s. Computing these stated frequencies requires the record of the whole collection process, which is not necessarily available in ecological studies. Furthermore, the techniques do not always take into account the time at which each new identified class appears or is discovered, so an explicit relationship between effort and capture is not addressed.

Soberón & Llorente (1993) (henceforth denoted SL) relate the concepts of effort and number of accumulated species by analytical methods. Their model is based on pure birth processes, where parametric transition probabilities (called collecting functions) are postulated on the basis of biological interpretation. The expected value of the ensuing stochastic process, A(n), is fit to the observed accumulation curve, using regression techniques to estimate the parameters of the process. Because comparisons will be made in a later section, we here state the three mean functions which SL consider. Their exponential model is A(n) = (a/b) [1 − exp(−bn)]. The Clench model (Clench, 1979), which is based on the Michelis-Menten equation, A (n) = an/ (1 + bn) is also a special case of the SL approach. Here A(n) → a/b as n → ∞, so a/b is interpreted as the total number of species present. The logarithmic model is A (n) = ln (1 + zan) / z. In this case, A (n) diverges as n → ∞.

Coldwell & Coddington (1994) discuss objectives in the modelling of species accumulation curves, and state that one of the main goals is often to predict the expected increase in species richness per additional unit of effort. The shape of the accumulation depends on the relative abundance of the sampled species. Some authors have adopted this fact implicitly (see Engen, 1978, for example). Although the notion of accumulation curve is not used (because attention is focused on estimating total number of species), the methods are based on abundance specified through explicit forms or probability distributions.

In this paper, we first propose a model using a different approach. We will postulate that the form of the accumulation curve depends on the capture probabilities of each species (that is, the probability that a given species is observed on a given trial). These probabilities depend, amongst other things, on the relative abundance present and the sampling method used. Additionally, the model will take into account the discrete nature of units of effort, the possibility of multiple occurrences of new species at a single unit of time, and the fact that capture may depend on random factors.

A few working hypotheses will be required in order to obtain analytical results. We will assume that the collection processes are independent in time; that each species is captured independently from all others, each with a certain probability; and that capturabilities are a realization of a random variable. This idea of random capturabilities is not new. However, no distribution is here justifiable on empirical means because no observed frequencies are available. We resort to a {Beta(α, β), α > 0, β > 0} family on grounds of richness and because it yields analytical tractability. The beta family has previously been used by some authors (Seber 1982; Otis et al. 1978) for modelling capture probabilities.

In Section 2 we establish that the vector of new species captured in the first n times has a multinomial distribution, with parameters which are functions of (k, α, β), where k is the total number of species in the area and α, β are the parameters of the beta distribution. Maximum likelihood and profile likelihood methods are explored for the estimation of parameters.

An additional issue considered here is when collecting lists are incomplete. In many cases, a complete record of the collection process is not available. This may be considered a case of censoring, since there exist time lapses where observation was performed but accumulation detail has been lost, due for instance to a compounding effect (historical reporting of lists where units are months instead of days, for example). For estimation in these cases, we will use the EM algorithm (Dempster et al., 1977).

To finalize, in Section 4 we analyze two collecting lists which illustrate the techniques developed. Results will be discussed and comparisons will be made with estimates based on other approaches.

2 Model formulation

2.1 Basic model

Let k denote the total number of different species present in an area. Its value is fixed, but unknown. The collection process consists in a sequence of trials indexed by n, n = 1, 2, …, which may conveniently be referred to as a time index and represents the amount of effort . The probability of observing a species in unit time is a function of various factors, such as abundance, diversity, detectability, field worker’s experience, and the method of observation. For modelling purposes, the totality of these factors may be considered as random, due to weather conditions or specific characteristics of the particular land area selected for sampling. For this reason, it is reasonable to consider that the probabilities that each species is caught or observed in a time unit are a realization of a random vector P_{∼} = (P_{1}, P_{2}, …, P_{k}) , having density supported on [ 0, 1]^{k}.

For i = 1, 2, ..., k, in any time unit let

p_{i} = Probability of capturing the species labeled i,

where (p_{1}, …, p_{k}) is a single realization of P_{∼}.

We will adopt the following assumptions as a starting point:

(a) Collection processes are independent in time.
(b) At each time n, each species is observed independently
from all others, with the same probabilities p_{1}, p_{2}, …, p_{k}.
(c) P_{1}, P_{2}, …, P_{k} are independent and identically distributed.

(1)

Assumption (1b) is the most critical, because it indicates that capture probabilities remain constant in time. It is well-known that abundance in many species is related directly to periodicity in yearly weather, so this assumption is definitely not applicable in all cases. It will be valid to the extent in which the sampling period is not affected by the seasons or by the biological age of specimens. The role of assumption (1c) is to model capture probabilities as a group. We are not interested in the probability of capturing the i-th species of the population, p_{i}, but in their joint behavior, since it is this which determines the shape of the accumulation curve. Capture probabilities may be related to the notion of abundance (see Engen, 1978), through a mechanism which can be as simple as proportionality, or a far more complex process. The point here is that we are interested only in the final capture probabilities, disregarding any relationships which may be present with abundance. The Beta(α, β) distribution will be considered here because of its flexibility and analytical convenience. Seber (1982) and Otis et al. (1978) also mention the use of the Beta(α, β) distribution for modelling capture probabilities.

The fact of finding or not finding the i-th species at time n, will be described through a Bernoulli random variable. We use the notation Y_{i}^{(n)} for such random variable, that is

Y_{i}^{(n)} ∼ Ber(p_{i}), for n = 1, 2, ...; i = 1, 2, ..., k.

We also use the notation

X_{n} = number of new species collected at time n.

Variables X_{n} are then described in terms of Y_{i}^{(j)} as follows:

X_{0} = 0,  (2)
X_{1} = ∑_{i=1}^{k} Y_{i}^{(1)},
⋮
X_{n} = ∑_{i=1}^{k} Y_{i}^{(n)}(1 − Y_{i}^{(n−1)}) ⋯ (1 − Y_{i}^{(1)}).

In these terms, the accumulation curve is the plot of S(n) vs. n, where S(0) = 0, and S(n) = S(n − 1) + X_{n}, for n = 1, 2, . . .. The random variable X_{n} is the size of the jump in the collection process at time n.

Under assumptions (1), the joint distribution of (X_{1}, X_{2}, …, X_{n}) is multinomial,

MN_{n} (k, E(P_{1}), E (P_{1} (1 − P_{1})) , …, E(P_{1} (1 − P_{1})^{n−1})) .

That is, its probability density is

Pr (X_{1} = x_{1}, X_{2} = x_{2}, …, X_{n} = x_{n})
= \frac{k!}{x_{1}! x_{2}! ⋯ x_{n}!(k − ∑_{i=1}^{n} x_{i})!} [E(1 − P_{1})^{n}]^{k−∑_{i=1}^{n} x_{i}} ∏_{i=1}^{n} E^{x_{i}} [P_{1}(1 − P_{1})^{i−1}] , (3)

where 0 ≤ x_{i} ≤ k for all i and satisfy 0 ≤ ∑_{i=1}^{n} x_{i} ≤ k. The proof is sketched in Appendix A. It is interesting to note that the collection process described is equivalent to a multinomial experiment. This interpretation considers each species as performing a sequence of k independent trials, each with n+1 mutually exclusive outcomes A_{0}, A_{1}, …, A_{n}—where A_{i} is the event “species is observed for the first time at time i” for i ≥ 1 and A_{0} is the event “species is not yet observed in n time units”—such that the probability of A_{i} occurring in each independent trial is E [P_{1} (1 − P_{1})^{i−1}] for i ≥ 1.

Assuming that the capture probabilities are a realization of a random vector having independent and identically distributed Beta(α, β) entries —a situation which will hereafter be referred to as the “beta case”—, the joint density of (X_{1}, X_{2}, …, X_{n}) is reduced to

Pr (X_{1} = x_{1}, X_{2} = x_{2}, …, X_{n} = x_{n})  (4)
= \frac{k!}{x_{1}! x_{2}! ⋯ x_{n}!(k − ∑_{i=1}^{n} x_{i})!} [\frac{b (α, β + n)}{b (α, β)}]^{k−∑_{i=1}^{n} x_{i}} ∏_{i=1}^{n} [\frac{b (α + 1, β + i − 1)}{b (α, β)}]^{x_{i}} ,

where b (·, ·) is the beta function. This density function is the basis for inferential procedures which will be addressed in Section 3. A few useful consequences of these results in the beta case follow.

For fixed n, the value of the accumulation curve S(n) has a binomial distribution,

S(n) ∼ Bin (k, 1 − \frac{b (α, β + n)}{b (α, β)}) ,

and thus

E (S(n)) = k (1 − \frac{b (α, β + n)}{b (α, β)}) , and

V (S(n)) = k (\frac{b (α, β + n)}{b (α, β)}) (1 − \frac{b (α, β + n)}{b (α, β)}) .

(5)

The value of E (S(n)) as a function of n is comparable in meaning to the values of A(n) which are produced by the SL approach.

Let s be the number of species observed until time N. If Z_{M} = ∑_{i=N+1}^{N+M} X_{i} is the number of species which would be collected in time units N + 1 to N + M, then

Z_{M} | N_{1}, …, N_{N} ∼ Bin (k − s, 1 − \frac{b(α, β + N + M − 1)}{b(α, β + N)}) .

A plot of

E [Z_{M} | N_{1}, …, N_{N}] = (k − s) (1 − \frac{b(α, β + N + M − 1)}{b(α, β + N)})  (6)

versus M would be a complementary device useful in deciding if a given value of M is rewarding in terms of collection effort.

3 Inference

We now establish inference procedures for the stated model. Data is X_{1} = x_{1}, X_{2} = x_{2}, …, X_{n} = x_{n}, where the vector (X_{1}, X_{2}, …, X_{n}) has density given by (4). Since this density is given explicitly, we will consider estimation via maximum likelihood, and because interest lies mainly in the parameter k, we also turn to the notion of profile likelihood. When censored data is present, the EM algorithm is a plausible way of obtaining parameter estimates.

All of these methods require optimization of non-linear functions in dimension three, at most. In our experience, the method of conjugate gradients (Press et al., 1988) was adequate, although careful grid selections had to be considered in searching for starting points. Derivatives of the log-likelihood are useful for this purpose and are stated in Appendix B.

3.1 Likelihood and profile likelihood

Let x_{1}, x_{2}, …, x_{n} be a realization of the random variables X_{1}, X_{2}, …, X_{n} whose density is given by (4), and let s = ∑ x_{i}. The log-likelihood (disregarding a constant additive term) is given by

l (k, α, β) = log (k!) − log ((k − s)!) − (k − s) log [\frac{b (α, β + n)}{b (α, β)}]
+ ∑_{i=1}^{n} x_{i} log [\frac{b (α + 1, β + i − 1)}{b (α, β)}] .  (7)

Maximization of this quantity provides point estimates of parameter values and functions of them. In particular, the quantities described towards the end of Section 2.1, which are of special interest in the biological context, may be estimated (by a substitution principle). However, precision and confidence intervals are not as straightforward to obtain, because the problem is highly non-regular. In Monte Carlo experimentation, we have found conclusive evidence that the non-regularity in this setting is such that normal approximations are not really useful for inferential procedures, and bootstrap intervals are preferred instead.

Since α and β may be considered nuisance parameters, and k the parameter of interest, profile likelihood seems a plausible tool to use for inference (Hinkley et al., 1991; Barndorff-Nielsen & Cox, 1994). We here only state material necessary for constructing confidence intervals for k. The profile likelihood is

L_{P} (k) = L (k, α̂_{k}, β̂_{k}) = sup_{α, β} L(k, α, β),

where L(k, α, β) is the likelihood and α̂_{k}, β̂_{k} are the maximum likelihood estimators of α and β for fixed k. Inference on k is based on treating L_{P} (k) as if it were a legitimate likelihood. The log-profile likelihood is

l_{P} (k) = sup_{α, β} l(k, α, β)

where l(k, α, β) is defined by (7). An approximate (1 − a) × 100% confidence interval for k based on profile likelihood is

R_{P} = \{k : −2 log \frac{L_{P} (k, α̂_{k}, β̂_{k})}{L_{P} (k, α̂k̂, β̂k̂)} ≤ χ_{1,1−a}^{2}\} ,  (8)

(see Barndorff-Nielsen & Cox, 1994). Correct coverage probabilities rely on the fact that the profile likelihood ratio statistic,

r_{P} (k) = −2 log \frac{L_{P} (k, α̂_{k}, β̂_{k})}{L_{P} (k, α̂k̂, β̂k̂)}

has an asymptotic χ_{1}^{2} distribution. We have performed a few simulations (reported in Section 4.1) to examine the distribution of r_{P} (k) . The conclusion is that the distribution of r_{P} (k) is very similar to a χ_{1}^{2} distribution, and more so when α, β are large for fixed k and for large k for fixed (α, β). We believe these results suggest that the χ_{1}^{2} approximation is reasonably good for practical applications. Future research may consider the possibility of modifying the profile likelihood, as in the conditioning argument of Cox & Reid (1987), the modification of Barndorff-Nielsen (1986), or the adjustment of McCullagh & Tibshirani (1990).

3.2 Bootstrapping methods

An alternative technique for estimating the covariance matrix and other measures of precision of estimators is the bootstrap (Efron, 1982). Here, we present the implementation of the method for building confidence intervals for functions of the parameters in (4).

Let θ = (k, α, β) and let F (θ) denote the multinomial distribution (4); the data vector y = (x_{1}, x_{2}, …, x_{n}) is a sample from F (θ) . Suppose interest is in g(θ) for some real-valued function g. Based on y, one computes θ̂, the maximum likelihood estimator, and interest lies in evaluating the precision of g(θ̂) as an estimator of the true value of g(θ). Computation of θ̂ is achieved numerically, and bootstrap intervals are found as follows (Efron & Tibshirani, 1986):

(i) Find θ̂ using data y;

(ii) Simulate B independent vectors y^{∗} (1) , y^{∗} (2) , …, y^{∗} (B) from the multinomial distribution F \binom{}{θ̂}, and for each of them compute g(θ̂^{∗} (b)), the maximum likelihood estimator of g(θ) based on y^{∗} (b) , b = 1, 2, …, B.

(iii) Having produced g(θ̂^{∗} (1)), g(θ̂^{∗} (2)), …, g(θ̂^{∗} (B)), define Ĝ, the empirical distribution function of {g(θ̂^{∗}(b)), b = 1, 2, 3, …, B}, as

Ĝ (s) = \frac{# (g(θ̂^{∗}(b)) ≤ s)}{B}.

This last empirical distribution is regarded as an approximation to Pr{g(θ̂) ≤ s}, and motivates the following approximate (1 − 2a) × 100% confidence interval for g(θ) (Efron, 1982, Sec. 10.4):

[ Ĝ^{−1} (a) , Ĝ^{−1} (1 − a)] ,

where a, 0 ≤ a ≤ 0.5, is fixed.

Under model (4), the simulation in step (ii) above consists in generating p_{1}^{∗}, p_{2}^{∗}, …, p∗ ; k̂

as a random sample from a Beta (α̂, β̂) distribution, and then simulating the process (2). Equivalently, B independent draws of vectors y^{∗} = (X_{1}^{∗}, X_{2}^{∗}, …, X_{n}^{∗}) may be made from the

MN_{n} (k̂, \frac{b(α̂ + 1, β̂)}{b(α̂, β̂)}, \frac{b(α̂ + 1, β̂ + 1)}{b(α̂, β̂)}, …, \frac{b(α̂ + 1, β̂ + n − 1)}{b(α̂, β̂)})

distribution.

3.3 Censored observations: EM algorithm

When censored or missing values occur, likelihoods may be maximized using the EM algorithm, formalized by Dempster, Laird, & Rubin (1977). Using their notation, let x = (y, z) denote observed data y, augmented with missing data z. Let f (x; θ) and g (y; θ) be families of densities indexed by the parameter θ. The relationship between f and g is given by g (y; θ) = ∫_{X (y)} f (x; θ) dx, where X (y) is the set of all possible values for missing data given a fixed value of y. The aim of the EM algorithm is to maximize g (y; θ) as a function of θ, but using f (x; θ) through a different type of relationship.

Define the function

Q (θ′, θ) = E [log f (x; θ′) | y, θ] .

The iteration θ^{(p)} → θ^{(p+1)}, which is called an EM step, is defined as follows:

E Step: Compute Q(θ, θ^{(p)}).

M Step: Choose θ^{(p+1)} as the value of θ which maximizes Q(θ, θ^{(p)}).

The heuristic idea behind this is that we would like to choose θ^{∗} to maximize log f (x; θ). Since we cannot compute log f (x; θ) because we do not observe all of x, we instead maximize its expected value given y and the current value of θ. Iteration produces convergence to the maximum likelihood estimate.

3.3.1 Implementation for the beta case

Let (m_{1}, s(m_{1})), (m_{2}, s (m_{2})) , …, (m_{r}, s (m_{r})) be the observed accumulation curve. The iteration prescribed by the EM algorithm is described as follows. A brief justification of these steps is included in Appendix C.

E Step: First define the quantities x_{1}, …, x_{n} as follows.

For i = 1, 2, … r, do:

If m_{i} = m_{i−1} + 1, set

x_{m_{i}} = s(m_{i}) − s(m_{i−1}).

Otherwise, for j = m_{i−1} + 1, m_{i−1} + 2, …, m_{i}, set

x_{j} = [s (m_{i}) − s (m_{i−1})] \frac{e_{j}}{∑_{m_{i−1}<j≤m_{i}} e_{j}} ,

where

e_{j} = \frac{b (α^{p} + 1, β^{p} + j − 1)}{b (α^{p}, β^{p})}.

Then define

Q^{∗} ((k, α, β) , (k^{p}, α^{p}, β^{p}))

= log (k!) − log ((k − s (m_{r}))!) − (k − s (m_{r})) log [\frac{b (α, β + m_{r})}{b (α, β)}]

+ ∑_{i=1}^{m_{r}} x_{i} log [\frac{b (α + 1, β + i − 1)}{b (α, β)}] .

M Step: Choose (k^{p+1}, α^{p+1}, β^{p+1}) as the value of (k, α, β) which maximizes

Q^{∗} ((k, α, β) , (k^{p}, α^{p}, β^{p})) .

4 Simulations and examples

This final section is intended to show distributional properties of the profile likelihood ratio statistic via simulation, and two examples previously studied by SL. The first example regards a non-censored collecting list of bats, and the second is a censored list of butterfly data which will therefore be analyzed using the EM algorithm.

4.1 Distribution of r_{P} and confidence intervals for k

This Monte-Carlo experiment shows the distribution of r_{P} and the empirical coverage probability of confidence intervals based on this statistic. Consider the eight cases obtained from combinations of n = 100, 150; k = 100, 150; and (α, β) = (.5, .8) , (5, 250) . Each of 1000 replications consists in simulating an accumulation curve, computing r_{P} , and finding the confidence region R_{P} defined by (8) with a = 0.05. The empirical coverage probability is then the proportion of times R_{P} contains k.

Histograms for values of r_{P} prove to be quite close to superimposed plots of χ_{1}^{2} densities. The empirical distributions are fairly close to a χ_{1}^{2} for all configurations, with Beta (5, 250) giving the best agreements. There is a slight tendency to heavier tails than expected, but this effect is not serious. Table 1 sheds information in terms of sample means. The mean of r_{P} is closer to 1 (the mean of the χ_{1}^{2} distribution) for (α, β) = (5, 250) for the same value of k. For these same cases, coverage is closer to 95%, which is the intended level.

Table 1: Summary results of 1000 simulations of the profile likelihood ratio statistic under eight parameter configurations.

Parameters  n mean Coverage Probability
α = 0.5 β = 0.8 100 1.2429  91.4
k = 100 150 1.2491 91.9
α = 0.5 β = 0.8 100 1.1348  93.6
k = 150 150 1.1805 93.4
α = 5 β = 250 100 1.0771  94.7
k = 100 150 1.0074 94.6
α = 5 β = 250 100 1.085  93.6
k = 150  150 1.1387  94.6

4.2 Examples

4.2.1 Chajul bat data

This example considers the list of bat species at Chajul Biological Station in southeastern Mexico. The list was analyzed previously by SL; these authors report 50 observed species after setting nets on 49 nights. The unit of effort in this example is sampled night.

We use maximum likelihood for parameter estimation. Estimated total number of species in the area is 54, while the estimated values of the beta parameters describing capture probabilities are α̂ = 24.00 and β̂ = 441.49. Figure 1 shows the Beta(α̂, β̂) distribution. It is almost symmetric and concentrated about 0.06, which suggests that all bat species roughly share a common capture probability and are fairly prone to be captured. This is, there do not exist species which are either extremely rare or extremely hard to capture.

Figure 2 shows the profile likelihood function, l_{p}(k), and the 95% confidence interval for the total number of species, k. The 95% interval is (50, 68]. Because of empirical evidence and since estimates are suggesting large values of α and β, the actual coverage probability of this interval is prone to be near 95%.

It is interesting to compare these results with those found by SL. Their “Clench” and “exponential” models produce estimates 85 and 56 for k, respectively. Using biological arguments, they favor the exponential model, so the estimate of k which we find using profile likelihood agrees with the previous analysis, with the added feature of confidence interval being stated.

Figure 3 shows the estimated expected value (5) and the conditional expectation (6) for M ≥ 50, using maximum likelihood parameter estimates. In this example, both of these are very close to each other, the reason being that the estimated value of the accumulation curve at time 49 is very close to the actual value observed. The following statements refer to estimated expected value and estimated conditional expected value (in that order), and show the type of answers which can be obtained: at time 79, estimated values are 52.86 and 53.00, that is, one expects to add 3 new species to the list after 30 additional units of effort.

At time 99, estimated values are 53.53 and 53.59, that is, 50 additional units are required to achieve (in the sense of an average) the estimated total number of species present.

4.2.2 Pakitza butterfly data

SL analyzed butterfly data obtained in a set of 200 person-hours (in September 1989) at Pakitza Biological Station, Peru. Considering the Clench, exponential and logarithmic models (see Section 1), produces indistinguishable fits over the range of observed data, but extrapolations are strikingly different (see Figure 4). SL argue in favor of the logarithmic fit because of the biological complexity and the large extent of the area.

In this example, the complete history of the collection process is not available. There are 14 points in time (person-hours) where the height of the collecting curve is known. This data set provides an example of a censored list, in the sense of Section 3.3. The parameter estimates using the EM algorithm (Section 3.3) in the multinomial model are k̂ = 981, α̂ = .68 and β̂ = 74.99. The estimated total number of species is thus 981, compared with the asymptotic values of 890 and 611 obtained in the Clench and exponential models, respectively.

The value 981 would seem, at first sight, to be very extreme when compared to either of the estimates obtained by the regression approach. In this respect, there exists an interesting posterior observation: after continuation for 565 person-hours, the species accumulation count stood at 979. This value is quite close to the estimate obtained by our model, and exceeds the values obtained by the Clench and exponential models. The logarithmic model lacks an asymptote, which makes it hard to interpret over the extrapolation range. Figure 4 shows the three models fitted by SL, as well as the estimated expected accumulation curve in our model computed with estimated parameter values.

In Figure 5 we plot the Beta(α̂, β̂) density function. The shape of this distribution, in contrast with the previous example regarding bats, suggests that the Pakitza ecosystem has more complex characteristics, that is, that it is a region with much greater biodiversity and geographical extension. Here, many species are hard to be sighted, and few species have a high probability of being observed. This produces an accumulation curve which grows quite rapidly in the beginning, but later requires very large amounts of effort to achieve further significant increase.

Figure 6 shows plots of estimated expected accumulation and conditional expectation. In this situation, these functions are notoriously different.

We observe that the estimated expected value of accumulation is 574.86 at the next time unit, while the estimated conditional is 568.82, that is, 0.82 species on the next sighting attempt. If the collection process were continued for 150 additional time units, we find expectation at 718.71, while the conditional expectation is 724.12 species. At time 565, at which future sampling reports 979 species, the model yields 751.97 as an estimate of the mean accumulation curve. This value seems too low. Although we do not have further information regarding the subsequent study, a possible reason for this is that the original study used a different sampling procedure. This suggests one possible extension of the model, of immediate application: that capture probabilities vary in time (in this case, it is possible that capture probabilities changed because of a change in collection crews and/or methods).

5 Conclusion

The model presented in this paper constitutes an alternative for solving the problem of modelling the relationship between effort and number of accumulated species. Taking into account the discrete nature of units of effort and considering that only capture probabilities are directly relevant, make the model more realistic and easily interpreted. The availability of a formula for conditional expected values of the collection process given the observed curve, is an advantage. This seems an adequate concept, in addition to the fact that the relationship between effort and capture is modelled.

On the other hand, it is important to evaluate precision of inferential procedures. Quantities such as estimators of conditional expectations, which are usually not considered, are easy to compute under this approach. Because censored collecting curves often appear, the solution in this case is also an important advantage. Results in a number of real and simulated data sets are very encouraging.

A disadvantage is that the bootstrap and EM methods here described require fairly sophisticated numerical methods which are not readily available in commercial software. The approach of SL, in contrast, uses simple least-squares procedures which are commonly available.

One interesting problem which remains open is the exploration of corrections to the profile likelihood procedures used for inference on the parameter k. However, at a modelling stage, a more compelling issue is the generalization to capture probabilities which are more general. The model we have stated in this paper is applicable in the event that capture probabilities remain approximately constant over the span of the collection, and in a large number of applications, this situation is not necessarily true.

Correlations between entries in the vector P give rise to intractable expressions, but perhaps the main reason for not considering them here is that postulating a correlation function seems to be a difficult task. On the other hand, it is certainly natural to conceive that capture probabilities may vary in time. Possible reasons for this are seasonal variations in abundance and diversity due to weather and immigration, or the field worker gathering experience or changing a method of capture. A possible generalization of the present model is to consider that probabilities are realizations of random variables with time-dependent densities. This suggests to model non-stationary capture distributions through a family indexed by time, perhaps with the inclusion of additional parameters. A brief argument would show that the resulting distribution is again multinomial.

Acknowledgments

Miguel Nakamura’s work was partially supported by CONACYT grant number 4224E9405. David Sprott suggested the shorter and more transparent proof of the main multinomial result. We also thank Jorge Soberón for relevant discussion at the modeling and interpretation stages.

A Proof of main result

Let Z_{ij} = 1 if species i is observed for the first time on trial j, 1 ≤ j ≤ n, and Z_{ij} = 0 otherwise. Let Z_{i0} = 1 if species i is not yet caught by the n-th trial, 1 ≤ i ≤ k and 0 otherwise. For 1 ≤ i ≤ k and 1 ≤ j ≤ n we define the conditional probabilities

π_{ij} = Pr(Z_{ij} = 1 | P_{i} = p_{i}) = p_{i}(1 − p_{i})^{j−1}, and
π_{i0} = Pr(Z_{i0} = 1 | P_{i} = p_{i}) = (1 − p_{i})^{n}.

Let P have density g. The marginal probabilities are Pr(Z_{ij} = 1) = ∫_{p=0}^{1} π_{ij} g(p) dp = E(π_{ij}), for 1 ≤ i ≤ k and 0 ≤ j ≤ n. Because for fixed i only one of the variables Z_{i0}, Z_{i1}, …, Z_{in} takes on the value 1 while the rest are all 0, it can be easily checked that the joint probability generating function (pgf) of this set of random variables is

φ_{Z_{ij}}(s_{0}, …, s_{n}) = E(π_{i0})s_{0} + ⋯ + E(π_{in})s_{n}.

The numbers of new species caught on each trial are X_{j} = ∑_{i=1}^{k} Z_{ij}, 1 ≤ j ≤ n. Assuming independence of species, the joint pgf of these variables is

φ_{X_{i}}(s_{0}, …, s_{n}) = φ_{Z_{ij}}^{k}(s_{0}, …, s_{n}) = [E(π_{i0})s_{0} + ⋯ + E(π_{in})s_{n}]^{k},

which is identified as the pgf of the multinomial distribution with index k and n + 1 classes having probability parameters E(π_{ij}). When each P_{i} has the same distribution, these are E(π_{ij}) = E[P (1 − P)^{j−1}], for j ≥ 1 and E(π_{i0}) = E[(1 − p)^{n}].

B Expressions for derivatives

The gradient of the log-likelihood (7) is

\frac{∂l (k, α, β)}{∂k} = ∑_{i=0}^{s−1} \frac{1}{k − i} + ln Γ (α + β) + ln Γ (β + n) − ln Γ (β) − ln Γ (α + β + n) ,
\frac{∂l (k, α, β)}{∂α} = − (k − s) ∑_{j=0}^{n−1} \frac{1}{j + (α + β)} − ∑_{i=1}^{n} x_{i} ∑_{j=0}^{i−1} \frac{1}{j + (α + β)} + \frac{s}{α},

\frac{∂l (k, α, β)}{∂β} = (k − s) ∑_{j=0}^{n} (\frac{1}{j + β} − \frac{1}{j + (α + β)}) + ∑_{i=2}^{n} x_{i} ∑_{j=0}^{i−2} \frac{1}{j + β}

− ∑_{i=1}^{n} x_{i} ∑_{j=0}^{i−1} \frac{1}{j + (α + β)}.

C Steps in the EM algorithm

The E step is based on a decomposition of the multinomial distribution. Let (m_{1}, s(m_{1})), (m_{2}, s (m_{2})) , …, (m_{r}, s (m_{r})) be the observed accumulation curve and let x_{1}, x_{2}, …, x_{n} constitute “complete data”. We define the following notation:

For i = 1, 2, …, r and m_{i−1} < j ≤ m_{i}, let x_{ij} = x and

p_{ij} = \frac{b (α + 1, β + j − 1)}{b (α, β)}.

In addition, x_{i+} = ∑_{j} x_{ij} = s(m_{i}) − s (m_{i−1}) ; x_{++} = ∑_{i} x_{i+}; p_{i+} = ∑_{j} p_{ij}; and p_{++} =

∑_{i} p_{i+} .

With this notation, the likelihood for complete data is a product of multinomial densities:

\frac{k!}{(k − x_{++})! ∏_{i,j} x_{ij}!}p_{++}^{k−x_{++}} ∏_{i,j} p_{ij}^{x_{ij}}

= ∏_{i} [\frac{x_{ij}!}{∏_{j} x_{ij}!} ∏_{j} (\frac{p_{ij}}{p_{i+}} )^{x_{ij}} ] × \frac{k!}{(k − x_{++})! ∏_{i} x_{i+}!}p_{++}^{k−x_{++}} ∏_{i} p_{i+}^{x_{i+}}.

Therefore,

Q (θ′, θ) = E [log f (x | θ′) | y, θ]
= − ∑_{i} ∑_{j} E [ln (x_{ij}!) | y, θ] + ∑_{i} ∑_{j} ln (p_{ij}) E [x_{ij} | y, θ] + ln k!
−E [ln(k − x_{++})! | y, θ] + ln (p_{++}) E [(k − x_{++}) | y, θ] .

Because x_{i+} and x_{++} are actually observed, we have

Q (θ′, θ) = ln k! − ln(k − x_{++})! + ln (p_{++}) (k − x_{++})

+ ∑_{i} ∑_{j} ln (p_{ij}) E [x_{ij} | y, θ] − ∑_{i} ∑_{j} E [ln (x_{ij}!) | y, θ] .

On the other hand, x_{ij} has the Bin (x_{i+}, p_{ij}/p_{i+}) distribution. Therefore,

x_{ij}^{∗} := E [x_{ij} | y, θ] = \{ \begin{matrix}x_{ij}, if x_{ij} was observed \\ x_{i+} \frac{p_{ij}}{p_{i+}} , otherwise \end{matrix} .

Since the term ∑_{i} ∑_{j} E[ln (x_{ij}!) | y, θ] is only a function of observations and current parameter values, it can be ignored in the M step of the algorithm. In other words, to maximize Q (θ′, θ) , is equivalent to maximizing

Q^{∗} (θ′, θ) = ln k! − ln(k − x_{++})! + (k − x_{++}) ln (p_{++}) + ∑_{i} ∑_{j} ln (p_{ij}) x_{ij}^{∗}.

It is interesting to note that this Q^{∗} (θ′, θ) is very much related to equation (4), which is the log-likelihood of the complete model. Here, non-observed values were replaced by expected values.

References

[1] Barndorff-Nielsen, O. E. and Cox, D. R. (1994). Inference and Asymptotics. Monographs on Statistics and Applied Probability 52. Chapman & Hall.

[2] Bunge, J. and Fitzpatrick, M. (1993). Estimating the Number of Species: A Review. Journal of the American Statistical Association, 88, 364–373.

[3] Coldwell, R. K. and Coddington, J. A. (1994) . Estimating terrestrial biodiversity trough extrapolation. Phil. Trans. R. Soc. Lond. B. Special volume on “The quantification and estimation or organismal biodiversity”, (In press). D. L. Hawkswort, ed.

[4] Cox, D. R. and Reid, N. (1992). A note on then difference between profile and modified profile likelihood. Biometrika, 79, 2, 408–411.

[5] Dempster, A. P., Laird, N. M. and Rubin, D. B. (1977) . Maximum Likelihood from Incomplete Data via EM Algorithm. Journal of the Royal Statistical Society, B 39, 1–38.

[6] Efron, B. (1982). The Jackknife, the Bootstrap and Other Resampling Plans. Society for Industrial and Applied Mathematics. Philadelphia, Pennsylvania.

[7] Efron, B. and Tibshirani, R. (1986). Bootstrap Methods for Standard Errors, Confidence Intervals, and Other Measures of Statistical Accuracy. Statistical Science, 1, 1, 54–77.

[8] Efron, B. and Thisted, R. (1976). Estimating the number of unseen species: How many words did Shakespeare know? Biometrika, 63, 3, 435–447.

[9] Engen, S. (1978). Stochastic Abundance Models. Monographs on Applied Probability and Statistics. Chapman and Hall.

[10] Hinkley, D. V., Reid, N. and Snell, E. J. (1991) . Statistical Theory and Modelling. In Honour of Sir David Cox, FRS. Chapman and Hall.

[11] McCullagh, P. and Tibshirani, R. (1990). A Simple Method for the Adjustment of Profile Likelihoods. Journal of the Royal Statistical Society, B 52, 2, 325–344.

[12] Otis, D. L., Burnham, K. P., White, G. C. and Anderson, D.R. (1978) . Wildlife Monographs: Statistical Inference from Capture Data on Closed Animal Populations. The Wildlife Society, Inc.

[13] Press, W. H., Flannery, B. P., Teukolsky, S. A., and Vetterling, W. T. (1988). Numerical Recipes in C. Cambridge University Press.

[14] Seber, G. A. F. (1982). The Estimation of Animal Abundance and related parameters (Second Edition). Charles Griffin & Company Limited.

[15] Soberón, J. and Llorente, J. (1993). The Use of Species Accumulation Functions for the Prediction of Species Richness. Conservation Biology, 7, 3, 480–488.