A Penalized Empirical Likelihood Approach for Estimating Population Sizes under the Negative Binomial Regression Model

Abstract

In capture–recapture experiments, the presence of overdispersion and heterogeneity necessitates the use of the negative binomial regression model for inferring population sizes. However, within this model, existing methods based on likelihood and ratio regression for estimating the dispersion parameter often face boundary and nonidentifiability issues. These problems can result in nonsensically large point estimates and unbounded upper limits of confidence intervals for the population size. We present a penalized empirical likelihood technique for solving these two problems by imposing a half-normal prior on the population size. Based on the proposed approach, a maximum penalized empirical likelihood estimator with asymptotic normality and a penalized empirical likelihood ratio statistic with asymptotic chi-square distribution are derived. To improve numerical performance, we present an effective expectation-maximization (EM) algorithm. In the M-step, optimization for the model parameters could be achieved by fitting a standard negative binomial regression model via the R basic function glm.nb(). This approach ensures the convergence and reliability of the numerical algorithm. Using simulations, we analyze several synthetic datasets to illustrate three advantages of our methods in finite-sample cases: complete mitigation of the boundary problem, more efficient maximum penalized empirical likelihood estimates, and more precise penalized empirical likelihood ratio interval estimates compared to the estimates obtained without penalty. These advantages are further demonstrated in a case study estimating the abundance of black bears (Ursus americanus) at the U.S. Army’s Fort Drum Military Installation in northern New York.

1. Introduction

The use of capture–recapture data to infer population sizes is widely sought after across multiple fields, such as ecology and epidemiology. The capture–recapture sampling dates back to fisheries [1,2] and is used to estimate the abundance of species in the ecosystem [3]. Additionally, the number of populations of drug or alcohol addicts can be inferred by considering each visit to institutions as one capture [4,5]. Similarly, epidemiologists may leverage multiple incomplete lists of patients to estimate disease prevalence [6].

In traditional practice, the Poisson regression model is commonly used for modeling the capture–recapture data. However, the literature has demonstrated that the Poisson regression model may perform inadequately when the number of captures is subject to overdispersion, which means that the corresponding expectation is larger than variance. In capture–recapture experiments, the overdispersion often arises from individual heterogeneity. For example, animal species may exhibit flocking and spatiotemporal aggregation behaviors [7]. As argued in ref. [8], the population size may be underestimated if the model used ignores heterogeneity and overdispersion. Without any adjustments, using the Poisson regression model can yield unreliable statistical inferences, leading to incorrect interpretations. In this situation, ref. [9] argued that the negative binomial regression model might be more suitable.

Various estimation methods have been investigated for capture–recapture data within the framework of the negative binomial regression model. Based on conditional likelihood, ref. [8] introduced a maximum likelihood estimator for the regression parameter. Additionally, they derived a Horvitz–Thompson type estimator and a corresponding Wald-type interval estimator for estimating the population size. Building on the recapture probability ratios, refs. [10,11,12] investigated issues related to model diagnostics and population size estimation. Their simulations and empirical investigations revealed that the negative binomial regression model outperforms the Poisson regression model in the presence of overdispersion. However, these methods face challenges, and at least two specific issues arise. First, the likelihood function may suffer from identification and boundary problems for the dispersion parameter, leading to nonsensically large values of the Horvitz–Thompson type estimator associated with large standard errors; see, for example, refs. [10,11,13]. Second, the Newton–Raphson procedure may not be reliable and even fail to converge for maximizing the conditional likelihood function, as illustrated in refs. [14,15].

Empirical likelihood offers an alternative approach to mitigating the two issues mentioned earlier. Using empirical likelihood techniques, ref. [16] investigated the full semiparametric likelihood inference for population sizes in capture–recapture studies. Additionally, this method was extended by refs. [17,18,19,20] to include the continuous time case, one inflation and missing covariates. Experience shows that the semiparametric empirical-likelihood-based method usually outperforms the conditional-likelihood-based method in numerical experiments.

However, our simulation studies indicate that the point and interval estimators derived using maximum empirical likelihood estimation are prone to boundary issues, particularly when a small proportion of individuals are captured and severe overdispersion is present. For instance, as shown in our simulation studies with ${\mathsf{\beta }}_0={(-0.5,0.3)}^{\top }$, $k_0=0.5$, and $N_0=250$, where the capture probability is as low as 40%, nearly 67% of simulated cases yield unreasonable empirical likelihood ratio interval estimates whose upper limits exceed the number of captured individuals by over 100 times. This boundary problem is similarly observed in Section 4 for a case study. A potential cause of this issue may be that in such scenarios, the empirical log-likelihood ratio function flattens out as the population size increases.

The following sections are arranged as follows. Section 2 revisits the semiparametric empirical likelihood method within the negative binomial regression model and introduces a penalized empirical likelihood estimation approach to addressing the boundary problem. Furthermore, the EM algorithm is designed to enhance the reliability of the estimation process. Section 3 presents a number of simulations to illustrate how the proposed methods perform in the finite-sample settings. These methods are put into practice for analyzing the black bear data from New York in Section 4. Finally, a discussion is presented in Section 5.

2. Methodology

2.1. Model and Data

Suppose there are N individuals in a population. Each individual is characterized by the number of captures, denoted by Y, a nonnegative interger-valued variable. A naive approach to modeling Y is to assume a Poisson distribution, $Y\sim \mathrm{Poisson}\left(\theta \right)$, where $\theta >0$ represents the rate or expected number of captures. The Poisson model inherently assumes that all individuals in the population are homogeneous, meaning they share the same rate parameter $\theta$. In addition, as noted in the introduction, a key limitation of the Poisson model is its inability to handle overdispersed data, where the variance exceeds the mean.

To address the heterogeneity and overdispersion, one can assume that individuals have varying rates. In practice, a common approach is to model the distribution of these rates using a gamma distribution as a prior. Specifically, $\theta \sim \mathrm{Gamma}(r,p/(1-p\left)\right)$, where $r>0$ is the shape parameter and $p/(1-p)$ is the scale parameter with $p\in (0,1)$. With this prior, the marginal probability mass function of Y can be derived in a closed form:

$$\begin{array}{c}\hfill p_Y\left(y\right)={\int }_0^{\infty }p(y\mid \theta )p\left(\theta \right)d\theta ={\displaystyle \frac{\Gamma (y+r)}{\Gamma (y+1)\Gamma \left(r\right)}}{p}^y{(1-p)}^r,y=0,1,2,\cdots .\end{array}$$

which corresponds to the Poisson–Gamma or negative binomial distribution in the probability textbook. This distribution models the number of successes before the rth failure occurs in a sequence of independent Bernoulli trials, with p representing the probability of success in each trial.

As highlighted in ref. [9], a common reparameterization of $p_Y\left(y\right)$ is often used to interpret the counting processes in ecological and biodiversity studies. This reparameterization expresses $p_Y\left(y\right)$ in terms of its mean $\mu$ and a dispersion or aggregation parameter k, which controls the variation in counts. By setting $p=\mu /(k+\mu )$ and $k=r$, $p_Y\left(y\right)$ can be reformulated as:

$$\begin{array}{c}\hfill p_Y\left(y\right)={\displaystyle \frac{\Gamma (y+k)}{\Gamma (y+1)\Gamma \left(k\right)}}{\left({\displaystyle \frac{\mu }{k+\mu }}\right)}^y{\left({\displaystyle \frac{k}{k+\mu }}\right)}^k.\end{array}$$

When the individual covariates, denoted by $\mathit{X}$, are available, it becomes necessary to account for the heterogeneity induced by these covariates. To do so, a parametric model is used, specifically $\mu (\mathit{x};\beta )=exp\left({\mathsf{\beta }}^{\top }\mathit{x}\right)$, which relates the mean parameter to the individual covariates. Thus, given $\mathit{X}=\mathit{x}$, the conditional probability mass function of Y is expressed by:

$$\begin{array}{c}\hfill P(Y=y\mid \mathit{X}=\mathit{x})={\displaystyle \frac{\Gamma (y+k)}{\Gamma (y+1)\Gamma \left(k\right)}}{\left\{{\displaystyle \frac{\mu (\mathit{x};\mathsf{\beta })}{k+\mu (\mathit{x};\mathsf{\beta })}}\right\}}^y{\left\{{\displaystyle \frac{k}{k+\mu (\mathit{x};\mathsf{\beta })}}\right\}}^k=:f(y,\mathit{x};\mathsf{\beta },k),\end{array}$$

(1)

where $\mathsf{\beta }$ represents the unknown regression coefficients and $k>0$ represents the dispersion parameter. This formulation is referred to as the negative binomial regression model. As a special case, Equation (1) reduces to $p_Y\left(y\right)$ when all coefficients, except the intercept, are zero. The negative binomial regression model also includes the geometric regression model when $k=1$ and reduces to the Poisson regression model as $k\to \infty$.

Given $\mathit{X}=\mathit{x}$, the conditional expectation of Y is equal to $\mu (\mathit{x};\mathsf{\beta })$, while the conditional variance is expressed as $\mu (\mathit{x};\mathsf{\beta })+{\left\{\mu (\mathit{x};\mathsf{\beta })\right\}}^2/k$, indicating a quadratic relationship. The parameter k controls the degree of overdispersion: as k decreases, the variance increases, leading to greater overdispersion. Overdispersion is commonly observed in capture–recapture studies, where the variance significantly exceeds the mean. Consequently, the negative binomial regression model is often more appropriate for modeling capture–recapture data under conditions of severe overdispersion, as compared to the Poisson model.

Because the event $\{Y=0\}$ is unobservable in capture–recapture studies, the zero-truncated version of model (1) is considered:

$$\begin{array}{ccc}\hfill P(Y=y\mid Y>0,\mathit{X}=\mathit{x})& =& {\displaystyle \frac{P(Y=y\mid \mathit{X}=\mathit{x})}{P(Y>0\mid \mathit{X}=\mathit{x})}}\hfill \\ & =& {\displaystyle \frac{f(y,\mathit{x};\mathsf{\beta },k)}{1-\varphi (\mathit{x};\mathsf{\beta },k)}},y=1,2,\cdots ,\hfill \end{array}$$

(2)

where $\varphi (\mathit{x};\mathsf{\beta },k)={[k/\{k+\mu (\mathit{x};\mathsf{\beta })\}]}^k$ represents the conditional probability that an individual with a covariate $\mathit{x}$ is not captured at all.

Consider a study that captured n distinct individuals, with $({\mathit{x}}_1,\cdots ,{\mathit{x}}_n)$ and $(y_1,\cdots ,y_n)$ denoting their individual covariates and capture frequencies, respectively. Under the model (2), ref. [8] proposed a maximum conditional likelihood estimator $(\tilde{\mathsf{\beta }},\tilde{k})$ by maximizing:

$$\prod _{i=1}^n{\displaystyle \frac{f(y_i,{\mathit{x}}_i;\mathsf{\beta },k)}{1-\varphi ({\mathit{x}}_i;\mathsf{\beta },k)}}.$$

According to the principle of inverse probability weighting, the Horvitz–Thompson type estimator of N is defined as $\tilde{N}={\sum }_{i=1}^n{\{1-\varphi ({\mathit{x}}_i;\tilde{\mathsf{\beta }},\tilde{k})\}}^{-1}$. However, $\tilde{N}$ might be inflated due to small detection probabilities.

2.2. Semiparametric Empirical Likelihood

The semiparametric empirical likelihood, initially derived from ref. [16], is an appealing technique for implementing the full likelihood method when capture–recapture data contain individual covariates. Taking the negative binomial regression model as an example, we provide a brief introduction to this technique.

Considering that n distinct individuals out of a total of N individuals were captured, n follows a binomial distribution and the corresponding probability is as follows:

$$P\left(n\right)=\left({\displaystyle \genfrac{}{}{0pt}{}{N}{n}}\right){(1-\alpha )}^n{\alpha }^{N-n},$$

where $\alpha =P(Y=0)$ represents the probability that a generic individual was not captured at all. For the given n individuals, the conditional probability of their covariates and capture counts is as follows:

$$\prod _{i=1}^n{\displaystyle \frac{P(Y=y_i,\mathit{X}={\mathit{x}}_i)}{P(Y>0)}}=\prod _{i=1}^n{\displaystyle \frac{f(y_i,{\mathit{x}}_i;\mathsf{\beta },k)P(\mathit{X}={\mathit{x}}_i)}{1-\alpha }}.$$

Multiplying these two expressions yields the full likelihood function:

$$\begin{array}{c}\hfill \left({\displaystyle \genfrac{}{}{0pt}{}{N}{n}}\right){\alpha }^{N-n}\times \prod _{i=1}^n\left\{f(y_i,{\mathit{x}}_i;\mathsf{\beta },k)P(\mathit{X}={\mathit{x}}_i)\right\}.\end{array}$$

(3)

In Equation (3), the marginal probability $P(\mathit{X}={\mathit{x}}_i)$ is unknown and shall be addressed by the empirical likelihood method; see refs. [21,22] for more details. Technically, we assume that $P(\mathit{X}={\mathit{x}}_i)=p_i$ for $i=1,2,\cdots ,n$, where $p_i\in (0,1)$ is subject to the constraint ${\sum }_{i=1}^n{p}_i=1$. With this substitution, we call the full likelihood the semiparametric empirical likelihood and refer to its logarithm as the empirical log-likelihood function, namely:

$$\begin{array}{cc}\hfill \tilde{\mathsf{\ell }}(N,\mathsf{\beta },k,\alpha ,\left\{p_i\right\})& =log\left({\displaystyle \genfrac{}{}{0pt}{}{N}{n}}\right)+(N-n)log\left(\alpha \right)+\sum _{i=1}^{n}log\left\{f(y_i,{\mathit{x}}_i;\mathsf{\beta },k)\right\}\hfill \\ & +\sum _{i=1}^{n}log\left(p_i\right).\hfill \end{array}$$

By the definition of $\alpha$ and the iterated expectation theorem, it follows that $\alpha =E\left\{P\right(Y=0\mid \mathit{X}\left)\right\}$, or equivalently:

$$\sum _{i=1}^n\{\varphi ({\mathit{x}}_i;\mathsf{\beta },k)-\alpha \}{p}_i=0.$$

With the constraints for $p_i$’s, the profile empirical log-likelihood function can be derived using the Lagrange multiplier method on Equation (3):

$$\begin{array}{cc}\hfill \mathsf{\ell }(N,\mathsf{\beta },k,\alpha )=& log\left({\displaystyle \genfrac{}{}{0pt}{}{N}{n}}\right)+(N-n)log\left(\alpha \right)+\sum _{i=1}^{n}log\left({\displaystyle \genfrac{}{}{0pt}{}{y_i+k-1}{y_i}}\right)\hfill \\ & +\sum _{i=1}^n\left[klog\left\{{\displaystyle \frac{k}{k+\mu ({\mathit{x}}_i;\mathsf{\beta })}}\right\}+y_{i}log\left\{{\displaystyle \frac{\mu ({\mathit{x}}_i;\mathsf{\beta })}{k+\mu ({\mathit{x}}_i;\mathsf{\beta })}}\right\}\right]\hfill \\ & -\sum _{i=1}^{n}log[1+\xi \{\varphi ({\mathit{x}}_i;\mathsf{\beta },k)-\alpha \}]-nlog\left(n\right),\hfill \end{array}$$

where $\xi$ is the Lagrange multiplier, satisfying:

$$\sum _{i=1}^n{\displaystyle \frac{\varphi ({\mathit{x}}_i;\mathsf{\beta },k)-\alpha }{1+\xi \{\varphi ({\mathit{x}}_i;\mathsf{\beta },k)-\alpha \}}}=0.$$

Notice that there are a finite number of unknown parameters in the profile empirical log-likelihood function. By maximizing this function, we obtain the maximum empirical likelihood estimator, expressed as $(\widehat{N},\widehat{\mathsf{\beta }},\widehat{k},\widehat{\alpha })=argmax\left\{\mathsf{\ell }(N,\mathsf{\beta },k,\alpha )\right\}.$

2.3. Penalized Empirical Likelihood Inference

When the number of captures exhibits severe overdispersion, both the estimators $\tilde{N}$ and $\widehat{N}$, proposed in Section 2.1 and Section 2.2, respectively, may exhibit spuriously large values, potentially leading to misleading conclusions. This issue has been addressed in ref. [13] (p. 84) for the Horvitz–Thompson type estimator. Our simulation studies further confirm that the empirical likelihood estimators may also suffer from the boundary problem. This issue may arise due to the limited information available about the population size, causing the profile empirical log-likelihood to fail in distinguishing between different values of large N.

To mitigate this problem, we intuitively incorporate prior information on the population size to reduce the probability of large values. We achieve this by augmenting the likelihood functions with an appropriate penalty term. Correspondingly, the penalized empirical log-likelihood function and its profile version are defined as:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\tilde{\mathsf{\ell }}}_p(N,\mathsf{\beta },k,\alpha ,\left\{p_i\right\})& =\tilde{\mathsf{\ell }}(N,\mathsf{\beta },k,\alpha ,\left\{p_i\right\})+f_p\left(N\right),\hfill \\ \hfill {\mathsf{\ell }}_p(N,\mathsf{\beta },k,\alpha )& =\mathsf{\ell }(N,\mathsf{\beta },k,\alpha )+f_p\left(N\right).\hfill \end{array}\end{array}$$

(4)

where the penalty term $f_p\left(N\right)$ takes the form of $-C{(N-\nu )}^{2}I(N>\nu )$, where $\nu$ is a lower bound of N, $C\ge 0$ is a tuning parameter, and $I(·)$ is the indicator function. For specific values of $\nu$ and C, a maximum penalized empirical likelihood estimator is proposed, namely, $({\widehat{N}}_p,{\widehat{\mathsf{\beta }}}_p,{\widehat{k}}_p,{\widehat{\alpha }}_p)=argmax\left\{{\mathsf{\ell }}_p(N,\mathsf{\beta },k,\alpha )\right\}.$

From a Bayesian perspective, adding the penalty term $f_p\left(N\right)$ into the log-likelihood is equivalent to imposing a prior for N that has a mixture of the half-normal distribution $\mathrm{N}(\nu ,{\left(2C\right)}^{-1})$ for $N>\nu$ and a uniform distribution $\mathrm{U}(n,\nu )$ for $n\le N\le \nu$. In other words, this penalty has no effect on the likelihood when $n\le N\le \nu$ and gradually decreases the likelihood when $N>\nu$. The larger the population size, the more pronounced the decrease. Consequently, the penalized method encourages large values of $\widehat{N}$ to shrink towards $\nu$. In practice, we recommend using the Chao estimator as the lower bound $\nu$; see ref. [23] for details about this estimator. Alternatively, the generalized Chao estimator proposed in ref. [24] can also be considered.

To derive the large-sample properties of the estimator $({\widehat{N}}_p,{\widehat{\mathsf{\beta }}}_p,{\widehat{k}}_p,{\widehat{\alpha }}_p)$, we define some notation when the parameter vector $(N,\mathsf{\beta },k,\alpha )$ takes its true value, namely, $(N_0,{\mathsf{\beta }}_0,k_0,{\alpha }_0)$. Define $\psi (\mathit{x};{\mathsf{\beta }}_0,k_0)=-log\left\{1+\mu (\mathit{x};{\mathsf{\beta }}_0)/k_0\right\}+{\left\{1+k_0/\mu (\mathit{x};{\mathsf{\beta }}_0)\right\}}^{-1}$, $\phi =\mathbb{E}{\{1-\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)\}}^{-1}$, and:

$$\begin{array}{c}\hfill S_2(y+k_0-1,y)=\left\{\begin{array}{cc}-\sum _{k=k_0}^{y+k_0-1}{\displaystyle \frac{1}{k^2}},\hfill & y=1,2,\cdots \hfill \\ 0,\hfill & y=0\hfill \end{array}\right..\end{array}$$

Let:

$$\begin{array}{c}\hfill \mathbf{W}=\left[\begin{array}{cccc}-{\mathrm{V}}_{11}& {\mathbf{0}}_{s\times 1}^{\top }& 0& -{\mathrm{V}}_{14}\\ {\mathbf{0}}_{s\times 1}& -{\mathbf{V}}_{22}+{\mathbf{V}}_{25}{\mathrm{V}}_{55}^{-1}{\mathbf{V}}_{52}& -{\mathbf{V}}_{23}+{\mathbf{V}}_{25}{\mathrm{V}}_{55}^{-1}{\mathrm{V}}_{53}& -{\mathbf{V}}_{24}+{\mathbf{V}}_{25}{\mathrm{V}}_{55}^{-1}{\mathrm{V}}_{54}\\ 0& -{\mathbf{V}}_{32}+{\mathrm{V}}_{35}{\mathrm{V}}_{55}^{-1}{\mathbf{V}}_{52}& -{\mathrm{V}}_{33}+{\mathrm{V}}_{35}{\mathrm{V}}_{55}^{-1}{\mathrm{V}}_{53}& -{\mathrm{V}}_{34}+{\mathrm{V}}_{35}{\mathrm{V}}_{55}^{-1}{\mathrm{V}}_{54}\\ -{\mathrm{V}}_{41}& -{\mathbf{V}}_{42}+{\mathrm{V}}_{45}{\mathrm{V}}_{55}^{-1}{\mathbf{V}}_{52}& -{\mathrm{V}}_{43}+{\mathrm{V}}_{45}{\mathrm{V}}_{55}^{-1}{\mathrm{V}}_{53}& -{\mathrm{V}}_{44}+{\mathrm{V}}_{45}{\mathrm{V}}_{55}^{-1}{\mathrm{V}}_{54}\end{array}\right],\end{array}$$

(5)

where ${\mathrm{V}}_{11}=1-{\alpha }_0^{-1}$, ${\mathrm{V}}_{41}={\mathrm{V}}_{14}={\alpha }_0^{-1}$, ${\mathrm{V}}_{44}=\phi -{\alpha }_0^{-1}$, ${\mathrm{V}}_{45}={\mathrm{V}}_{54}={(1-{\alpha }_0)}^2\phi$, ${\mathrm{V}}_{55}={(1-{\alpha }_0)}^4\phi -{(1-{\alpha }_0)}^3$, and:

$$\begin{array}{ccc}\hfill {\mathbf{V}}_{22}& =& \mathbb{E}\left[\left\{{\displaystyle \frac{k_0\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)\mu (\mathit{X};{\mathsf{\beta }}_0)}{1-\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)}}-k_0-\mu (\mathit{X};{\mathsf{\beta }}_0,k_0)\right\}{\displaystyle \frac{k_0\mu (\mathit{X};{\mathsf{\beta }}_0){\mathit{X}}^{\otimes 2}}{{\{k_0+\mu (\mathit{X};{\mathsf{\beta }}_0)\}}^2}}\right],\hfill \\ \hfill {\mathbf{V}}_{23}& =& {\mathbf{V}}_{32}^{\top }=-\mathbb{E}\left[{\displaystyle \frac{k_0\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)\psi (\mathit{X};{\mathsf{\beta }}_0,k_0)\mu (\mathit{X};{\mathsf{\beta }}_0)\mathit{X}}{\{1-\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)\}\{k_0+\mu (\mathit{X};{\mathsf{\beta }}_0)\}}}\right],\hfill \\ \hfill {\mathrm{V}}_{33}& =& \mathbb{E}\left[{\displaystyle \frac{\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0){\psi }^2(\mathit{X};{\mathsf{\beta }}_0,k_0)}{1-\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)}}+{\displaystyle \frac{k_0\mu (\mathit{X};{\mathsf{\beta }}_0)+{\mu }^2(\mathit{X};{\mathsf{\beta }}_0)}{k_0{\{k_0+\mu (\mathit{X};{\mathsf{\beta }}_0)\}}^2}}\right]\hfill \\ & & +\mathbb{E}\left\{S_2(Y+k_0-1,Y)\right\},\hfill \\ \hfill {\mathbf{V}}_{24}& =& \mathbb{E}\left[{\displaystyle \frac{\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)}{1-\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)}}·{\displaystyle \frac{k_0\mu (\mathit{X};{\mathsf{\beta }}_0)\mathit{X}}{k_0+\mu (\mathit{X};{\mathsf{\beta }}_0)}}\right],{\mathbf{V}}_{25}={\mathbf{V}}_{52}^{\top }={(1-{\alpha }_0)}^2{\mathbf{V}}_{24},\hfill \\ \hfill {\mathrm{V}}_{34}& =& {\mathrm{V}}_{43}=-\mathbb{E}\left\{{\displaystyle \frac{\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)\psi (\mathit{X};{\mathsf{\beta }}_0,k_0)}{1-\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)}}\right\},{\mathrm{V}}_{35}={\mathrm{V}}_{53}={(1-{\alpha }_0)}^2{\mathrm{V}}_{34}.\hfill \end{array}$$

The following theorem presents the large-sample properties of the maximum penalized empirical likelihood estimator $({\widehat{N}}_p,{\widehat{\mathsf{\beta }}}_p,{\widehat{k}}_p,{\widehat{\alpha }}_p)$ associated with the penalized empirical likelihood ratio statistic of N.

Theorem 1.

Suppose that $k_0>0$, $0<{\alpha }_0<1$, and the tuning parameter satisfies $C=O\left(N^{-2}\right)$. If $\mathbf{W}$ is nonsingular and $\phi <\infty$, as $N_0\to \infty$:

(a) 

$\sqrt{N_0}{\{log\left({\widehat{N}}_p\right)-log\left(N_0\right),{({\widehat{\mathsf{\beta }}}_p-{\mathsf{\beta }}_0)}^{\top },{\widehat{k}}_p-k_0,{\widehat{\alpha }}_p-{\alpha }_0\}}^{\top }\stackrel{d}{⟶}\mathrm{N}(\mathbf{0},{\mathbf{W}}^{-1})$;

(b) 

$\sqrt{N_0}({\widehat{N}}_p/N_0-1)\stackrel{d}{⟶}N(0,{\sigma }^2)$, where

$$\begin{array}{c}\hfill {\sigma }^2=\phi -1-\left[\begin{array}{cc}{\mathbf{V}}_{42}& {\mathrm{V}}_{43}\end{array}\right]{\left[\begin{array}{cc}{\mathbf{V}}_{22}& {\mathbf{V}}_{23}\\ {\mathbf{V}}_{32}& {\mathrm{V}}_{33}\end{array}\right]}^{-1}\left[\begin{array}{c}{\mathbf{V}}_{24}\\ {\mathrm{V}}_{34}\end{array}\right];\end{array}$$

(c) 

$2\{{\mathsf{\ell }}_p({\widehat{N}}_p,{\widehat{\mathsf{\beta }}}_p,{\widehat{k}}_p,{\widehat{\alpha }}_p)-{max}_{(\mathsf{\beta },k,\alpha )}{\mathsf{\ell }}_p(N,\mathsf{\beta },k,\alpha )\}\stackrel{d}{⟶}{\chi }_1^2$, where ${\chi }_1^2$ denotes the standard chi-square distribution.

Proof. 

As the proposed semiparametric empirical likelihood approach can be seen as an extension of the EL method of ref. [16] to the negative binomial regression model, the proof of Theorem 1 is very similar to those of Theorem 1 and Corollary 1 in ref. [16]. Here, we only highlight the difference and the formulae of $\mathbf{V}$ and $\mathsf{\Sigma }$ in our framework.

We first argue that ${\xi }_0=-{(1-{\alpha }_0)}^{-1}$ is the limit of $\widehat{\xi }$, where $\widehat{\xi }$ is the solution to:

$$\sum _{i=1}^n{\displaystyle \frac{\varphi ({\mathit{x}}_i;{\widehat{\mathsf{\beta }}}_p,{\widehat{k}}_p)-{\widehat{\alpha }}_p}{1+\xi \{\varphi ({\mathit{x}}_i;{\widehat{\mathsf{\beta }}}_p,{\widehat{k}}_p)-{\widehat{\alpha }}_p\}}}=0.$$

For this purpose, we define a function:

$$\begin{array}{ccc}\hfill \hslash (N,\mathsf{\beta },k,\alpha ,\xi )& =& log\left({\displaystyle \genfrac{}{}{0pt}{}{N}{n}}\right)+(N-n)log\left(\alpha \right)-\sum _{i=1}^{n}log[1+\xi \{\varphi ({\mathit{x}}_i;\mathsf{\beta },k)-\alpha \}]\hfill \\ & & +\sum _{i=1}^n\left[log\left({\displaystyle \genfrac{}{}{0pt}{}{y_i+k-1}{y_i}}\right)+klog\left\{{\displaystyle \frac{k}{k+\mu ({\mathit{x}}_i;\mathsf{\beta })}}\right\}+y_{i}log\left\{{\displaystyle \frac{\mu ({\mathit{x}}_i;\mathsf{\beta })}{k+\mu ({\mathit{x}}_i;\mathsf{\beta })}}\right\}\right]\hfill \\ & & -nlog\left(n\right)+f_p\left(N\right).\hfill \end{array}$$

It can be seen that ${\mathsf{\ell }}_p(N,\mathsf{\beta },k,\alpha )=\hslash (N,\mathsf{\beta },k,\alpha ,{\xi }_{*}),$ where ${\xi }_{*}$ is the solution to $\mathsf{\partial }\hslash /\mathsf{\partial }\xi =0$. The first partial derivatives of $\hslash (N,\mathsf{\beta },k,\alpha ,\xi )$ with respect to $\alpha$ and $\xi$ are:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\mathsf{\partial }\hslash }{\mathsf{\partial }\alpha }}& =& {\displaystyle \frac{N-n}{\alpha }}+\xi \sum _{i=1}^n{\displaystyle \frac{1}{1+\xi \{\varphi ({\mathit{x}}_i;\mathsf{\beta },k)-\alpha \}}},\hfill \\ \hfill {\displaystyle \frac{\mathsf{\partial }\hslash }{\mathsf{\partial }\xi }}& =& -\sum _{i=1}^n{\displaystyle \frac{\varphi ({\mathit{x}}_i;\mathsf{\beta },k)-\alpha }{1+\xi \{\varphi ({\mathit{x}}_i;\mathsf{\beta },k)-\alpha \}}}.\hfill \end{array}$$

Setting the above equations to zero gives $\widehat{\xi }=-({\widehat{N}}_p-n)/\left(n{\widehat{\alpha }}_p\right).$ Since ${\alpha }_0$ is the probability of never being captured, n follows a Binomial distribution Bi($N_0,1-{\alpha }_0$). When $({\widehat{N}}_p,{\widehat{\alpha }}_p)$ is consistent, it follows from the strong law of large numbers that as $N_0\to \infty$:

$$\widehat{\xi }=-{\displaystyle \frac{\widehat{N}/N_0-n/N_0}{(n/N_0){\widehat{\alpha }}_p}}\stackrel{p}{⟶}-{\displaystyle \frac{1-(1-{\alpha }_0)}{(1-{\alpha }_0){\alpha }_0}}=-{\displaystyle \frac{1}{1-{\alpha }_0}}={\xi }_0,$$

where $\stackrel{p}{⟶}$ denotes convergence in probability.

Below, we derive the formulae of $\mathbf{V}$ and $\mathsf{\Sigma }$. Let $\mathsf{\gamma }={({\gamma }_1,{\mathsf{\gamma }}_2^{\top },{\gamma }_3,{\gamma }_4,{\gamma }_5)}^{\top }$, with ${\gamma }_1=\sqrt{N_0}(N/N_0-1),{\mathsf{\gamma }}_2=\sqrt{N_0}(\mathsf{\beta }-{\mathsf{\beta }}_0),{\gamma }_3=\sqrt{N_0}(k-k_0),{\gamma }_4=\sqrt{N_0}(\alpha -{\alpha }_0)$,

${\gamma }_5=\sqrt{N_0}(\xi -{\xi }_0)$. Define

$$\mathcal{H}\left(\mathsf{\gamma }\right)=\hslash (N_0+N_0^{1/2}{\gamma }_1,{\mathsf{\beta }}_0+N_0^{-1/2}{\mathsf{\gamma }}_2,k_0+N_0^{-1/2}{\gamma }_3,{\alpha }_0+N_0^{-1/2}{\gamma }_4,{\xi }_0+N_0^{-1/2}{\gamma }_5).$$

According to Lemma 2 in the Supplementary Material of [16], deriving the formula of $\mathbf{V}$ is equivalent to calculating the first two partial derivatives of $\mathcal{H}\left(\mathsf{\gamma }\right)$ with respect to $\mathsf{\gamma }$ at $\mathbf{0}$. It follows from the law of large numbers and the central limit theorem that:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\mathsf{\partial }\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_1}}& =& N_0^{1/2}\{S_1(N_0,n)+log\left({\alpha }_0\right)+f_p^{\prime }\left(N_0\right)\}=N_0^{1/2}\left({\displaystyle \frac{n/N_0-1}{{\alpha }_0}}+1\right)+O_p\left(N_0^{-1/2}\right),\hfill \\ \hfill {\displaystyle \frac{\mathsf{\partial }\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\mathsf{\gamma }}_2}}& =& N_0^{-1/2}\sum _{i=1}^n\left[\left\{-{\displaystyle \frac{\mu ({\mathit{x}}_i;{\mathsf{\beta }}_0)}{1-\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)}}+y_i\right\}{\displaystyle \frac{k_0{\mathit{x}}_i}{k_0+\mu ({\mathit{x}}_i;{\mathsf{\beta }}_0)}}\right],\hfill \\ \hfill {\displaystyle \frac{\mathsf{\partial }\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_3}}& =& N_0^{-1/2}\sum _{i=1}^n\left\{{\displaystyle \frac{\psi (\mathit{x};{\mathsf{\beta }}_0,k_0)}{1-\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)}}-{\displaystyle \frac{y_i}{k_0+\mu ({\mathit{x}}_i;{\mathsf{\beta }}_0)}}+S_1(y_i+k_0-1,y_i)\right\},\hfill \\ \hfill {\displaystyle \frac{\mathsf{\partial }\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_4}}& =& N_0^{-1/2}\left\{{\displaystyle \frac{N_0-n}{{\alpha }_0}}-\sum _{i=1}^n{\displaystyle \frac{1}{1-\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)}}\right\},\hfill \\ \hfill {\displaystyle \frac{\mathsf{\partial }\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_5}}& =& -N_0^{-1/2}(1-{\alpha }_0)\sum _{i=1}^n{\displaystyle \frac{\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)-{\alpha }_0}{1-\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)}},\hfill \end{array}$$

where the first equation uses the result (a) of Lemma A1 and Equation (A1) of Lemma A2 in the Appendix A. Using the result (b) of Lemma A1 and Equation (A2) of Lemma A2, we have:

$${\displaystyle \frac{{\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_1^2}}=N_0{S}_2(N_0,n)+N_0{f}_p^{″}\left(N_0\right)=1-{\alpha }_0^{-1}+O_p\left(N_0^{-1/2}\right).$$

Similarly, it can be verified that ${\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)/\left(\mathsf{\partial }{\gamma }_4\mathsf{\partial }{\gamma }_1\right)={\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)/\left(\mathsf{\partial }{\gamma }_1\mathsf{\partial }{\gamma }_4\right)={\alpha }_0^{-1}+O_p\left(N_0^{-1/2}\right)$ and:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{{\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\mathsf{\gamma }}_2\mathsf{\partial }{\mathsf{\gamma }}_2^{\top }}}& =& N_0^{-1}\sum _{i=1}^n\left(\left[{\displaystyle \frac{k_0\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)\mu ({\mathit{x}}_i;{\mathsf{\beta }}_0)}{{\{1-\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)\}}^2}}-{\displaystyle \frac{k_0}{1-\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)}}-y_i\right]{\displaystyle \frac{k_0\mu ({\mathit{x}}_i;{\mathsf{\beta }}_0){\mathit{x}}_i^{\otimes 2}}{{\{k_0+\mu ({\mathit{x}}_i;{\mathsf{\beta }}_0)\}}^2}}\right),\hfill \\ \hfill {\displaystyle \frac{{\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\mathsf{\gamma }}_2\mathsf{\partial }{\gamma }_3}}& =& {\left\{{\displaystyle \frac{{\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_3\mathsf{\partial }{\mathsf{\gamma }}_2^{\top }}}\right\}}^{\top }=N_0^{-1}\sum _{i=1}^n[-{\displaystyle \frac{k_0\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)\psi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)}{{\left\{1-\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)\right\}}^2}}+{\displaystyle \frac{y_i}{k_0+\mu ({\mathit{x}}_i;{\mathsf{\beta }}_0)}}\hfill \\ & & -{\displaystyle \frac{1}{1-\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)}}·{\displaystyle \frac{\mu ({\mathit{x}}_i;{\mathsf{\beta }}_0)}{k_0+\mu ({\mathit{x}}_i;{\mathsf{\beta }}_0)}}]{\displaystyle \frac{\mu ({\mathit{x}}_i;{\mathsf{\beta }}_0){\mathit{x}}_i}{k_0+\mu ({\mathit{x}}_i;{\mathsf{\beta }}_0)}},\hfill \\ \hfill {\displaystyle \frac{{\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\mathsf{\gamma }}_2\mathsf{\partial }{\gamma }_4}}& =& {\left\{{\displaystyle \frac{{\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_4\mathsf{\partial }{\mathsf{\gamma }}_2^{\top }}}\right\}}^{\top }=N_0^{-1}\sum _{i=1}^n{\displaystyle \frac{\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)}{{\{1-\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)\}}^2}}·{\displaystyle \frac{k_0\mu ({\mathit{x}}_i;{\mathsf{\beta }}_0){\mathit{x}}_i}{k_0+\mu ({\mathit{x}}_i;{\mathsf{\beta }}_0)}},\hfill \\ \hfill {\displaystyle \frac{{\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\mathsf{\gamma }}_2\mathsf{\partial }{\gamma }_5}}& =& {\left\{{\displaystyle \frac{{\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_5\mathsf{\partial }{\mathsf{\gamma }}_2^{\top }}}\right\}}^{\top }={(1-{\alpha }_0)}^2{N}_0^{-1}\sum _{i=1}^n{\displaystyle \frac{\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)}{{\{1-\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)\}}^2}}·{\displaystyle \frac{k_0\mu ({\mathit{x}}_i;{\mathsf{\beta }}_0){\mathit{x}}_i}{k_0+\mu ({\mathit{x}}_i;{\mathsf{\beta }}_0)}},\hfill \\ \hfill {\displaystyle \frac{{\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_3^2}}& =& N_0^{-1}\sum _{i=1}^n[{\displaystyle \frac{\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0){\psi }^2({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)}{{\{1-\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)\}}^2}}+{\displaystyle \frac{1}{1-\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)}}·{\displaystyle \frac{{\mu }^2({\mathit{x}}_i;{\mathsf{\beta }}_0)}{k_0{\{k_0+\mu ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)\}}^2}}\hfill \\ & & +{\displaystyle \frac{y_i}{{\{k_0+\mu ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)\}}^2}}+S_2(y_i+k_0-1,y_i)],\hfill \\ \hfill {\displaystyle \frac{{\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_3\mathsf{\partial }{\gamma }_4}}& =& {\displaystyle \frac{{\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_4\mathsf{\partial }{\gamma }_3}}=-N_0^{-1}\sum _{i=1}^n{\displaystyle \frac{\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)\psi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)}{{\{1-\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)\}}^2}},\hfill \\ \hfill {\displaystyle \frac{{\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_3\mathsf{\partial }{\gamma }_5}}& =& {\displaystyle \frac{{\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_5\mathsf{\partial }{\gamma }_3}}=-{(1-{\alpha }_0)}^2{N}_0^{-1}\sum _{i=1}^n{\displaystyle \frac{\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)\psi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)}{{\{1-\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)\}}^2}},\hfill \\ \hfill {\displaystyle \frac{{\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_4^2}}& =& N_0^{-1}\sum _{i=1}^n{\displaystyle \frac{1}{{\{1-\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)\}}^2}}-{\displaystyle \frac{1-n/N_0}{{\alpha }_0^2}},\hfill \\ \hfill {\displaystyle \frac{{\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_4\mathsf{\partial }{\gamma }_5}}& =& {\displaystyle \frac{{\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_5\mathsf{\partial }{\gamma }_4}}={(1-{\alpha }_0)}^2{N}_0^{-1}\sum _{i=1}^n{\displaystyle \frac{1}{{\{1-\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)\}}^2}},\hfill \\ \hfill {\displaystyle \frac{{\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_5^2}}& =& {(1-{\alpha }_0)}^4{N}_0^{-1}\sum _{i=1}^n{\displaystyle \frac{1}{{\{1-\varphi ({\mathit{x}}_i;{\mathsf{\beta }}_0,k_0)\}}^2}}-{(1-{\alpha }_0)}^3.\hfill \end{array}$$

With the arguments of Lemma A3 in the Appendix A, we have that the leading term of ${\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)/\left(\mathsf{\partial }\mathsf{\gamma }\mathsf{\partial }{\mathsf{\gamma }}^{\top }\right)$ is as follows:

$$\begin{array}{c}\hfill {\displaystyle \frac{{\mathsf{\partial }}^2\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }\mathsf{\gamma }\mathsf{\partial }{\mathsf{\gamma }}^{\top }}}=\mathbf{V}+O_p\left(N_0^{-1/2}\right),\mathbf{V}=\left[\begin{array}{ccccc}{\mathrm{V}}_{11}& {\mathbf{0}}^{\top }& 0& {\mathrm{V}}_{14}& 0\\ \mathbf{0}& {\mathbf{V}}_{22}& {\mathbf{V}}_{23}& {\mathbf{V}}_{24}& {\mathbf{V}}_{25}\\ 0& {\mathbf{V}}_{32}& {\mathrm{V}}_{33}& {\mathrm{V}}_{34}& {\mathrm{V}}_{35}\\ {\mathrm{V}}_{41}& {\mathbf{V}}_{42}& {\mathrm{V}}_{43}& {\mathrm{V}}_{44}& {\mathrm{V}}_{45}\\ 0& {\mathbf{V}}_{52}& {\mathrm{V}}_{53}& {\mathrm{V}}_{54}& {\mathrm{V}}_{55}\end{array}\right],\end{array}$$

where ${\mathrm{V}}_{ij}$s are the same as those defined in Equation (5).

Next, we derive the formula of $\mathsf{\Sigma }$. Define ${\mathbf{u}}_n={(u_{n1},{\mathbf{u}}_{n2}^{\top },u_{n3},u_{n4},u_{n5})}^{\top }$, where:

$$u_{n1}=N_0^{1/2}\left({\displaystyle \frac{n/N_0-1}{{\alpha }_0}}+1\right),{\mathbf{u}}_{n2}={\displaystyle \frac{\mathsf{\partial }\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\mathsf{\gamma }}_2}},u_{n3}={\displaystyle \frac{\mathsf{\partial }\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_3}},u_{n4}={\displaystyle \frac{\mathsf{\partial }\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_4}},u_{n5}={\displaystyle \frac{\mathsf{\partial }\mathcal{H}\left(\mathbf{0}\right)}{\mathsf{\partial }{\gamma }_5}}.$$

It can be verified that $\mathsf{\partial }\mathcal{H}\left(\mathbf{0}\right)/\mathsf{\partial }\mathsf{\gamma }={\mathbf{u}}_n+O_p\left(N_0^{-1/2}\right)$ and $\mathbb{E}\left({\mathbf{u}}_n\right)=\mathbf{0}$. Here, we only verify that $\mathbb{E}\left({\mathbf{u}}_{n2}\right)=\mathbf{0}$ and $\mathbb{E}\left(u_{n3}\right)=0$. In fact, it follows Lemma A3 that:

$$\begin{array}{ccc}\hfill \mathbb{E}\left({\mathbf{u}}_{n2}\right)& =& N_0^{1/2}\mathbb{E}\left\{-{\displaystyle \frac{k_0\mathit{X}\mu (\mathit{X};{\mathsf{\beta }}_0)}{k_0+\mu (\mathit{X};{\mathsf{\beta }}_0)}}+{\displaystyle \frac{k_0\mathit{X}\mathbb{E}(Y\mid \mathit{X})}{k_0+\mu (\mathit{X};{\mathsf{\beta }}_0)}}\right\}=\mathbf{0},\hfill \\ \hfill \mathbb{E}\left(u_{n3}\right)& =& N_0^{1/2}\mathbb{E}\left[\psi (\mathit{X};{\mathsf{\beta }}_0,k_0)-{\displaystyle \frac{\mathbb{E}(Y\mid \mathit{X})}{k_0+\mu (\mathit{X};{\mathsf{\beta }}_0)}}+\mathbb{E}\{S_1(Y+k_0-1,Y)\mid \mathit{X}\}\right]\hfill \\ & =& N_0^{1/2}\mathbb{E}\left[\psi (\mathit{X};{\mathsf{\beta }}_0,k_0)-{\displaystyle \frac{\mu (\mathit{X};{\mathsf{\beta }}_0)}{k_0+\mu (\mathit{X};{\mathsf{\beta }}_0)}}-log\left\{{\displaystyle \frac{k_0}{k_0+\mu (\mathit{X};{\mathsf{\beta }}_0)}}\right\}\right]=0,\hfill \end{array}$$

where the last second equation uses Equation (A3) of Lemma A4 in the Appendix A.

For the covariance matrix of ${\mathbf{u}}_n$, it follows Lemma A4 that:

$$\begin{array}{ccc}\hfill \mathbb{V}\mathrm{ar}\left(u_{n1}\right)& =& {\alpha }_0^{-1}-1,\mathbb{C}\mathrm{ov}(u_{n4},u_{n1})=-{\alpha }_0^{-1},\hfill \\ \hfill \mathbb{C}\mathrm{ov}({\mathbf{u}}_{n2},u_{n1})& =& \mathbb{C}\mathrm{ov}({\mathbf{u}}_{n2},u_{n4})=\mathbb{C}\mathrm{ov}({\mathbf{u}}_{n2},u_{n5})=\mathbf{0},\hfill \\ \hfill \mathbb{C}\mathrm{ov}(u_{n3},u_{n1})& =& \mathbb{C}\mathrm{ov}(u_{n3},u_{n4})=\mathbb{C}\mathrm{ov}(u_{n3},u_{n5})=\mathbb{C}\mathrm{ov}(u_{n5},u_{n1})=0,\hfill \\ \hfill \mathbb{V}\mathrm{ar}\left({\mathbf{u}}_{n2}\right)& =& \mathbb{E}{\left[{\displaystyle \frac{1}{N_0}}\sum _{i=1}^{N_0}\left\{-{\displaystyle \frac{\mu ({\mathit{X}}_i;{\mathsf{\beta }}_0)}{1-\varphi ({\mathit{X}}_i;{\mathsf{\beta }}_0,k_0)}}+Y_i\right\}{\displaystyle \frac{k_0{\mathit{X}}_{i}I(Y_i>0)}{k_0+\mu ({\mathit{X}}_i;{\mathsf{\beta }}_0)}}\right]}^{\otimes 2}\hfill \\ & =& \mathbb{E}\left[{\left\{-{\displaystyle \frac{\mu (\mathit{X};{\mathsf{\beta }}_0)}{1-\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)}}+Y\right\}}^2{\displaystyle \frac{k_0^{2}I(Y>0){\mathit{X}}^{\otimes 2}}{{\{k_0+\mu (\mathit{X};{\mathsf{\beta }}_0)\}}^2}}\right]\hfill \\ & =& \mathbb{E}\left[\left\{{\displaystyle \frac{{\mu }^2(\mathit{X};{\mathsf{\beta }}_0)}{{\{1-\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)\}}^2}}+Y^2-{\displaystyle \frac{2Y\mu (\mathit{X};{\mathsf{\beta }}_0)}{1-\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)}}\right\}{\displaystyle \frac{k_0^{2}I(Y>0){\mathit{X}}^{\otimes 2}}{{\{k_0+\mu (\mathit{X};{\mathsf{\beta }}_0)\}}^2}}\right]\hfill \\ & =& \mathbb{E}\left[\left\{-{\displaystyle \frac{k_0\mu (\mathit{X};{\mathsf{\beta }}_0)\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)}{1-\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)}}+k_0+\mu (\mathit{X};{\mathsf{\beta }}_0)\right\}{\displaystyle \frac{k_0\mu (\mathit{X};{\mathsf{\beta }}_0){\mathit{X}}^{\otimes 2}}{{\{k_0+\mu (\mathit{X};{\mathsf{\beta }}_0)\}}^2}}\right]\hfill \\ & =& -{\mathbf{V}}_{22},\hfill \\ \hfill \mathbb{C}\mathrm{ov}({\mathbf{u}}_{n2},u_{n3})& =& \mathbb{E}\left[\left\{-{\displaystyle \frac{\mu (\mathit{X};{\mathsf{\beta }}_0)}{1-\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)}}+Y\right\}\right\{{\displaystyle \frac{\psi (\mathit{x};{\mathsf{\beta }}_0,k_0)}{1-\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)}}-{\displaystyle \frac{Y}{k_0+\mu (\mathit{X};{\mathsf{\beta }}_0)}}\hfill \\ & & +S_1(Y+k_0-1,Y)\}·{\displaystyle \frac{k_0\mathit{X}I(Y>0)}{k_0+\mu (\mathit{X};{\mathsf{\beta }}_0)}}]\hfill \\ & =& \mathbb{E}\left[\right\{{\displaystyle \frac{\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)\psi (\mathit{x};{\mathsf{\beta }}_0,k_0)}{1-\varphi (\mathit{x};{\mathsf{\beta }}_0,k_0)}}-{\displaystyle \frac{1}{k_0}}+{\displaystyle \frac{\mathbb{E}[Y{S}_1(Y+k_0-1,Y)\mid \mathit{X}]}{\mu (\mathit{X};{\mathsf{\beta }}_0)}}\hfill \\ & & +log{\displaystyle \frac{k_0}{k_0+\mu (\mathit{X};{\mathsf{\beta }}_0)}}\left\}{\displaystyle \frac{k_0\mu (\mathit{X};{\mathsf{\beta }}_0)\mathit{X}}{k_0+\mu (\mathit{X};{\mathsf{\beta }}_0)}}\right]=-{\mathbf{V}}_{23},\hfill \\ \hfill \mathbb{V}\mathrm{ar}\left(u_{n3}\right)& =& \mathbb{E}\left[{\left\{{\displaystyle \frac{\psi (\mathit{X};{\mathsf{\beta }}_0,k_0)}{1-\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)}}-{\displaystyle \frac{Y}{k_0+\mu (\mathit{X};{\mathsf{\beta }}_0)}}+S_1(Y+k_0-1,Y)\right\}}^{2}I(Y>0)\right]\hfill \\ & =& \mathbb{E}[-{\displaystyle \frac{{\left\{\psi (\mathit{X};{\mathsf{\beta }}_0,k_0)\right\}}^2}{1-\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)}}+{\displaystyle \frac{(k_0+1){\left\{\mu (\mathit{X};{\mathsf{\beta }}_0,k_0)\right\}}^2+k_0\mu (\mathit{X};{\mathsf{\beta }}_0,k_0)}{k_0{\{k_0+\mu (\mathit{X};{\mathsf{\beta }}_0,k_0)\}}^2}}\hfill \\ & & +\mathbb{E}\left[{\left\{S_1(Y+k-1,Y)\right\}}^2\mid \mathit{X}\right]-{\displaystyle \frac{2\mathbb{E}\{Y{S}_1(Y+k-1,Y)\mid \mathit{X}\}}{k+\mu (\mathit{X};{\mathsf{\beta }}_0,k_0)}}]\hfill \\ & =& -{\mathbf{V}}_{33},\hfill \\ \hfill \mathbb{V}\mathrm{ar}\left(u_{n4}\right)& =& N_0^{-1}\mathbb{E}{\left[\sum _{i=1}^{N_0}\left\{{\displaystyle \frac{I(Y_i=0)}{{\alpha }_0}}-{\displaystyle \frac{I(Y_i>0)}{1-\varphi ({\mathit{X}}_i;{\mathsf{\beta }}_0,k_0)}}\right\}\right]}^2\hfill \\ & =& \mathbb{E}{\left\{{\displaystyle \frac{I(Y=0)}{{\alpha }_0}}-{\displaystyle \frac{I(Y>0)}{1-\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)}}\right\}}^2={\alpha }_0^{-1}+\mathbb{E}\left\{{\displaystyle \frac{1}{1-\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)}}\right\},\hfill \end{array}$$

$$\begin{array}{ccc}\hfill \mathbb{C}\mathrm{ov}(u_{n4},u_{n5})& =& \mathbb{E}\left[-{\displaystyle \frac{1-{\alpha }_0}{N_0}}\sum _{i=1}^{N_0}\left\{{\displaystyle \frac{I(Y_i=0)}{{\alpha }_0}}-{\displaystyle \frac{I(Y_i>0)}{1-\varphi ({\mathit{X}}_i;{\mathsf{\beta }}_0,k_0)}}\right\}\right.\hfill \\ & & \left.·{\displaystyle \frac{\{\varphi ({\mathit{X}}_i;{\mathsf{\beta }}_0,k_0)-{\alpha }_0\}I(Y_i>0)}{1-\varphi ({\mathit{X}}_i;{\mathsf{\beta }}_0,k_0)}}\right]\hfill \\ & =& -(1-{\alpha }_0)+{(1-{\alpha }_0)}^2\mathbb{E}\left\{{\displaystyle \frac{1}{1-\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)}}\right\},\hfill \\ \hfill \mathbb{V}\mathrm{ar}\left(u_{n5}\right)& =& N_0^{-1}{(1-{\alpha }_0)}^2\mathbb{E}\left[\sum _{i=1}^{N_0}{\displaystyle \frac{{\{\varphi ({\mathit{X}}_i;{\mathsf{\beta }}_0,k_0)-{\alpha }_0\}}^{2}I(Y_i>0)}{{\{1-\varphi ({\mathit{X}}_i;{\mathsf{\beta }}_0,k_0)\}}^2}}\right]\hfill \\ & =& {(1-{\alpha }_0)}^2\mathbb{E}\left[{\displaystyle \frac{{\{\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)-{\alpha }_0\}}^2}{1-\varphi (\mathit{X};{\mathsf{\beta }}_0,k_0)}}\right]={\mathrm{V}}_{55}.\hfill \end{array}$$

By the central limit theorem, as $N_0\to \infty$ we have ${\mathbf{u}}_n\stackrel{d}{⟶}\mathrm{N}(\mathbf{0},\mathsf{\Sigma })$, where:

$$\begin{array}{ccc}\hfill \mathsf{\Sigma }& =& \left(\begin{array}{ccccc}-{\mathrm{V}}_{11}& {\mathbf{0}}^{\top }& 0& -{\mathrm{V}}_{14}& 0\\ \mathbf{0}& -{\mathbf{V}}_{22}& -{\mathbf{V}}_{23}& \mathbf{0}& \mathbf{0}\\ 0& -{\mathbf{V}}_{32}& -{\mathrm{V}}_{33}& 0& 0\\ -{\mathrm{V}}_{41}& {\mathbf{0}}^{\top }& 0& 2{\mathrm{V}}_{45}{(1-{\alpha }_0)}^{-2}-{\mathrm{V}}_{44}& {\mathrm{V}}_{55}{(1-{\alpha }_0)}^{-2}\\ 0& {\mathbf{0}}^{\top }& 0& {\mathrm{V}}_{55}{(1-{\alpha }_0)}^{-2}& {\mathrm{V}}_{55}\end{array}\right).\hfill \end{array}$$

Since the covariance matrix $\mathsf{\Sigma }$ has the same form as the $\mathsf{\Sigma }$ in Lemma 3 of the Supplementary Material in ref. [16], so does the matrix $\mathbf{W}$. The rest of the proof is similar and omitted. This completes the proof of Theorem 1. □

When $C=0$, there is no penalty term and the likelihood functions with and without penalty coincide. This implies that the asymptotic results in Theorem 1 hold for the empirical likelihood estimators without penalty. Utilizing the result (c), a penalized empirical likelihood ratio interval estimator can be constructed, namely:

$$\begin{array}{c}\hfill {\mathcal{I}}_p=\left\{N:2\left\{{\mathsf{\ell }}_p({\widehat{N}}_p,{\widehat{\mathsf{\beta }}}_p,{\widehat{k}}_p,{\widehat{\alpha }}_p)-\underset{(\mathsf{\beta },k,\alpha )}{max}{\mathsf{\ell }}_p(N,\mathsf{\beta },k,\alpha )\right\}\le {\chi }_1^2(1-a)\right\},\end{array}$$

where ${\chi }_1^2(1-a)$ stands for the $(1-a)$th quantile of ${\chi }_1^2$. Correspondingly, the empirical likelihood ratio interval estimator derived without penalty is as follows:

$$\begin{array}{c}\hfill \mathcal{I}=\left\{N:2\left\{\mathsf{\ell }(\widehat{N},\widehat{\mathsf{\beta }},\widehat{k},\widehat{\alpha })-\underset{(\mathsf{\beta },k,\alpha )}{max}\mathsf{\ell }(N,\mathsf{\beta },k,\alpha )\right\}\le {\chi }_1^2(1-a)\right\}.\end{array}$$

Despite both interval estimators asymptotically yielding correct coverage probability $(1-a)$, our simulation studies indicate that ${\mathcal{I}}_p$ generally outperforms $\mathcal{I}$ in terms of interval width.

Remark 1.

One might question whether overdispersion exists, or equivalently, whether the zero-truncated Poisson regression model adequately fit the data. Various methods have been proposed to address this question. See, for instances, refs. [8,11,25,26].

2.4. Numerical Implementation

In this section, we aim to develop an EM algorithm to facilitate the proposed estimation method described in Section 2.3. For better presentation, we begin by considering a special case when N is fixed. Our primary objective is to maximize the profile penalized empirical log-likelihood function for a given N, as specified in Equation (4). In other words, we shall design an EM algorithm to calculate the maximum penalized empirical likelihood estimator of $(\mathsf{\beta },k,\alpha )$ when N is fixed.

In this case, the observed data can be represented as $\mathit{O}=\{(y_1,{\mathit{x}}_1),(y_2,{\mathit{x}}_2),\cdots ,(y_n,{\mathit{x}}_n)\}$, where each $y_i$ is positive. Additionally, the observed data include the counts $(y_{n+1},y_{n+2},\cdots ,y_N)$, all of which are zero. For these individuals not captured, their covariate information is missing and represented as ${\mathit{O}}^{*}=({\mathit{x}}_{n+1}^{*},{\mathit{x}}_{n+2}^{*},\cdots ,{\mathit{x}}_N^{*})$. According to the principle of empirical likelihood, the potential values of the ${\mathit{x}}_i^{*}$’s are drawn from $({\mathit{x}}_1,{\mathit{x}}_2,\cdots ,{\mathit{x}}_n)$, where the associated probabilities are $(p_1,p_2,\cdots ,p_n)$.

The observed and missing data constitute the complete data. The likelihood is as follows:

$$\begin{array}{c}\hfill \prod _{i=1}^n\left\{P(\mathit{X}={\mathit{x}}_i)P(Y=y_i\mid \mathit{X}={\mathit{x}}_i)\right\}\times \prod _{i=n+1}^N\left\{P(\mathit{X}={\mathit{x}}_i^{*})P(Y=0\mid \mathit{X}={\mathit{x}}_i^{*})\right\}.\end{array}$$

Correspondingly, the log-likelihood of $\mathsf{\theta }={({\mathsf{\beta }}^{\top },k,\alpha ,\left\{p_i\right\})}^{\top }$ becomes:

$$\begin{array}{ccc}\hfill {\mathsf{\ell }}_c\left(\mathsf{\theta }\right)& =& \sum _{i=1}^n[log\left\{f(y_i,{\mathit{x}}_i;\mathsf{\beta },k)\right\}+log\left(p_i\right)]\hfill \\ & & +\sum _{j=n+1}^N\sum _{i=1}^{n}I({\mathit{x}}_j^{*}={\mathit{x}}_i)\left[log\left\{f(0,{\mathit{x}}_i;\mathsf{\beta },k)\right\}+log\left(p_i\right)\right].\hfill \end{array}$$

The core of the EM-algorithm is its iterative process, which consists of an expectation step (E-step) followed by a maximization step (M-step) in each iteration. Before these two steps, we use ${\mathsf{\theta }}^{old}=({\mathsf{\beta }}^{old},k^{old},{\alpha }^{old},\left\{p_i^{old}\right\})$ to denote the current value of parameters. In the E-step, we need to compute the expectation of ${\mathsf{\ell }}_c\left(\mathsf{\theta }\right)$ conditional on $\mathit{O}$ and ${\mathsf{\theta }}^{old}$. For this purpose, we calculate the conditional expectation of the indicator $I({\mathit{x}}_j^{*}={\mathit{x}}_i)$, which is equal to:

$$\begin{array}{c}\hfill P(\mathit{X}={\mathit{x}}_i\mid Y=0,\mathsf{\theta }={\mathsf{\theta }}^{old})={\displaystyle \frac{\varphi ({\mathit{x}}_i;{\mathsf{\beta }}^{old},k^{old})p_i^{old}}{{\alpha }^{old}}},\end{array}$$

where ${\alpha }^{old}={\sum }_{i=1}^n\varphi ({\mathit{x}}_i;{\mathsf{\beta }}^{old},k^{old})p_i^{old}$ denotes the current value of $\alpha$. Correspondingly, the conditional expectation of the log-likelihood ${\mathsf{\ell }}_c\left(\mathsf{\theta }\right)$ is equal to:

$$\begin{array}{ccc}\hfill Q\left(\mathsf{\theta }\mid {\mathsf{\theta }}^{old}\right)& =& \sum _{i=1}^n\left[log\left\{f(y_i,{\mathit{x}}_i;\mathsf{\beta })\right\}+u_i^{old}log\left\{f(0,{\mathit{x}}_i;\mathsf{\beta },k)\right\}\right]\hfill \\ & & +\sum _{i=1}^n\left(1+u_i^{old}\right)log\left(p_i\right)\hfill \\ & =:& {\mathsf{\ell }}_1(\mathsf{\beta },k)+{\mathsf{\ell }}_2\left(\left\{p_i\right\}\right).\hfill \end{array}$$

where $u_i^{old}=(N-n)\varphi ({\mathit{x}}_i;{\mathsf{\beta }}^{old},k^{old})p_i^{old}/{\alpha }^{old}$ represents the weight for $i=1,\cdots ,n$.

The M-step consists of maximizing the function $Q(\mathsf{\theta }\mid {\mathsf{\theta }}^{old})$. The separation of parameters $(\mathsf{\beta },k)$ and $p_i$’s makes the maximization procedure much more elegant, which can be implemented using the following steps.

Step 1. 

Update $(\mathsf{\beta },k)$ to $({\mathsf{\beta }}^{new},k^{new})$ by maximizing ${\mathsf{\ell }}_1(\mathsf{\beta },k)$. Given that ${\mathsf{\ell }}_1(\mathsf{\beta },k)$ can be interpreted as a weighted log-likelihood function, we propose that maximizing ${\mathsf{\ell }}_1(\mathsf{\beta },k)$ is analogous to fitting a negative binomial regression model to the observed counts $(y_1,y_2,\cdots ,y_n)$ and the n-dimensional zero vector with covariates $({\mathit{x}}_1,{\mathit{x}}_2,\cdots ,{\mathit{x}}_n,{\mathit{x}}_1,{\mathit{x}}_2,\cdots ,{\mathit{x}}_n)$ and weights $(1,1,\cdots ,1,u_1^{old},u_2^{old},\cdots ,u_n^{old})$. This step can be readily implemented through the glm.nb() function from the MASS package in R.

Step 2. 

Update $p_i$ values by maximizing ${\mathsf{\ell }}_2(p_1,\cdots ,p_n)$ under the positive and sum-to-one constraints. This step yields a closed form, namely, $p_i^{new}=(u_i^{old}+1)/{\sum }_{i=1}^n(u_i^{old}+1)$ for $i=1,\cdots ,n$.

Step 3. 

Update $\alpha$ by calculating ${\alpha }^{new}={\sum }_{i=1}^n\varphi ({\mathit{x}}_i;{\mathsf{\beta }}^{new},k^{new})p_i^{new}$.

The E- and M-steps are repeated until the sequence of $({\mathsf{\beta }}^{new},k^{new},{\alpha }^{new},\left\{p_i^{new}\right\})$ or ${\tilde{\mathsf{\ell }}}_p(N,{\mathsf{\beta }}^{new},k^{new},{\alpha }^{new},\left\{p_i^{new}\right\})$ converges. The EM algorithm outlined above exhibits a desirable property under very general circumstances: the penalized empirical likelihood does not decrease with successive iterations. Given that the penalized empirical log-likelihood is bounded above by zero, the convergence of the sequence $({\mathsf{\beta }}^{new},k^{new},{\alpha }^{new},\left\{p_i^{new}\right\})$ to a local maximum of ${\tilde{\mathsf{\ell }}}_p(N,\mathsf{\beta },k,\alpha ,\left\{p_i\right\})$ is always guaranteed.

To compute the maximum penalized empirical likelihood estimator $({\widehat{N}}_p,{\widehat{\mathsf{\beta }}}_p,{\widehat{k}}_p,{\widehat{\alpha }}_p)$, the aforementioned EM algorithm remains applicable after some modifications. In this scenario, the current parameter is denoted by ${\mathsf{\theta }}^{old}=(N^{old},{\mathsf{\beta }}^{old},k^{old},{\alpha }^{old},\left\{p_i^{old}\right\})$ and the weight is $u_i=(N^{old}-n)\varphi ({\mathit{x}}_i;{\mathsf{\beta }}^{old},k^{old})p_i^{old}/{\alpha }^{old}$ in the E-step. In addition, the M-step incorporates a maximization step for the population size parameter.

Step 4. 

Calculate the updated value $N^{new}$, by maximizing the partial log-likelihood function relavent on N, expressed as $log\left({\displaystyle \genfrac{}{}{0pt}{}{N}{n}}\right)+(N-n)log\left({\alpha }^{new}\right)+f_p\left(N\right)$. This optimization can be efficiently performed using the optimize() function available in the R software (version 4.3.1, https://www.r-project.org/).

The penalized empirical likelihood ratio confidence interval for N is computed by identifying the two zeros of the modified penalized likelihood ratio function:

$$\begin{array}{c}\hfill {\mathsf{\ell }}_m\left(N\right)=2\left\{{\mathsf{\ell }}_p({\widehat{N}}_p,{\widehat{\mathsf{\beta }}}_p,{\widehat{k}}_p,{\widehat{\alpha }}_p)-\underset{(\mathsf{\beta },k,\alpha )}{max}{\mathsf{\ell }}_p(N,\mathsf{\beta },k,\alpha )\right\}-{\chi }_1^2(1-a),\end{array}$$

(6)

where the search for these zeros is conducted within the intervals $[n,{\widehat{N}}_p]$ and $[{\widehat{N}}_p,M]$, and M is a sufficiently large user-specified value ensuring that ${\mathsf{\ell }}_m\left(M\right)>0$. This can be implemented via the uniroot() function available in the R software. In summary, the pseudocodes outlined in Algorithms A1–A3 (Appendix B) offer the procedures for calculating both the maximum penalized empirical likelihood estimator and the corresponding penalized empirical likelihood ratio confidence interval for N.

3. Simulations

To demonstrate the efficiency of penalized empirical likelihood estimators, several simulations are conducted and multiple synthetic datasets are analyzed. In the simulation settings, we fix the abundance $N_0$ at 250, 500, and 1000. We consider two different scenarios for generating the covariate X:

(A)

A binary variable, $X\sim \mathrm{Bi}(1,0.3)$, is used to represent a discrete-valued covariate, as in the case study presented in Section 4.

(B)

Alternatively, a continuous variable, $X\sim U(0,1)$, is considered.

Given $\mathit{X}={(1,X)}^{\top }$, we simulate the count response Y from a negative binomial regression model (1), where the regression coefficient ${\mathsf{\beta }}_0$ is set at ${(-0.5,0.3)}^{\top }$ or ${(0.1,0.3)}^{\top }$ and the dispersion parameter $k_0$ is set at 0.5, 1, or 5.

Under each of the 18 ($3\times 2\times 3$) parameter combinations, we simulate 2000 random samples for each scenario. Subsequently, we calculate the maximum empirical likelihood estimators ($\widehat{N}$ and ${\widehat{N}}_p$) as well as the empirical likelihood ratio interval estimators ($\mathcal{I}$ and ${\mathcal{I}}_p$), where the estimators with the subscript p are derived via the penalized method. When the penalty is applied, we recommend an adaptive tuning parameter value of $C={\left\{2n{(\nu -n)}^2\right\}}^{-1}$, which is proven effective in our numerical studies.

3.1. Evaluation of Point Estimators

We evaluate the performance of two point estimators $\widehat{N}$ and ${\widehat{N}}_p$, by assessing their relative bias in percent (%Rbias) and relative mean squared error (RMSE). For a generic estimator $\check{N}$, the %Rbias and RMSE are as follows:

$$\begin{array}{c}\hfill \%\mathrm{Rbias}={\displaystyle \frac{1}{2000}}\sum _{j=1}^{2000}{\displaystyle \frac{{\check{N}}_j-N_0}{N_0}}\times 100\mathrm{and}\mathrm{RMSE}={\displaystyle \frac{1}{2000}}\sum _{j=1}^{2000}{\displaystyle \frac{{({\check{N}}_j-N_0)}^2}{N_0}},\end{array}$$

respectively, where ${\check{N}}_j$ represents the estimate derived from the jth random sample.

Table 1. Relative biases in percent (%Rbiases) and relative mean squared errors (RMSEs) in simulation studies. The RMSE values are rounded to the nearest integer.

		${\mathsf{\beta }}_{\mathbf{0}}={(-\mathbf{0.5},\mathbf{0.3})}^{\mathsf{\top }}$						${\mathsf{\beta }}_{\mathbf{0}}={(\mathbf{0.1},\mathbf{0.3})}^{\mathsf{\top }}$
		${\mathit{k}}_{\mathbf{0}}=\mathbf{0.5}$		${\mathit{k}}_{\mathbf{0}}=\mathbf{1}$		${\mathit{k}}_{\mathbf{0}}=\mathbf{5}$		${\mathit{k}}_{\mathbf{0}}=\mathbf{0.5}$		${\mathit{k}}_{\mathbf{0}}=\mathbf{1}$		${\mathit{k}}_{\mathbf{0}}=\mathbf{5}$
	${\mathit{N}}_{\mathbf{0}}$	$\widehat{\mathit{N}}$	${\widehat{\mathit{N}}}_{\mathit{p}}$	$\widehat{\mathit{N}}$	${\widehat{\mathit{N}}}_{\mathit{p}}$	$\widehat{\mathit{N}}$	${\widehat{\mathit{N}}}_{\mathit{p}}$	$\widehat{\mathit{N}}$	${\widehat{\mathit{N}}}_{\mathit{p}}$	$\widehat{\mathit{N}}$	${\widehat{\mathit{N}}}_{\mathit{p}}$	$\widehat{\mathit{N}}$	${\widehat{\mathit{N}}}_{\mathit{p}}$
Scenario A
%Rbias	250	82.6	−2.92	34.9	4.24	16.7	10.1	29.8	0.59	5.64	1.86	0.46	0.44
	500	37.8	3.81	14.6	5.53	4.31	4.02	6.65	2.08	1.36	1.14	0.24	0.24
	1000	19.8	7.90	3.27	2.73	1.44	1.44	1.97	1.62	0.55	0.53	−0.06	−0.06
RMSE	250	4006	73	1567	73	330	57	1768	40	136	24	5	4
	500	2558	130	881	114	51	41	271	53	25	23	4	4
	1000	2336	246	123	96	28	28	61	54	17	16	4	4
Scenario B
%Rbias	250	83.14	−1.75	30.38	4.17	15.08	9.02	24.23	1.00	6.07	2.20	0.28	0.25
	500	33.18	4.30	11.01	4.66	3.72	3.57	6.48	2.03	1.01	0.86	0.19	0.19
	1000	15.12	6.95	3.02	2.56	1.19	1.18	1.55	1.24	0.48	0.47	−0.07	−0.07
RMSE	250	4222	73	1228	68	297	56	1143	38	128	24	4	4
	500	2474	121	720	96	34	31	571	48	19	17	4	4
	1000	1168	210	115	84	24	24	55	48	15	15	3	3

reports the simulated %Rbiases and RMSEs of both estimators $\widehat{N}$ and ${\widehat{N}}_p$.

We first examine the simulation results when ${\mathsf{\beta }}_0={(0.1,0.3)}^{\top }$ with the average capture probability of 56%. From the last two columns of Table 1, we see that the %Rbiases and RMSEs of both estimators are comparable when $k_0=5$, where the variance-to-mean ratio (VMR) is as low as 1.3. However, the results differ significantly when $k_0$ decreases to 1 and further to 0.5, with VMRs of 2.3 and 3.5, respectively. As the degree of overdispersion increases, the maximum empirical likelihood estimator $\widehat{N}$ uniformly exhibits larger biases and RMSEs than the maximum penalized empirical likelihood estimator ${\widehat{N}}_p$. For example, in Scenario A when $k_0=0.5$ and $N_0=250$, the relative bias of $\widehat{N}$ is as large as 29.8%, and the RMSE of $\widehat{N}$ is about 44 times (1768/40) the RMSE of ${\widehat{N}}_p$.

Next, we examine the results when ${\mathsf{\beta }}_0={(-0.5,0.3)}^{\top }$. In this scenario, the %Rbiases and RMSEs of both estimators are larger than those when ${\mathsf{\beta }}_0={(0.1,0.3)}^{\top }$. This is expected because the average capture probability reduces from 56% to 40%, making the sample provide less information about the model parameters. As in the aforementioned scenario, similar conclusions can be drawn, with the advantages of ${\widehat{N}}_p$ over $\widehat{N}$ being even clear. All relative biases of ${\widehat{N}}_p$ are smaller than 10%. In contrast, $\widehat{N}$ generally overestimates the population size, with the largest relative bias approaching 83% in Scenario B when $k_0=0.5$ and $N_0=250$. In this case, its RMSE is about 58 times (4222/73) the RMSE of ${\widehat{N}}_p$.

3.2. Evaluation of Interval Estimators

We evaluate and contrast the performance of two empirical likelihood ratio interval estimators: ${\mathcal{I}}_p$ (with penalty) and $\mathcal{I}$ (without penalty). This comparison focuses on their coverage probabilities and interval widths. As discussed in the introduction, the interval estimators may have unbounded upper limits. In the simulation, the interval estimator’s upper limit was set to the minimum of the right endpoint and 100 times the sample size.

Table 2. Simulated coverage probabilities (CPs, unit: %) and average widths (AWs) of interval estimators $\mathcal{I}$ and ${\mathcal{I}}_p$ at the 95% level, along with the proportion of bounded cases (PBC, unit: %) whose upper limits are less than 100 times the observed sample sizes. The AW and PBC values are rounded to the nearest integer.

		${\mathsf{\beta }}_0={(-\mathbf{0.5},\mathbf{0.3})}^{\mathsf{\top }}$						${\mathsf{\beta }}_0={(\mathbf{0.1},\mathbf{0.3})}^{\mathsf{\top }}$
		${\mathit{k}}_{\mathbf{0}}=\mathbf{0.5}$		${\mathit{k}}_{\mathbf{0}}=\mathit{1}$		${\mathit{k}}_{\mathbf{0}}=\mathbf{5}$		${\mathit{k}}_{\mathbf{0}}=\mathbf{0.5}$		${\mathit{k}}_{\mathbf{0}}=\mathbf{1}$		${\mathit{k}}_{\mathbf{0}}=\mathbf{5}$
	${\mathit{N}}_{\mathbf{0}}$	$\mathcal{I}$	${\mathcal{I}}_{\mathit{p}}$	$\mathcal{I}$	${\mathcal{I}}_{\mathit{p}}$	$\mathcal{I}$	${\mathcal{I}}_{\mathit{p}}$	$\mathcal{I}$	${\mathcal{I}}_{\mathit{p}}$	$\mathcal{I}$	${\mathcal{I}}_{\mathit{p}}$	$\mathcal{I}$	${\mathcal{I}}_{\mathit{p}}$
Scenario A
CP	250	94.00	93.90	93.05	93.10	97.75	97.80	93.45	93.60	94.35	94.35	95.60	95.60
	500	94.70	94.70	94.00	94.05	97.15	97.15	93.75	93.90	94.55	94.55	95.25	95.25
	1000	94.60	94.80	95.25	95.25	95.80	95.80	94.65	94.65	95.20	95.20	96.10	96.10
AW	250	742	290	836	363	647	410	776	317	577	330	185	171
	500	1532	745	1294	833	929	745	1104	696	576	498	200	200
	1000	2721	1763	2046	1584	971	909	1426	1152	591	582	258	258
PBC	250	33	100	46	100	67	100	58	100	86	100	100	100
	500	46	100	66	100	90	100	80	100	98	100	100	100
	1000	67	100	89	100	99	100	97	100	100	100	100	100
Scenario B
CP	250	93.85	93.6	93.50	93.55	97.15	97.20	93.95	93.65	93.90	93.95	95.15	95.15
	500	94.85	94.8	94.60	94.60	96.65	96.65	94.30	94.40	94.85	94.90	95.10	95.10
	1000	94.15	94.2	95.35	95.35	96.10	96.10	94.85	94.85	95.65	95.65	95.75	95.75
AW	250	6002	820	5646	973	3654	947	5112	605	2122	458	174	154
	500	9773	1804	6999	1744	2462	1075	4817	1010	887	504	182	181
	1000	12250	3343	4981	2157	896	820	2084	1195	545	536	237	237
PBC	250	34	100	48	100	73	100	60	100	88	100	100	100
	500	48	100	70	100	93	100	83	100	99	100	100	100
	1000	70	100	92	100	100	100	98	100	100	100	100	100

presents the simulation results at the 95% confidence level. Overall, the coverage accuracy of the two interval estimators is similar, with coverage probabilities either matching or slightly exceeding the nominal 95% level. However, the widths of ${\mathcal{I}}_p$ are always shorter than those of $\mathcal{I}$, indicating that the penalized empirical likelihood ratio interval estimator offers greater precision. Specifically, the width of ${\mathcal{I}}_p$ is 12% (605/5112) of that of $\mathcal{I}$ in Scenario B when $N_0=250$, $k_0=0.5$, and ${\mathsf{\beta }}_0={(0.1,0.3)}^{\top }$.

Finally, we explore the possible explanations for the significant differences in the widths of the confidence intervals. For this purpose, Table 2 also reports the proportion of bounded cases, in which the upper limits are less than 100 times the observed sample sizes. In the other cases, the upper limits can be seen unbounded. From Table 2, it is clear that there are many unbounded cases for $\mathcal{I}$. The smaller $N_0$ or k, the more likely this boundary problem occurs. In contrast, the upper limits of ${\mathcal{I}}_p$ are always bounded. This demonstrates that the penalized empirical likelihood method effectively reduces the occurrence of boundary problems, leading to more precise intervals.

4. Case Study

Among the various bear species found in North America, the black bear (Ursus americanus) stands out as the most widely distributed game species. However, their populations have become fragmented, owing to over-harvesting. Understanding black bear demography is crucial for informing effective management and conservation strategies [27].

This study aims to make statistical inferences for the abundance of black bears at the U.S. Army’s Fort Drum Military Installation in northern New York. For this purpose, the capture–recapture experiment was conducted using an array of 38 baited “hair snares” during June and July 2006. The study area and the locations of the 38 hair snares are illustrated in Figure 1.4 in ref. [28]. Each week, for eight consecutive weeks, barbed wire traps were baited and checked for hair samples. Over the 8-week period, a total of 47 black bears were captured at least once, with their sex recorded upon capture. The original data structure also includes the trap locations and whether each bear was captured at each location. Although capture status was recorded for each bear across all 38 hair snares, we treat this as a standard capture–recapture dataset. Specifically, a bear was considered captured if it was caught at least once in any given week. The analyzed data include the number of weeks that bears were captured and the sex of each bear; see

Table 3. Description of the black bear data.

Number of Captures	1	2	3	4	5	6	7
Male	11	7	5	1	2	1	1
Female	8	4	2	1	0	1	3

for the frequencies.

We first examine the estimation results of abundance in the framework of the Poisson regression model. With this model, we find that the maximum empirical likelihood estimate equals 52. This point estimate is smaller than Chao lower bound estimate of 63, indicating that the Poisson regression model may be not adequate to fit this dataset. This inadequacy is further illustrated via the ratio plot (see Figure 1A), indicating a nonconstant relationship between the ratios of frequencies and the number of captures.

Now, we examine the estimation results in the framework of the negative binomial regression model. Applying the maximum empirical likelihood methods with and without penalty, we find that both approaches yield a point estimate of 72 with a standard error of 8.9 for the abundance, significantly higher than the estimate of 52 derived from the Poisson regression model. This difference can be attributed to overdispersion present in the negative binomial regression model, as indicated by an overdispersion parameter estimate of 1.2 with a standard error of 0.2. For completeness, we also calculate the ratio regression estimate proposed in ref. [12], which is 74 with a standard error of 10.2. Our point estimate is close to the ratio regression estimate but more efficient, as indicated by the smaller standard error.

Applying the empirical likelihood method with and without penalty, we obtain significantly different interval estimates of the abundance. At the 95% confidence level, the empirical likelihood ratio interval estimate ($\mathcal{I}$) is [52, 1343], whereas the penalized empirical likelihood ratio interval estimate (${\mathcal{I}}_p$) is [52, 201]. Clearly, $\mathcal{I}$ exhibits a significantly larger upper limit in comparison to ${\mathcal{I}}_p$. In order to understand this discrepancy, we plot the profile empirical log-likelihood ratio functions with and without penalty of the abundance (see Figure 1B). The two functions are very close when the abundance N is less than 100 but diverge when $N>100$. The empirical log-likelihood ratio function increases gradually without penalty, flattening as N increases. In contrast, the penalized log-likelihood ratio function exhibits a rapid ascent with increasing N, thereby enhancing the precision of the confidence interval. This demonstrates how the penalized empirical likelihood method effectively mitigates the boundary problem caused by overdispersion.

As suggested by one peer reviewer, we also applied the classical Lincoln–Petersen method [1,2] to the black bear data. For performing the Lincoln–Petersen method, the eight weeks of data were divided into two samples: the initial weeks constituted sample 1, while the remaining weeks formed sample 2. The Lincoln–Petersen estimator is calculated as $n_1{n}_2/n_{12}$, where $n_1$, $n_2$, and $n_{12}$ represent the number of individuals in sample 1, sample 2, and the overlap between the two samples, respectively. The variance of this estimator is calculated as:

$${\displaystyle \frac{(n_1+1)(n_2+1)(n_1-n_{12})(n_2-n_{12})}{{(n_{12}+1)}^2(n_{12}+2)}}.$$

Table 4. Lincoln–Petersen estimates and standard errors (SEs).

Sample 1	Sample 2	${\mathit{n}}_{\mathbf{1}}$	${\mathit{n}}_{\mathbf{2}}$	${\mathit{n}}_{\mathbf{12}}$	Estimate	SE
week 1	weeks 2–8	9	45	7	58	7.8
weeks 1–2	weeks 3–8	14	43	10	60	7.7
weeks 1–3	weeks 4–8	27	38	18	57	5.2
weeks 1–4	weeks 5–8	33	38	24	52	3.2
weeks 1–5	weeks 6–8	38	31	22	54	3.8
weeks 1–6	weeks 7-8	40	31	22	54	3.8
weeks 1–7	week 8	45	11	9	55	6.0

presents the Lincoln–Petersen estimates along with their standard errors. It is clear that all estimates are consistently lower than the maximum penalized empirical likelihood estimate of 72. This discrepancy is likely due to the assumption of sample independence inherent in the Lincoln–Petersen method, which may not hold true in this context. Specifically, bears captured in the first sample may be more likely to be recaptured due to the lure, leading to a positive correlation between the two samples and an overestimation of capture probability. Additionally, the standard errors for the Lincoln–Petersen estimates are uniformly smaller than the standard error of 8.9 obtained from our proposed method. This difference might arise because the Lincoln–Petersen estimator assumes equal capture probability across all individuals, failing to account for heterogeneity, such as sex-based variation. In contrast, our proposed method explicitly incorporates sex as a factor, addressing this limitation.

Given that the true population size is unknown, we aim to validate the reliability of our method’s conclusions by assessing whether the negative binomial regression model accurately fits the observed frequency data. To this end, we have included a goodness-of-fit test in the revision. Specifically, for the black bear data, we define the ${\chi }^2$ test statistic as:

$${\chi }^2={\displaystyle \frac{{(e_1-m_1)}^2}{e_1}}+{\displaystyle \frac{{(e_2-m_2)}^2}{e_2}}+{\displaystyle \frac{{(e_3-m_3)}^2}{e_3}}+{\displaystyle \frac{{(e_4-m_4)}^2}{e_4}}+{\displaystyle \frac{{(e_5-m_5)}^2}{e_5}},$$

where $m_1=19$, $m_2=11$, $m_3=7$, $m_4=2$, and $m_5=8$ represent the observed frequencies of captures occurring once, twice, three times, four times, and at least five times, respectively. The expected frequencies are computed as:

$$e_k=n{\displaystyle \frac{{\sum }_{i=1}^n{\widehat{p}}_{i}f(k,{\mathit{x}}_i;{\widehat{\mathsf{\beta }}}_p,{\widehat{k}}_p)}{{\sum }_{i=1}^n{\widehat{p}}_i\{1-\varphi ({\mathit{x}}_i;{\widehat{\mathsf{\beta }}}_p,{\widehat{k}}_p)\}}},k=1,2,3,4,$$

with $e_5=n-{\sum }_{k=1}^4{e}_k$, where the ${\widehat{p}}_i$’s are those convergence values of the $p_i$’s obtained from the EM algorithm. The ${\chi }^2$ test statistic, calibrated with a ${\chi }_4^2$ distribution, yields a value of 1.66 with a p-value of 0.8. This result suggests that the negative binomial regression model fits the black bear data well, thereby supporting the reliability of our estimation results with this model.

5. Conclusions and Discussion

The use of the negative binomial regression model is prevalent to address heterogeneity and dispersion related to capture–recapture frequency data. As noted in refs. [10,11], fitting this model often leads to identification and boundary problems for the dispersion parameter and then yields unbounded population size estimates. To address this boundary problem, we proposed imposing a half-normal prior on the population size or equivalently decreasing the empirical log-likelihood function for large population sizes by adding a suitable penalty function. This penalized technique could improve robustness of the maximum penalized empirical likelihood estimator, ensuring the derivation of a consistently bounded interval estimator. This penalized empirical likelihood approach constitutes a significant contribution of our paper. Additionally, we introduced an efficient EM algorithm to maximize the penalized empirical likelihood function. Unlike the classical Newton-type optimization method, the EM algorithm guarantees local convergence of the numerical procedure and yields stable estimates of population size. Compared to estimators without penalty, the proposed maximum penalized empirical likelihood estimator exhibits a higher efficiency and the penalized empirical likelihood ratio interval estimator is more precise.

This paper assumes that the capture–recapture experiment is conducted at a single site. However, in ecological studies, experiments are often conducted at multiple sites or follow a spatial pattern [29]. For example, the case study in Section 4 belongs to this case. To accommodate these general application scenarios, the count data can be considered as copies of $(Y_{\left(1\right)},Y_{\left(2\right)},\cdots ,Y_{\left(J\right)})$, where, for $j=1,2,\cdots ,J$, $Y_{\left(j\right)}$ represents the number of times a generic individual was captured at the jth site. In such cases, the multivariate mixed negative binomial regression model used in ref. [30] can be adopted to fit the count data. Specifically, the conditional probability model (1) becomes:

$$\begin{array}{cc}& P(Y_{\left(1\right)}=y_{\left(1\right)},Y_{\left(2\right)}=y_{\left(2\right)},\cdots ,Y_{\left(J\right)}=y_{\left(J\right)}\mid \mathit{X}=\mathit{x})\hfill \\ \hfill =& {\int }_0^{\infty }\cdots {\int }_0^{\infty }\prod _{j=1}^J{\displaystyle \frac{\Gamma (y_{\left(j\right)}+k_j)}{\Gamma (y_{\left(j\right)}+1)\Gamma \left(k_j\right)}}{\left\{{\displaystyle \frac{k_j}{k_j+{\lambda }_j\mu (\mathit{x};\mathsf{\beta })}}\right\}}^{k_j}{\left\{{\displaystyle \frac{{\lambda }_j\mu (\mathit{x};\mathsf{\beta })}{k_j+{\lambda }_j\mu (\mathit{x};\mathsf{\beta })}}\right\}}^{y_{\left(j\right)}}\hfill \\ & h({\lambda }_1,\cdots ,{\lambda }_J;\gamma )d{\lambda }_1\cdots d{\lambda }_J,\hfill \end{array}$$

where $({\lambda }_1,\cdots ,{\lambda }_J)$ describes the sites’ random effects and $h({\lambda }_1,\cdots ,{\lambda }_J;\gamma )$ describes the spatial variation of capture intensities across sites. Following the approach used in ref. [30], one can opt for a gamma distribution as the mixing distribution, with a probability density function described by:

$$h({\lambda }_1,\cdots ,{\lambda }_J;\gamma )={\displaystyle \frac{{\gamma }^{\gamma }}{\Gamma \left(\gamma \right)}}exp(-\gamma \lambda ){\lambda }^{\gamma -1}I({\lambda }_1=\cdots ={\lambda }_J=\lambda ),$$

where $\gamma >0$. When spatial auto-correlation arises, indicating that random effects are partially correlated across different sites, $(log\left({\lambda }_1\right),\cdots ,log\left({\lambda }_J\right))$ can be modeled by a Gaussian random field. The covariance matrix of this field can be determined by flexible spatial covariance functions [31].

Throughout the paper, the conditional mean of the count data are assumed to have a log-linear relationship with individual covariates, as depicted in model (1). In practice, the relationship between the conditional mean and covariates may be more complex. To adapt our approach for practical application scenarios, one can employ the negative binomial additive model with flexible smoothing functions for each covariate [32]. In addition, applying the proposed methods to a more generalized Poisson mixture model, where the Tweedie distribution is used as the mixing distribution [33], is a straightforward yet challenging task.