ARTICLE IN PRESS Deep-Sea Research II 56 (2009) 1812–1815
Contents lists available at ScienceDirect
Deep-Sea Research II journal homepage: www.elsevier.com/locate/dsr2
Species number with confidence Woollcott K. Smith a,, Andrew R. Solow b, Nancy J. Maciolek c a b c
Statistics Department, Temple University, Philadelphia, PA 19122, United States Woods Hole Oceanographic Institution, Woods Hole, MA 02543, United States AECOM Environment, Marine & Coastal Center, Woods Hole, MA 02543, United States
a r t i c l e in f o
a b s t r a c t
Available online 31 May 2009
Benthic sampling by J. Frederick Grassle and Nancy J. Maciolek during the 1980s contributed to a recognition of the enormous diversity of biota in the deep sea. The problem of predicting on the basis of such sampling, the number of species in a community is a classical one in statistical ecology. Here, we show how to construct a lower prediction bound on species number using the sequential broken-stick model of relative abundances. We illustrate the method using some data from Grassle and Maciolek’s original sampling effort. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Benthic sampling Sequential broken-stick model Species richness
1. Introduction Over the past 40 years, recognition has grown of the enormous diversity of life in the deep sea. A signal contribution to this recognition was the benthic sampling conducted by Grassle and Maciolek (1992) during the 1980s at depths of 1500–2500 m off Delaware and New Jersey on the east coast of the US. In a total sample area of around 20 m2, Grassle and Maciolek collected approximately 90,000 individual organisms representing nearly 800 species, of which almost 60% were new to science. While this remains one of the most extensive quantitative sampling efforts in the deep sea, it represents the merest glimpse of deep-sea diversity. The problem of predicting the number of species in a community based on a small sample of individuals is a classical problem in statistical ecology. Excellent reviews are provided in Bunge and Fitzpatrick (1993) and Chao (2005). Briefly, under the standard formulation of this problem, a total of n individuals are sampled at random from a community of an unknown number S of species with unknown relative abundance vector p ¼ (p1,p2,y,pS). Importantly, the species represented in the sample are unordered, in the sense that the correspondence to the elements of p is unknown. Thus, the observed counts of different species represent an unordered multinomial sample of known size n with an unknown probability vector p. Our interest centers on using the results of this sampling to predict S. Both parametric and nonparametric approaches to prediction have been proposed in this setting. Here, an approach is parametric if it is based on a parametric model of p. One popular Corresponding author.
E-mail addresses:
[email protected] (W.K. Smith),
[email protected] (A.R. Solow),
[email protected] (N.J. Maciolek). 0967-0645/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.dsr2.2009.05.030
approach to predicting S is to fit and extrapolate a model to the sample species accumulation curve of Smith and Grassle (1977)—see, for example, Colwell and Coddington (1994). As appealing as this seems, the approach has no basis in statistical theory. Moreover, as a practical matter, its results are highly sensitive to the choice of extrapolating models, many of which fit the sample species accumulation curve equally well. Chao (1984) and Chao and Lee (1992) proposed nonparametric estimators of a lower bound on S. Not surprisingly, these tend to under-predict S, in some cases by a wide margin. Parametric methods based on the Dirichlet distribution for p have been proposed by Lewins and Joanes (1984) and Boender and Rinnooy Kan (1987). Recently, Solow and Smith (2009) proposed a method based on the sequential broken-stick model for relative abundances (Sugihara, 1980). We discuss this model in more detail below. In many practical situations, n is relatively small, S is large, and the elements of p are uneven. This is certainly the case for the deep-sea sampling of Grassle and Maciolek (1992). In such situations, no point predictor of S will be precise, and it may be more useful to have a prediction interval for S. The purpose of this paper is to describe and illustrate the construction of a parametric prediction interval under the sequential broken-stick model. The basic idea was described in Solow and Smith (2009), but the main focus in that paper was on point prediction.
2. Method In this section, we outline the construction of a one-sided lower 1a prediction interval for S of the form (SL, N) under the sequential broken-stick model. The lower bound SL represents the smallest value below which we can be confident at level 1a that
ARTICLE IN PRESS W.K. Smith et al. / Deep-Sea Research II 56 (2009) 1812–1815
1813
S does not fall. The same general approach can be used to construct a one-sided upper or two-sided prediction interval. The bound SL is given by the smallest value of So for which the null hypothesis Ho:S ¼ SL cannot be rejected at significance level a against the one-sided alternative hypothesis H1:S4So. Consider a test statistic T of which large values favor H1 and let Tobs be the observed value of T. Specific choices of T are discussed below. It follows that the bound SL is given by the smallest value of So that satisfies probðT T obs jS ¼ So Þ4a
(1)
Under the multinomial sampling model outlined in the previous section Z probðT T obs jS ¼ So Þ ¼ probðT T obs jpÞ pðpjS ¼ So Þ dp (2) where p(p|S ¼ So) is the probability density function of the relative abundance vector p when S ¼ So and the integral is over all possible such vectors. Along with n, the specification of p determines the distribution of T. Even in the absence of an analytical result, it is straightforward to approximate prob(TrTobs|p) by repeatedly simulating a multinomial sample, calculating the corresponding value of T, and determining whether or not it exceeds Tobs. Thus, the key element in (2) is the abundance distribution p(p|S ¼ So). Following Solow and Smith (2009), we adopt the sequential broken-stick model to which we now turn. Sugihara (1980) proposed the sequential broken-stick model for p(p|S ¼ So) based on a model of resource apportionment. This model has received some empirical support (e.g., Nee et al., 1991; Sugihara et al., 2003; Solow and Smith, 2009). The model, which is also called the Kolmogorov breakage model, is defined mechanistically in the following way. Break the unit interval into two pieces at random. Select one of the two pieces at random and break it into two pieces at random. There are now three pieces. Select one of them at random and break it into two pieces at random. Continue the process until there are S pieces. The relative abundances are given by the lengths of these S pieces. Siegel and Sugihara (1983) gave some theoretical results for this model, but the induced probability density p(p|S ¼ So) is not available. However, it is a simple matter to simulate relative abundance vectors from it. This suggests the following algorithm for approximating prob(TrTobs|S ¼ So): (i) simulate a relative abundance vector p from the sequential broken-stick model with S ¼ So; (ii) using the simulated value of p, evaluate prob(TrTobs|p), if necessary by simulation as outlined above; (iii) repeat the procedure a large number of times and approximate prob(TrTobs|S ¼ So) by the average of the values of prob(TrTobs|p) found in step (ii). To implement this approach, it is necessary to specify a test statistic T. The lack of theoretical results for the sequential broken-stick model makes this choice somewhat ad hoc. One natural choice is the observed number Sobs of species in the sample. The conditional distribution of Sobs given p is given in Emigh (1983), so that prob(TrTobs|S ¼ So) can in principle be found exactly. Unfortunately, its calculation is numerically P challenging. Let l ¼ Sj¼1 ð1 pj Þn Emigh (1983) noted that the conditional distribution of the number SSobs of unobserved species given p is approximately Poisson with mean l, thereby providing an approximation to the conditional distribution of Sobs. This approximation works well only if n is large and the largest element of p is relatively small. Here, we will approximate prob(TrTobs|p) by direct simulation.
Fig. 1. Lower bound for the 90% prediction interval for species number as a function of Sobs for selected values of n.
In Fig. 1, the lower 0.90 prediction bound SL is plotted against Sobs for selected values of n. So, for example, if, in a sample of 400 individuals, 100 species are observed, then with confidence level 0.90 the true number of species in the community is at least 620. On the other hand, if 100 species are observed in a sample of 800 individuals, this bound falls to around 300. Another natural choice of T is the bias-corrected version of the predictor of Chao (1984). This is given by f ðf 1Þ S^ Chao ¼ Sobs þ 1 1 2ðf 2 þ 1Þ
(3)
where f1 is the number of species in the sample represented by a single individual and f2 is the number represented by two individuals. As in the case of Sobs, it is straightforward to approximate the sampling distribution of S^ Chao by the simulation method outlined above. Fig. 2 repeats Fig. 1 using S^ Chao instead of Sobs. It is notable that, in all cases, the lower 0.90 confidence bound is actually larger than S^ Chao itself: that is, under the sequential broken-stick model, S exceeds S^ Chao with confidence greater than 0.90. This underscores the tendency of S^ Chao to underpredict S. It is not possible to determine from these figures which test statistic provides the larger lower bound. In ongoing work, we are studying the issue of optimizing the choice of test statistic. All species richness estimation methods depend on model structure. This is clearly the case for the sequential broken-stick model, where species with small or large relative abundances, pi, are equally likely to be split. Suppose we modify the sequential broken-stick model in the following way: at the n-th step in the sequential model the probability that species i is split is Prðspecies i is split at break nÞ ¼
pyi n P j¼1
.
(4)
pyj
The new parameter, 0ryr1, controls the ‘‘evenness’’ of the relative abundance distribution. For y ¼ 0 this is the sequential broken-stick model and for y ¼ 1 the probability that species i is split is just equal to the relative abundance of species i, pi. The original broken-stick model of MacArthur (1957) is exactly this model with y ¼ 1. For fixed y, the 90% confidence bound can be
ARTICLE IN PRESS 1814
W.K. Smith et al. / Deep-Sea Research II 56 (2009) 1812–1815
Fig. 2. As in Fig. 1 when the prediction bound is based on S^ Chao .
found using the method we have outlined for the sequential broken-stick model. Fig. 3 plots the lower limit for 90% confidence interval for a sample of size for n ¼ 500 and y ¼ 0.0, 0.1, 0.2, 0.5 and 1.0. This clearly indicates that the proposed confidence interval method is not robust to changes in y. In principle, one can use the evenness of the observed abundance data to estimate y before constructing a confidence interval for species richness. The details of this approach are complex and will be covered in a subsequent paper.
Fig. 3. Modified broken-stick model using species-break selection probabilities given in Eq. (4). Plotted is the lower bound for the 90% prediction interval for different values of y and a sample size of n ¼ 500. Note that y ¼ 1 corresponds to the original MacArthur (1957) broken-stick model and that y ¼ 0 corresponds to the sequential broken-stick model, Sugihara (1980).
Table 1 Summary statistics for a single box core and for the pooled box cores at station M9.
Single box core Pooled box cores
n
Sobs
f1
f2
S^ Chao
320 5956
93 278
44 74
13 33
161 357
Table 2 Lower 0.90 prediction bound for species number based on a single box core and on the pooled box cores at station M9 and on Sobs and S^ . Chao
3. Application Beyond the light it has shed on biological diversity in the deep sea, the sampling effort of Grassle and Maciolek (1992) is a rich source of data for the development of statistical methods. Grassle and Maciolek (1992) collected samples at a total of 10 stations. At each station, three replicate box cores were taken at each of six different times. Finally, nine sub-cores were taken within each box core and used for infaunal analysis. As an illustration of the method presented here, we analyzed two sets of data from Station M9. The first set of data was taken from a single box core from May 1984 and the second consisted of the data pooled over all times and replicates. Summary statistics for these data sets are presented in Table 1. Based on these summary statistics, we found lower 0.90 prediction bounds SL for species richness based on the observed number Sobs and on the bias-corrected predictor S^ Chao. In each case, the conditional distribution of the test statistic was approximated from 10,000 relative abundance vectors simulated from the sequential broken-stick model. The results are presented in Table 2. All of the bounds cover the overall total of 798 species collected by Grassle and Maciolek (1992). Although the bounds based on Sobs and S^ Chao are not terribly different, in both cases the former are larger than the latter. As before, it is notable that the lower bound based on S^ Chao is rather substantially larger than S^ Chao itself, underlining the conservative nature of this predictor.
Single box core Pooled box cores
Sobs
S^ Chao
690 700
630 610
4. Discussion The construction of a prediction interval—as opposed to a point prediction—for species number essentially requires a model of the distribution of relative abundances. In this paper, we have focused on the sequential broken-stick model for the purposes of illustrating our approach. In applied work, no model should be adopted uncritically. Solow and Smith (2009) described a method for assessing the goodness of fit of the sequential broken-stick model. However, the general approach outlined can be used in conjunction with any model of relative abundances. A minor technical issue concerns the choice of test statistic on which to base the prediction interval. This determination is difficult for the sequential broken-stick model because the model admits no analytical results. A more serious problem is the underlying multinomial model. This model arises if individual organisms are sampled at random. When area sampling is used, patchiness in the distributions of individual species can invalidate the multinomial model. This problem is shared with virtually all formal statistical approaches to inference about species number
ARTICLE IN PRESS W.K. Smith et al. / Deep-Sea Research II 56 (2009) 1812–1815
and the development of methods that accommodate patchiness, which is an area of ongoing research.
Acknowledgments We thank two anonymous referees for comments that led to the discussion of the sensitivity of our procedure to model assumptions and the results displayed in Fig. 3. References Boender, C.G.E., Rinnooy Kan, A.H.G., 1987. A multinomial Bayesian approach to the estimation of population and vocabulary size. Biometrika 74, 849–856. Bunge, J., Fitzpatrick, M., 1993. Estimating the number of species: a review. Journal of the American Statistical Association 88, 364–373. Chao, A., 1984. Nonparametric estimation of the number of classes in a population. Scandinavian Journal of Statistics 11, 265–270. Chao, A., Lee, S., 1992. Estimating the number of classes via sample coverage. Journal of the American Statistical Association 87, 210–217. Chao, A., 2005. Species estimation and applications. In: Balakrishnan, C.B., Read, C.B., Vidakovic, B. (Eds.), Encyclopedia of Statistical Sciences. Wiley, New York, pp. 7907–7916.
1815
Colwell, R.K., Coddington, J.A., 1994. Estimating terrestrial biodiversity through extrapolation. Philosophical Transactions of the Royal Society of London 345, 101–118. Emigh, T.H., 1983. On the number of observed classes from a multinomial distribution. Biometrics 39, 485–491. Grassle, J.F., Maciolek, N.J., 1992. Deep-sea species richness: regional and local diversity estimates from quantitative bottom samples. American Naturalist 139, 313–341. Lewins, W.A., Joanes, D.N., 1984. Bayesian estimation of the number of species. Biometrics 40, 323–328. MacArthur, R.H., 1957. On the relative abundance of bird species. Proceeding of the National Academy of Sciences 43, 293–295. Nee, S., Harvey, P.H., May, R.M., 1991. Lifting the veil on abundance patterns. Proceedings of the Royal Society of London B293, 161–163. Siegel, A.F., Sugihara, G., 1983. Moments of particle size distributions under sequential breakage with applications to species abundance. Journal of Applied Probability 20, 158–164. Smith, W.K., Grassle, J.F., 1977. Sampling properties of a family of diversity measures. Biometrics 33, 283–292. Solow, A.R., Smith, W.K., 2009. Estimating species number under an inconvenient abundance model. Journal of Agricultural, Biological, and Ecological Statistics 14, 242–252. Sugihara, G., 1980. Minimal community structure: an explanation of species abundance patterns. American Naturalist 116, 770–787. Sugihara, G., Bersier, L.-F., Southwood, T.R.E., Pimm, S.L., May, R.M., 2003. Predicted correspondences between species abundances and dendograms of niche similarities. Proceedings of the National Academy of Sciences 100, 5246–5251.