Fitting the truncated negative binomial distribution to count data

Abstract

Modeling empirical distributions of repeated counts with parametric probability distributions is a frequent problem when studying species abundance. One must choose a family of distributions which is flexible enough to take into account very diverse patterns and possess parameters with clear biological/ecological interpretations. The negative binomial distribution fulfills these criteria and was selected for modeling counts of marine fish and invertebrates. This distribution depends on a vector $\left( K,\mathfrak {P}\right)$ of parameters, and ranges from the Poisson distribution (when $K\rightarrow +\infty$) to Fisher’s log-series, when $K\rightarrow 0$. Moreover, these parameters have biological/ecological interpretations which are detailed in the literature and in this study. We compared three estimators of K, $\mathfrak {P}$ and the parameter $\alpha$ of Fisher’s log-series, following the work of Rao CR (Statistical ecology. Pennsylvania State University Press, University Park, 1971) on a three-parameter unstandardized variant of the negative binomial distribution. We further investigated the coherence underlying parameter values resulting from the different estimators, using both real count data collected in the Mauritanian Exclusive Economic Zone (MEEZ) during the period 1987–2010 and realistic simulations of these data. In the case of the MEEZ, we first built homogeneous lists of counts (replicates), by gathering observations of each species with respect to “typical environments” obtained by clustering the sampled stations. The best estimation of $\left( K,\mathfrak {P}\right)$ was generally obtained by penalized minimum Hellinger distance estimation. Interestingly, the parameters of most of the correctly sampled species seem compatible with the classical birth-and-dead model of population growth with immigration by Kendall (Biometrika 35:6–15, 1948).

Appendices

Appendix 1: Determination of a threshold for the truncated Hellinger distance

While results about the asymptotic distribution of our estimators abound, nothing is known about the distribution of the goodness-of-fit index [see Formula (9)]

$$\begin{aligned} d_{H}^{T}\left( TNBD\left( \hat{K},\mathfrak {\hat{P}}\right) ,\,TNBD\left( K,\mathfrak {P}\right) \right) \end{aligned}$$

(10)

where $\left( \hat{K},\mathfrak {\hat{P}}\right)$ is an estimate of $\left( K,\mathfrak {P}\right)$. In order to determine from $d_{H}^{T}$ the species which were correctly fitted, we performed a Monte Carlo study. It consisted in generating a sample of the statistics (10) for each one of the three estimators used, from a population of “negative binomial species” similar to the genuine population of the C4 habitat considered as a reference structure. This study is detailed hereunder.

1.1 The reference distribution of $\left( K,\mathfrak {P}\right)$

We plotted on Fig. 7 the minimum $PHD_{h}$ estimates of the vector $\left( K,\mathfrak {P}\right)$ associated with the species collected in the C4 habitat. About 35.6 % of the species were associated with very small values of the first parameter ($\hat{K}\le e^{-10}$); discarding these species, we could fit a bi-dimensional log-normal distribution of parameters $\left( \mu _{B},\Sigma _{B}\right)$ to the remaining vectors of estimates, whose confidence ellipsoids are also represented on Fig. 7. Neither $Log\left( \hat{K}\right)$ nor $Log\left( \mathfrak {\hat{P}}\right)$ strictly obeyed a normal distribution, but this model was retained for the sake of simplicity, since the corresponding 95 % confidence region widely covers the data (see Fig. 7). As for the discarded species, we postulated that $Log\left( \mathfrak {\hat{P}}\right)$ could be considered as obeying some Gaussian distribution, $\mathcal {N}\left( \mu _{S},\sigma _{S}\right)$.

Fig. 7

Fit of the estimated parameters for the C4 data. The ellipsoids correspond to 50 and 95 % confidence regions for the reference distribution, $\mathcal {N}\left( \mu _{B},\Sigma _{B}\right)$

1.2 Generating a “population” consistent with the reference distribution

To build a sample of $d_{H}^{T}\left( TNBD\left( \hat{K},\mathfrak {\hat{P}}\right) ,\,TNBD\left( K,\mathfrak {P}\right) \right)$ for counts having the same overall characteristics as the C4 data, we generated random counts of 300 “NB species”, whose random parameters obeyed the mixture distribution

$$\begin{aligned} \mathcal {M}:=0.356\,\mathcal {U}\left( \left[ e^{-12},e^{-8}\right] \right) \otimes \mathcal {LN}\left( \mu _{S},\sigma _{S}\right) +0.644\,\mathcal {LN}\left( \mu _{B},\Sigma _{B}\right) \end{aligned}$$

where $\left( \mu _{S},\sigma _{S}\right) =\left( -0.701062,2.68525\right)$ and $\left( \mu _{B},\Sigma _{B}\right)$ were estimated from the C4 data. In practice, the parameters $\left( k,p\right)$ of each species were first drawn according to $\mathcal {M}$; then a sample of $\beta =3000$ (or $\beta =6000$ when $k\le e^{-8}$) counts obeying $NBD\left( k,p\right)$ were drawn. The simulated data were then processed the same way as the MEEZ ones.

Fig. 8

In black: 50 and 95 % confidence ellipsoids  for the reference distribution $\mathcal {N}\left( \mu _{B},\Sigma _{B}\right)$; in gray: same confidence ellipsoids for the distribution $\mathcal {N}\left( \widehat{\mu _{B}},\widehat{\Sigma _{B}}\right)$ obtained from the mixture distribution $\mathcal {M}$. Dots correspond to the parameters of the “NB species”, estimated by PHD $_{h}$

On Fig. 8, we superimposed to the parameters of the species (estimated by $PHD_{h}$), confidence ellipsoids of the reference distribution $\mathcal {LN}\left( \mu _{B},\Sigma _{B}\right)$ and of the distribution $\mathcal {LN}\left( \widehat{\mu _{B}},\widehat{\Sigma _{B}}\right)$, whose parameters were estimated from the independent draws of $\mathcal {M}$. It is worth noting that in this case, there was no significant difference between the empirical distribution of $\left( \hat{K},\mathfrak {\hat{P}}\right)$ and the reference distribution $\mathcal {LN}\left( \mu _{B},\Sigma _{B}\right)$ (P values: Cramer–Von Mises =0.544087; Pearson $\chi ^{2}$ = 0.523489).

1.3 Results

About 30 % (93) of the “species” were observed ${<}6$ times, and discarded. The goodness-of-fit density estimates for the remaining ones are plotted on Fig. 9. The three estimators perform equally well in the case of genuine TNB distributions, while the goodness-of-fit by ULSD is very different. Among the 207 remaining species, only 3 % (6) were better fitted by ULSD, while the best estimator for TNBD was $PHD_{h}$ (185 species: 89 % of the total), followed by MLE (10 species) and pseudo-MLE (6 species). Thus almost all “aggregative species” were discarded due to their rarity.

Fig. 9

Simulations: density estimates of the four goodness-of-fit criteria

The quantile of order 0.95 of the goodness-of-fit associated with $PHD_{h}$ was 0.531096; consequently, we considered that $d_{H}^{T} = 0.53$ is an approriate threshold value that would not be exceeded by genuine negative binomial species. This threshold was used in Sect. 6.

Appendix 2: The Williams–Rao’s condition and the estimation of $\beta$

Notice that the equality $\alpha =K\,\beta$ is implicit in Eq. (4), but that it not a constraint in the root finding of the pseudo-MLE system (6). Consequently, this relationship must be considered as a sign of consistency of the estimation of $\left( K,\mathfrak {P,\alpha }\right)$. If in addition, the first condition of (5) is fulfilled, we should also have: $\dfrac{\alpha _{LS}}{K\,\beta }\approx \dfrac{\alpha _{Rao}}{K\,\beta }\approx 1$. We investigated the validity of this relationship by performing 50 Monte Carlo experiments with “negative binomial species” random draws, in each one of four typical cases:

• the “mean” one: $\left( K,\mathfrak {P}\right) =\left( 1.193,73.15\right)$ is the mean of the bivariate Log-normal distribution fitting the C4 habitat non-aggregative species (see Sect. 5.2 and Appendix 1)

• the “common” one: $\left( K,\mathfrak {P}\right) =\left( 0.7767,14.43\right)$ is the spatial median (Serfling 2004) of the parameters of the simulated “NB species”; in this case, we chose $\beta =10^{5}$ as the sample size of each one of the 50 simulations

• an “aggregative” case: $\left( K,\mathfrak {P}\right) =\left( 0.0001,14.43\right)$, with $\beta =10^{7}$

• a “bell-shaped” case: $\left( K,\mathfrak {P}\right) =\left( 10,14.43\right)$, with $\beta =10^{4}$.

In these four cases, the best fit was obtained with $PHD_{h}$, and we observed that $\alpha _{LS}$ and $\alpha _{Rao}$ could be considered as normally distributed (according to the Cramer-von Mises test), and that the mean of $\alpha _{Rao}$ was always close to $K\,\beta$ (T test), while the relationship $\dfrac{\alpha _{LS}}{K\,\beta }\approx 1$ was verified only in the aggregative case (see Fig. 12). The equality $\alpha _{LS}=K\,\beta$ was unacceptable in the “common” case (see Fig. 11), as well as in the bell-shaped one (Fig. 10) and in the “mean” case (not shown).

Fig. 10

The Williams–Rao’s condition in the bell-shaped case

Fig. 11

The Williams–Rao’s condition in the “common” case

Fig. 12

The Williams–Rao’s condition in the aggregative case—the blue density is $\dfrac{\alpha _{Rao}}{K\,\beta }$; the yellow one is $\dfrac{\alpha _{LS}}{K\,\beta }$ (Color figure online)

From another side, dealing with real data, we often have to face an additional estimation task: because of the dubious nature of the collected zeros (if there are), the true value of $\beta$ is unknown and we only know the number $\beta ^{+}$ of strictly positive counts. Then, it is classical to estimate $\beta$ by:

$$\begin{aligned} \beta ^{\left( \hat{K},\mathfrak {\hat{P}}\right) }:=\dfrac{\beta ^{+}}{\left( 1-1/\left( 1+\mathfrak {\hat{P}}\right) ^{\hat{K}}\right) } \end{aligned}$$

(11)

where $\left( \hat{K},\mathfrak {\hat{P}}\right)$ is an estimation of the TNBD parameters, obtained either from pseudo-MLE or $PHD_{h}$. The bell-shaped case is problem-free, because the probability of zero is extremely weak.

Let us now examine the most interesting case: the “common” one. On Fig. 13, we plotted on the left panel KDE estimates of the 50 values of the expression (11) obtained with $\left( \hat{K},\mathfrak {\hat{P},\hat{\alpha }}\right) =\left( K_{Rao},\mathfrak {P}_{Rao},\alpha _{Rao}\right)$ and divided by $\beta$: they were very close to 1. On the right panel, we plotted the values of (11) obtained with $\left( \hat{K},\mathfrak {\hat{P}}\right) =\left( K_{PHD_{h}},\mathfrak {P}_{PHD_{h}}\right)$. Thus, this figure shows that all estimates of $\beta$ are excellent in the “common” case; as a consequence, the results above concerning the Williams–Rao’s condition (see Fig. 11) stay valid in this case.

Fig. 13

Estimating $\beta$: the “common” case

Quite similar results were obtained in the “mean” case, as the reader can see on Fig. 14.

Fig. 14

Estimating $\beta$: the “mean” case

Things are very different in the aggregative situation, as the reader can see on Fig. 15: the estimator (11) based on $\left( K_{Rao},\mathfrak {P}_{Rao},\alpha _{Rao}\right)$ strongly underestimates or overestimates $\beta$. More precisely, in about 40 % of the samples, $\beta$ was highly overestimated, which is natural, since $\lim _{\hat{K}\rightarrow 0}\beta ^{\left( \hat{K},\mathfrak {\hat{P}}\right) }=\dfrac{\beta ^{+}}{\hat{K}\,\ln \left( 1+\mathfrak {\hat{P}}\right) }$. As for the estimator (11) associated with $\left( \hat{K},\mathfrak {\hat{P}}\right) =\left( K_{PHD_{h}},\mathfrak {P}_{PHD_{h}}\right)$, it always underestimated $\beta$. Consequently, when K is very small, the consistency of the estimations of the parameters of $UNBD\left( K,\mathfrak {P,\alpha }\right)$ is questionable and the log-series model is probably more sound—even if the fit is not quite as good as with $UNBD\left( K,\mathfrak {P,\alpha }\right)$. This is undoubtedly the meaning of Fig. 12.

Fig. 15

Estimating $\beta$ in the aggregative case; the 5 % upper values were dropped