A Statistical Model for Count Data Analysis and Population Size Estimation: Introducing a Mixed Poisson–Lindley Distribution and Its Zero Truncation

Abstract

Count data consists of both observed and unobserved events. The analysis of count data often encounters overdispersion, where traditional Poisson models may not be adequate. In this paper, we introduce a tractable one-parameter mixed Poisson distribution, which combines the Poisson distribution with the improved second-degree Lindley distribution. This distribution, called the Poisson-improved second-degree Lindley distribution, is capable of effectively modeling standard count data with overdispersion. However, if the frequency of the unobserved events is unknown, the proposed distribution cannot be directly used to describe the events. To address this limitation, we propose a modification by truncating the distribution to zero. This results in a tractable zero-truncated distribution that encompasses all types of dispersions. Due to the unknown frequency of unobserved events, the population size as a whole becomes unknown and requires estimation. To estimate the population size, we develop a Horvitz–Thompson-like estimator utilizing truncated distribution. Both the untruncated and truncated distributions exhibit desirable statistical properties. The estimators for both distributions, as well as the population size, are asymptotically unbiased and consistent. The current study demonstrates that both the truncated and untruncated distributions adequately explain the considered medical datasets, which are the number of dicentric chromosomes after being exposed to different doses of radiation and the number of positive Salmonella. Moreover, the proposed population size estimator yields reliable estimates.

1. Introduction

Unobserved events in count data are events that were not recorded. For example, the unobserved events in insurance claims refer to the events where the policyholders do not claim. In some cases, even these events are not identified. For example, the number of times a motorist gets stopped by the police. The number of motorists that did not get stopped cannot be identified, which is known as unobserved events. Modeling count data by examining the observed events or both the observed and unobserved events is typical in statistical modeling. This study would like to explore this area by introducing a mixed Poisson–Lindley distribution and extending it to its zero-truncated version, as well as population size estimation.

The Lindley distribution was originally introduced in Bayesian analysis as a mixture of an exponential distribution and a gamma distribution [1]. Ghitany et al. [2] conducted a comprehensive examination of its statistical characteristics. Subsequent research has expanded the Lindley distribution into both two-parameter [3,4,5,6,7,8,9,10,11] and three-parameter variants [12,13,14], each offering innovative enhancements and broader applications. These distributions have been particularly useful in the fields of survival analysis and reliability assessments.

Sankaran [15] pioneered the use of the Lindley distribution as a mixing distribution with the Poisson distribution, thus creating the Poisson–Lindley (PL) distribution. This model, and the methods to estimate its parameters, were studied extensively in later research [16]. Various mixed PL distributions have been proposed as alternatives to the traditional Poisson and negative binomial models for fitting count data [17,18,19]. These new distributions share the overdispersion characteristic common to mixed Poisson distributions [20]. Other notable mixed Poisson models include the Poisson Inverse Pareto distribution [21] and the Poisson-transmuted record-type exponential distribution [22], among others.

When count data exclusively contain positive numbers, zero truncation is a common technique used to adjust distributions accordingly. Examples include the zero-truncated Poisson [23], zero-truncated negative binomial [24], and zero-truncated PL [25] distributions. In fields such as criminology, accurately estimating the size of a population when the frequency of non-events is not observable is a significant challenge [26,27,28,29,30,31,32,33]. Rossmo and Routledge [30] highlighted the necessity of understanding the size of a criminal population to inform the creation of effective laws and policies. Such estimations are still relatively rare in medical research but are equally needed.

The current study aims to introduce a new, practical mixed PL distribution by combining an improved second-degree Lindley (ISDL) distribution [7] with a Poisson variable, resulting in what we name the Poisson-improved second-degree Lindley (PISDL) distribution. This constitutes the first aim of this research. Given the ISDL distribution’s superior modeling performance [7], we anticipate that the PISDL distribution may outperform the original PL distribution. To accommodate strictly positive data, we also introduce a zero-truncated version of the PISDL (ZTPISDL) distribution, which represents the second aim of this research. Furthermore, we propose an innovative estimator for population size in relation to the ZTPISDL distribution, detailed in Section 5.3, marking the third and most original objective of the paper. Understanding the population size is critical for developing comprehensive policies. This study uses data from epidemiology and cytogenetics to demonstrate the applications of the proposed distributions and estimator.

This paper is structured as follows. Section 1 outlines the study’s objectives, building upon previous research on mixed PL distributions, their truncations, and population size estimations. Section 2 presents the PISDL distribution, its statistical properties, and estimation techniques. Section 3 introduces and investigates the zero-truncated PISDL distribution. Section 4 describes the development of a new population size estimator, assuming a PISDL distribution, accompanied by simulation studies in Section 2, Section 3 and Section 4 to evaluate estimator performance. Section 5 applies the proposed models to medical dataset analyses. Finally, Section 6 concludes with a discussion of the implications, limitations, and future research directions of this study.

2. Poisson-Improved Second-Degree Lindley Distribution

2.1. Probability Mass Function of the PISDL Distribution

Definition 1.

A random variable $X$ is said to follow a Poisson-improved second-degree Lindley (PISDL) distribution with parameter θ if it obeys $X|\theta ~Poisson\left(\theta \right)$ and$\theta ~ISDL\left(\lambda \right)$, where $\theta , \lambda >0$.

Theorem 1.

Let $X$ be a random variable that follows PISDL with parameter $\lambda$; then, the probability mass function (pmf) of $X$ is given as

$$f\left(x\right)=\frac{{\lambda }^3}{{\lambda }^2+2\lambda +2}\left[\frac{x\left(x+5\right)+2\lambda \left(x+2\right)+{\lambda }^2+5}{{\left(\lambda +1\right)}^{x+3}}\right]; x=0,1,2,\dots , \lambda >0.$$

(1)

Proof. 

Using Definition 1, the pmf of $X|\theta$ is

$$f\left(x|\theta \right)=\frac{{\theta }^x\exp \left(-\theta \right)}{x!}; x=0,1,2,\dots , \theta >0.$$

The probability density function of $\theta |\lambda$ is

$$g\left(\theta \right)=\frac{{\lambda }^3}{{\lambda }^2+2\lambda +2}{\left(1+\theta \right)}^2\exp \left(-\lambda \theta \right); \lambda >0,\theta >0.$$

The resulting marginal distribution $f\left(x\right)$for the PISDL distribution with parameter $\lambda$ is obtained as follows:

$$\begin{array}{ll}f\left(x\right)& ={\displaystyle \underset{0}{\overset{\infty }{\int }}f\left(x|\theta \right)g\left(\theta \right)d\theta }\\ & =\frac{{\lambda }^3}{\left({\lambda }^2+2\lambda +2\right)x!}{\displaystyle \underset{0}{\overset{\infty }{\int }}{\theta }^x{\left(1+\theta \right)}^2\exp \left[-\theta \left(\lambda +1\right)\right]d\theta }\\ & =\frac{{\lambda }^3}{{\lambda }^2+2\lambda +2}\left[\frac{x\left(x+5\right)+2\lambda \left(x+2\right)+{\lambda }^2+5}{{\left(\lambda +1\right)}^{x+3}}\right].\end{array}$$

 □

Figure 1 shows the pmf plot for the PISDL distribution. In Figure 1, the distribution is skewed to the right, unimodal, and decreasing, which is further supported by the decreasing ratio of probability given as

$$\frac{f\left(x+1\right)}{f\left(x\right)}=\frac{1}{\lambda +1}\left[1+\frac{2\left(x+\lambda +3\right)}{x\left(x+5\right)+2\lambda \left(x+2\right)+{\lambda }^2+5}\right].$$

It is worth noting that the PISDL distribution is actually a three-component mixture distribution that can be written as $f\left(x\right)=p_1{f}_1\left(x\right)+p_2{f}_2\left(x\right)+p_3{f}_3\left(x\right)$, where $f_i\left(x\right)$ is the pmf of the negative binomial distribution with a parameter number of successes i and proportion $\lambda /\left(\lambda +1\right)$. When $i=1$, $f_i\left(x\right)$ is the pmf of the geometric distribution, which is a special case for the negative binomial distribution. The formulae for $p_i$ and $f_i\left(x\right)$ for $i=1,2,3$ are given as

$$p_1=\frac{{\lambda }^2}{{\lambda }^2+2\lambda +2}; p_2=\frac{2\lambda }{{\lambda }^2+2\lambda +2}; p_3=1-p_1-p_2;$$

$$f_1\left(x\right)=\frac{\lambda }{{\left(\lambda +1\right)}^{x+1}}; f_2\left(x\right)=\frac{{\lambda }^2\left(x+1\right)}{{\left(\lambda +1\right)}^{x+2}}; f_3\left(x\right)=\frac{{\lambda }^3\left(x+1\right)\left(x+2\right)}{2{\left(\lambda +1\right)}^{x+3}}.$$

Even though the PISDL distribution is a three-component mixture negative binomial distribution, the existence of three modes cannot be seen in any of the plots in Figure 1 for the selected values of $\lambda$. This insinuates that the three modes, which come from the three sub-populations, must be located very close to each other. As mentioned in [34], if the modes of the sub-populations are located very close to each other, then the population will have a single mode. As such, if the existence of the modes of the sub-populations, each with very close mode values, can be certain, then this distribution can be considered as one of the candidates for model fittings.

The cumulative distribution function (cdf) for the PISDL distribution is given in Equation (2) and visualized in Figure 2. In Figure 2, it is clear that the PISDL distribution has a valid cdf since $F\left(x\right)\to 1$ as $x\to \infty$.

$$F\left(x\right)=1-\frac{{\lambda }^4+2{\lambda }^3\left(x+3\right)+{\lambda }^2\left(x^2+7x+13\right)+2\lambda \left(x+4\right)+2}{\left({\lambda }^2+2\lambda +2\right){\left(\lambda +1\right)}^{x+3}}.$$

(2)

Based on Equation (2), the survival function $S\left(x\right)$ can be obtained and given as

$$S\left(x\right)=\frac{{\lambda }^4+2{\lambda }^3\left(x+3\right)+{\lambda }^2\left(x^2+7x+13\right)+2\lambda \left(x+4\right)+2}{\left({\lambda }^2+2\lambda +2\right){\left(\lambda +1\right)}^{x+3}}.$$

(3)

The plot for the survival function of the PISDL distribution is given in Figure 3.

The hazard rate function $h\left(x\right)$ is defined by taking the ratio of the pmf to the survival function, i.e., $h\left(x\right)=f\left(x\right)/S\left(x\right)$, and is given as

$$h\left(x\right)=\frac{{\lambda }^3\left[x\left(x+5\right)+2\lambda \left(x+2\right)+{\lambda }^2+5\right]}{{\lambda }^4+2{\lambda }^3\left(x+3\right)+{\lambda }^2\left(x^2+7x+13\right)+2\lambda \left(x+4\right)+2}.$$

(4)

The hazard rate function plot is given in Figure 4. In Figure 4, it can be noted that the hazard rate functions show an increasing pattern with a limiting value of $\lambda$, meaning $\underset{x\to \infty }{\lim }h\left(x\right)=\lambda$.

2.2. Some Statistical Properties of the PISDL Distribution

The $k^{th}$ moment in the origin of the PISDL distribution can be written by the following generic expression:

$$E\left(X^k\right)={\displaystyle \sum _{x=0}^{\infty }{x}^{k}f\left(x\right)}.$$

(5)

In particular, the first two moments in the origin using Equation (5) are obtained, respectively, as

$$E\left(X\right)=\frac{{\lambda }^2+4\lambda +6}{\lambda \left({\lambda }^2+2\lambda +2\right)}, E\left(X^2\right)=\frac{\left(\lambda +2\right)\left({\lambda }^2+4\lambda +6\right)+4\left(\lambda +3\right)}{{\lambda }^2\left({\lambda }^2+2\lambda +2\right)}.$$

Hence, the index of dispersion (${IOD}_X$) can be written as

$$\begin{array}{cc}\hfill IO{D}_X& =\frac{Var\left(X\right)}{E\left(X\right)}=\frac{E\left(X^2\right)-{\left[E\left(X\right)\right]}^2}{E\left(X\right)}=\frac{E\left(X^2\right)}{E\left(X\right)}-E\left(X\right)\hfill \\ & =1+\frac{1}{\lambda }\left(\frac{{\lambda }^4+2\left(4{\lambda }^3+12{\lambda }^2+12\lambda +6\right)}{{\lambda }^4+6{\lambda }^3+16{\lambda }^2+20\lambda +12}\right).\hfill \end{array}$$

(6)

Since the ${IOD}_X>1$, the PISDL distribution is overdispersed for all $\lambda$.

The mode of the PISDL distribution can be obtained by maximizing the log pmf of the PISDL distribution or equivalently by solving the quadratic equation $A{x_m}^2+B{x}_m+C=0$, where $x_m\ge 0$ is the mode of the distribution, where $A=\ln \left(\lambda +1\right)$, $B=\left(2\lambda +5\right)\ln \left(\lambda +1\right)-2$, and $C=\left(5+\lambda \left(\lambda +4\right)\right)\ln \left(\lambda +1\right)-2\lambda -5$. As a result, the solution for the equation is

$$x_m=\frac{2-\left(5+2\lambda \right)\ln \left(1+\lambda \right)+\sqrt{4+\left(5+4\lambda \right){\left[\ln \left(1+\lambda \right)\right]}^2}}{2\ln \left(1+\lambda \right)}.$$

It can be shown from Figure 1 that $x_m>0,$ for $0<\lambda <1$and when $\lambda \ge 1$, then $x_m=0.$ The moment, the probability and the cumulant generating functions are given, respectively, is

$$M_X\left(t\right)=\frac{{\lambda }^3}{{\lambda }^2+2\lambda +2}\left(s+2s^2+2s^3\right); s=s\left(t\right)={\left[\lambda +1-\exp \left(t\right)\right]}^{-1},$$

$$G_X\left(t\right)=\frac{{\lambda }^3}{{\lambda }^2+2\lambda +2}\left(r+2r^2+2r^3\right); r=s\left(\ln t\right)={\left(\lambda +1-t\right)}^{-1},$$

$$C_X\left(t\right)=\ln M_X\left(t\right)=3\ln \lambda -\ln \left({\lambda }^2+2\lambda +2\right)+\ln \left(s+2s^2+2s^3\right).$$

2.3. Parameter Estimation of the PISDL Distribution

The parameter of the PISDL distribution needs to be estimated before modeling real-world datasets. Here, parameter estimation is based on the two commonly used estimation methods, which are the methods of moments and maximum likelihood.

2.3.1. Method of Moments Estimator

The method of moments estimator can be obtained by equating the sample mean to the population mean. Therefore, the moment estimator of $\lambda$, hereby denoted as $\tilde{\lambda },$can be obtained by solving the following equation:

$$\overline{x}=\frac{{\tilde{\lambda }}^2+4\tilde{\lambda }+6}{\tilde{\lambda }\left({\tilde{\lambda }}^2+2\tilde{\lambda }+2\right)},$$

(7)

or equivalently solving the cubic equation $A{\tilde{\lambda }}^3+B{\tilde{\lambda }}^2+C\tilde{\lambda }+D=0$, where $A=\overline{x}$, $B=2\overline{x}-1$, $C=2\overline{x}-4$, and $D=-6$.

2.3.2. Maximum Likelihood Estimator

Note that the log-likelihood function can be written as$l=\ln L={\sum }_{i=1}^n\ln f\left(x_i\right)={\sum }_{x=0}^{\infty }{n}_x\ln \left[f\left(x\right)\right]$, where $n_x$ is the frequency of x-valued data. Hence, l is given as

$$\begin{array}{cc}\hfill l=& 3n\ln \lambda -n\ln \left({\lambda }^2+2\lambda +2\right)-n\left(\overline{x}+3\right)\ln \left(\lambda +1\right)\hfill \\ & +{\displaystyle \sum _{x=0}^{\infty }{n}_x\ln \left[x\left(x+5\right)+2\lambda \left(x+2\right)+{\lambda }^2+5\right]}.\hfill \end{array}$$

By differentiating l with respect to $\lambda$, we obtain

$$\frac{dl}{d\lambda }=\frac{3n}{\lambda }-\frac{2n\left(\lambda +1\right)}{{\lambda }^2+2\lambda +2}-\frac{n\left(\overline{x}+3\right)}{\lambda +1}+{\displaystyle \sum _{x=0}^{\infty }{n}_x\frac{2\left(\lambda +x+2\right)}{x\left(x+5\right)+2\lambda \left(x+2\right)+{\lambda }^2+5}},$$

(8)

in which equating it to zero and solving it yields the maximum likelihood estimator (MLE) for $\lambda$, which is abbreviated as$\widehat{\lambda }$. Equivalently, one can directly maximize the l to obtain a similar result. To estimate the variance of $\widehat{\lambda }$, Fisher’s information about $\lambda$ and $I\left(\lambda \right)$ needs to be obtained, and it is given as

$$\begin{array}{cc}\hfill I\left(\lambda \right)=& -\frac{4{\lambda }^6+7{\lambda }^5-7{\lambda }^4-40{\lambda }^3-52{\lambda }^2-36\lambda -12}{{\lambda }^2{\left(\lambda +1\right)}^2{\left({\lambda }^2+2\lambda +2\right)}^2}\hfill \\ & +\frac{4{\lambda }^3}{{\left(\lambda +1\right)}^3\left({\lambda }^2+2\lambda +2\right)}{\displaystyle \sum _{x=0}^{\infty }\frac{{\left(\lambda +x+2\right)}^2}{\left[x\left(x+5\right)+2\lambda \left(x+2\right)+{\lambda }^2+5\right]{\left(\lambda +1\right)}^x}.}\hfill \end{array}$$

(9)

The summation term cannot be written in a closed form, but one may use several Lerch Transcendent [35] functions for it. However, it is only practical to leave the summation term as is. So, the variance of $\widehat{\lambda }$ can be written as $I^{-1}\left(\lambda \right)/n$ and, subsequently, the $\left(1-\alpha \right)100\%$ confidence interval can be written as $\widehat{\lambda }\mp z_{\alpha }{I}^{-1/2}\left(\lambda \right)/\sqrt{n}$. If $\alpha =0.025$, and then the 95% confidence interval will be obtained. The summation term reduces considerably to a simple form when $\widehat{\lambda }$ is substituted, thus resulting in a constant. This will eventually give an estimated variance of $\widehat{\lambda }$. For example, for $n=100$ and when $\widehat{\lambda }$ = 1.8793, the summation equates to 1.4284, which then results in the estimated variance of 0.0145.

2.4. Simulation Study

A simulation study is conducted to assess the performance of the two earlier estimators in estimating the parameter of the PISDL distribution. The algorithm for the simulation study is as follows:

Step 1:

Generate $N=1000, 2000, \dots , \mathrm{10,000}$random data that follows the PISDL distribution with $\lambda =0.5$.

Step 2:

Obtain the estimated $\lambda$ using MLE and moment estimator.

Step 3:

Repeat Steps 1–2 for a total of 2000 iterations and obtain the estimates.

Step 4:

Calculate the mean absolute deviation, MAD, and the mean squared error values, MSEs, given, respectively, as $MAD={\sum }_{i=1}^{2000}\left|\check{\lambda }-\lambda \right|/2000$ and $MSE={\sum }_{i=1}^{2000}{\left(\check{\lambda }-\lambda \right)}^2/2000$, where $\check{\lambda }$ can be the MLE or moment estimator for $\lambda$.

Step 5:

Repeat Steps 3–4 for $\lambda =2.0, 5.0$.

Generally, if the MAD and MSE approach zero as N increases, then the estimate is asymptotically unbiased and consistent. For the simulation study, R software version 3.0.2 is used, and the estimated parameter using MLE is obtained using the optim command to optimize the log-likelihood value. The results of the simulation study are presented in Figure 5. In Figure 5, for any value of $\lambda$, as N increases, the MAD and the MSE of the MLE and the moment estimator decrease, suggesting that both estimates are asymptotically unbiased and consistent.

3. Zero-Truncated Poisson-Improved Second-Degree Lindley Distribution

Truncation is a widely used trait in real-world situations in a variety of domains, including industry [23,24,25], medicine [23,24], and many more. The progression of a disease that is not an increasing function but will stabilize after a certain period is an example of truncation. Therefore, a flexible truncated count data distribution is introduced by truncating the PISDL distribution at zero, yielding a zero-truncated Poisson-improved second-degree Lindley distribution (ZTPISDL). It was observed that the PISDL and the PL distributions are equally competent based on the two datasets considered in Section 5. Therefore, it is expected that the ZTPISDL distribution to be as competent as the zero-truncated PL (ZTPL) distribution. The development and the statistical properties of the ZTPISDL distribution are discussed in the following sections.

3.1. Probability Mass Function of the ZTPISDL Distribution

Definition 2.

A random variable Y is said to follow a ZTPISDL distribution if it obeys $X|\lambda ~PISDL\left(\lambda \right)$ and $p\left(y\right)=f\left(y\right)/\left[1-f\left(0\right)\right]$, where p(y) is the pmf of Y.

Theorem 2.

Let Y be a random variable that follows the ZTPISDL with parameter $\lambda$, and then the pmf of Y is given as

$$p\left(y\right)=\frac{{\lambda }^3}{{\lambda }^4+6{\lambda }^3+13{\lambda }^2+8\lambda +2}\left[\frac{y\left(y+5\right)+2\lambda \left(y+2\right)+{\lambda }^2+5}{{\left(\lambda +1\right)}^y}\right], y=1,2,3,\dots$$

(10)

Proof. 

Using Definition 2, the pmf of Y is obtained. □

Figure 6 shows the pmf plot for the ZTPISDL distribution that follows similar shapes as the pmf plot for the PISDL distribution, which is skewed to the right, unimodal, and decreasing. Using the log of $p\left(y\right)$ defined in (10), the mode $y_{mod}$ of the distribution can be obtained by solving the following formula:

$$y_{\mathrm{mod}}=\frac{-D+\sqrt{D^2-4AE}}{2A},$$

where $A=\ln \left(\lambda +1\right), D=A\left(2\lambda +5\right)-2$, and $E=A\left({\lambda }^2+4\lambda +5\right)-2\lambda -5$.

3.2. Some Statistical Properties of the ZTPISDL Distribution

If $Y~ZTPISDL$ for $y=1, 2, 3, \dots$, then the $k^{th}$moment for Y can be easily obtained because it satisfies $E\left(Y^k\right)=E\left(X^k\right)/\left[1-f\left(0\right)\right]$. The first two moments in the origin are obtained and can be, respectively, written as

$$\begin{array}{l}E\left(Y\right)=\frac{{\left(\lambda +1\right)}^3\left({\lambda }^2+4\lambda +6\right)}{\lambda \left({\lambda }^4+6{\lambda }^3+13{\lambda }^2+8\lambda +2\right)},\\ E\left(Y^2\right)=\frac{{\left(1+\lambda \right)}^3\left(\left(\lambda +2\right)\left({\lambda }^2+4\lambda +6\right)+4\left(3+\lambda \right)\right)}{{\lambda }^2\left(2+\lambda \left(1+\lambda \right)\left(8+\lambda \left(5+\lambda \right)\right)\right)}.\end{array}$$

Using a similar approach as the ${IOD}_X$for X, the index of dispersion for Y (${IOD}_Y$) can be written as

$$IO{D}_Y=\frac{Var\left(Y\right)}{E\left(Y\right)}=\frac{1}{\lambda }-\frac{2\left(2+\lambda \right)}{6+\lambda \left(4+\lambda \right)}+\frac{4+2\lambda \left(8+\lambda \left(5+\lambda \right)\right)}{2+\lambda \left(1+\lambda \right)\left(8+\lambda \left(5+\lambda \right)\right)}.$$

Based on the ${IOD}_Y$, the ZTPISDL distribution is underdispersed (overdispersed) when $\lambda >\left(<\right)1.51494$. The ZTPISDL distribution is only equidispersed when $\lambda =1.51494$. The recurrence probability for Y is similar for X, except the x is substituted with y because

$$\frac{p\left(y+1\right)}{p\left(y\right)}=\frac{f\left(y+1\right)/\left[1-f\left(0\right)\right]}{f\left(y\right)/\left[1-f\left(0\right)\right]}=\frac{f\left(y+1\right)}{f\left(y\right)}.$$

The generating functions for both ZTPISDL and PISDL distributions can be related since their pmfs are related as well. Therefore, their relationships are given, respectively, below so that the moment generating, the cumulant generating, and the probability generating functions for Y can be worked out as follows:

$$M_Y\left(t\right)=\frac{M_X\left(t\right)-f\left(0\right)}{1-f\left(0\right)},$$

$$C_Y\left(t\right)=\ln \left[M_Y\left(t\right)\right]=\ln \left[M_X\left(t\right)-f\left(0\right)\right]-\ln \left[1-f\left(0\right)\right],$$

$$G_Y\left(t\right)=M_Y\left(\ln t\right)=\frac{M_X\left(\ln t\right)-f\left(0\right)}{1-f\left(0\right)}=\frac{G_X\left(t\right)-f\left(0\right)}{1-f\left(0\right)}.$$

3.3. Parameter Estimation of the ZTPISDL Distribution

The parameter of the ZTPISDL distribution can be estimated using the moment and MLE techniques. The moment estimator of $\lambda$, hereby denoted as $\tilde{\lambda }$, can be obtained by solving the following equation:

$$\overline{y}=\frac{{\left(\tilde{\lambda }+1\right)}^3\left({\tilde{\lambda }}^2+4\tilde{\lambda }+6\right)}{\tilde{\lambda }\left({\tilde{\lambda }}^4+6{\tilde{\lambda }}^3+13{\tilde{\lambda }}^2+8\tilde{\lambda }+2\right)},$$

(11)

or equivalently solving the quintic equation $A{\tilde{\lambda }}^5+B{\tilde{\lambda }}^4+C{\tilde{\lambda }}^3+D{\tilde{\lambda }}^2+E\tilde{\lambda }+F=0$, such that $A=1-\overline{y}$, $B=7-6\overline{y}$, $C=21-13\overline{y}$, $D=31-8\overline{y}$, $E=22-2\overline{y}$, and $F=6$. For MLE, the log-likelihood function l is given as

$$\begin{array}{cc}\hfill l=& 3n\ln \lambda -n\ln \left({\lambda }^4+6{\lambda }^3+13{\lambda }^2+8\lambda +2\right)-n\overline{y}\ln \left(\lambda +1\right)\hfill \\ & +{\displaystyle \sum _{y=1}^{\infty }{n}_y\ln \left[y\left(y+5\right)+2\lambda \left(y+2\right)+{\lambda }^2+5\right]}.\hfill \end{array}$$

By differentiating l with respect to $\lambda$, we obtain

$$\frac{dl}{d\lambda }=\frac{3n}{\lambda }-\frac{n\left(4{\lambda }^3+18{\lambda }^2+26\lambda +8\right)}{{\lambda }^4+6{\lambda }^3+13{\lambda }^2+8\lambda +2}-\frac{n\overline{y}}{\lambda +1}+{\displaystyle \sum _{y=1}^{\infty }{n}_y\frac{2\left(\lambda +y+2\right)}{y\left(y+5\right)+2\lambda \left(y+2\right)+{\lambda }^2+5}},$$

(12)

in which equating it to zero and solving it implies the MLE for $\lambda$, which is denoted as $\widehat{\lambda }$. Equivalently, one can directly maximize the l directly to obtain a similar result.

3.4. Simulation Study

A simulation study is conducted to assess the performance of the obtained estimators for the parameter of the ZTPISDL distribution. The algorithm for the simulation study is similar to the one in Section 2.4, except that the data are generated using the ZTPISDL distribution. Similarly, R software is used, and the estimated parameter using MLE is obtained using the optim command to optimize the log-likelihood value. The findings of the simulation study are shown in Figure 7. Figure 7 shows that for any value of $\lambda$, as N increases, the MAD and the MSE of both estimates fall, indicating that they are both asymptotically unbiased and consistent.

When dealing with a truncated distribution, the population size is usually unknown and needs to be estimated. Assuming that a population follows the ZTPISDL distribution, a population size estimator is developed and studied. The discussion on the population size is provided in the next section.

4. Population Size Estimation

4.1. Horvitz–Thompson Estimator under ZTPISDL Distribution (HT-ZTPISDL)

A popular estimator for the population size is the Horvitz–Thompson estimator [36], which includes information from both truncated and untruncated distributions. A Horvitz–Thompson estimator has the following form:

$$\widehat{N}=\frac{n}{1-\mathrm{Pr}\left(X=0|\widehat{\omega }\right)},$$

where $\widehat{\omega }$ is the estimator for ω for a distribution and n is the sample size. Basically, the estimated parameter of the truncated distribution is substituted into the probability mass function of the untruncated distribution for the unobserved events. The estimator for the population size$\widehat{N}$, which follows a ZTPISDL distribution in the form of the Horvitz–Thompson estimator (HT-ZTPISDL), is given as

$$\widehat{N}=\frac{n{\left(\widehat{\lambda }+1\right)}^3\left({\widehat{\lambda }}^2+2\widehat{\lambda }+2\right)}{{\widehat{\lambda }}^4+6{\widehat{\lambda }}^3+13{\widehat{\lambda }}^2+8\widehat{\lambda }+2},$$

(13)

where $\widehat{\lambda }$ is the MLE for $\lambda$ in the ZTPISDL distribution.

4.2. Variance and Confidence Interval for HT-ZTPISDL

Böhning [37] has provided a simple yet informative method for obtaining a variance for any population size estimator using the conditional expectation technique. The variance of HT-ZTPISDL can be written as

$$Var\left(\widehat{N}\right)=Va{r}_{\widehat{\lambda },n}\left[\frac{n}{g\left(\widehat{\lambda }\right)}\right]=Va{r}_n\left\{E_{\widehat{\lambda }|n}\left[\frac{n}{g\left(\widehat{\lambda }\right)}\right]\right\}+E_n\left\{Va{r}_{\widehat{\lambda }|n}\left[\frac{n}{g\left(\widehat{\lambda }\right)}\right]\right\},$$

(14)

where $g\left(\lambda \right)=1-f\left(0\right)$. Observe that the variation in estimating the population size comes from two sources of variation. The first term in Equation (14) explains the binomial variation involved in sampling n units of data with population size $N$ and probability $g\left(\lambda \right)$ [37]. The second term in Equation (14) explains the variation that occurs when estimating parameter $\lambda$ using n observed data [37]. Using the delta method for the first term of Equation (14), we obtain

$$E_{\widehat{\lambda }|n}\left[\frac{n}{g\left(\widehat{\lambda }\right)}\right]=\frac{n}{g\left(\lambda \right)}.$$

As $n~binomial\left(N,g\left(\lambda \right)\right)$, we obtain

$$Va{r}_n\left\{E_{\widehat{\lambda }|n}\left[\frac{n}{g\left(\widehat{\lambda }\right)}\right]\right\}=Va{r}_n\left[\frac{n}{g\left(\lambda \right)}\right]=\frac{Ng\left(\lambda \right)\left[1-g\left(\lambda \right)\right]}{g{\left(\lambda \right)}^2}.$$

The equation above is further estimated by substituting $Ng\left(\lambda \right)$with n and$\lambda$ with $\widehat{\lambda }$, yielding

$${\widehat{Var}}_n\left\{E_{\widehat{\lambda }|n}\left[\frac{n}{g\left(\widehat{\lambda }\right)}\right]\right\}\approx \frac{n\left[1-g\left(\widehat{\lambda }\right)\right]}{g{\left(\widehat{\lambda }\right)}^2}.$$

For the ZTPISDL distribution, we obtain

$$g\left(\lambda \right)=\frac{{\lambda }^4+6{\lambda }^3+13{\lambda }^2+8\lambda +2}{{\left(\lambda +1\right)}^3\left({\lambda }^2+2\lambda +2\right)}.$$

Therefore,

$${\widehat{Var}}_n\left\{E_{\widehat{\lambda }|n}\left[\frac{n}{g\left(\widehat{\lambda }\right)}\right]\right\}\approx \frac{n{\widehat{\lambda }}^3{\left(\widehat{\lambda }+1\right)}^3\left({\widehat{\lambda }}^2+2\widehat{\lambda }+2\right)\left({\widehat{\lambda }}^2+4\widehat{\lambda }+5\right)}{{\left({\widehat{\lambda }}^4+6{\widehat{\lambda }}^3+13{\widehat{\lambda }}^2+8\widehat{\lambda }+2\right)}^2}.$$

(15)

Now, consider the second term of Equation (14), and assume that

$$E_n\left\{Va{r}_{\widehat{\lambda }|n}\left[\frac{n}{g\left(\widehat{\lambda }\right)}\right]\right\}\approx Va{r}_{\widehat{\lambda }|n}\left[\frac{n}{g\left(\widehat{\lambda }\right)}\right].$$

Using the delta method, we obtain

$$Va{r}_{\widehat{\lambda }|n}\left[\frac{n}{g\left(\widehat{\lambda }\right)}\right]=n^{2}Va{r}_{\widehat{\lambda }|n}\left[\frac{1}{g\left(\widehat{\lambda }\right)}\right]\approx n^2{\left[\frac{g^{\prime }\left(\lambda \right)}{g{\left(\lambda \right)}^2}\right]}^{2}Var\left(\widehat{\lambda }\right),$$

where $Var\left(\widehat{\lambda }\right)=I^{-1}\left(\lambda \right)/n$ and $I\left(\lambda \right)$is given in Equation (9). Therefore,

$$Va{r}_{\widehat{\lambda }|n}\left[\frac{n}{g\left(\widehat{\lambda }\right)}\right]=\frac{n{\lambda }^4{\left(\lambda +1\right)}^4{\left({\lambda }^4+10{\lambda }^3+37{\lambda }^2+52\lambda +30\right)}^2}{{\left({\lambda }^4+6{\lambda }^3+13{\lambda }^2+8\lambda +2\right)}^4}{I}^{-1}\left(\lambda \right).$$

By substituting $\lambda$ with $\widehat{\lambda }$, the second term of Equation (14) can be estimated as

$${\widehat{E}}_n\left\{Va{r}_{\widehat{\lambda }|n}\left[\frac{n}{g\left(\widehat{\lambda }\right)}\right]\right\}=\frac{n{\widehat{\lambda }}^4{\left(\widehat{\lambda }+1\right)}^4{\left({\widehat{\lambda }}^4+10{\widehat{\lambda }}^3+37{\widehat{\lambda }}^2+52\widehat{\lambda }+30\right)}^2}{{\left({\widehat{\lambda }}^4+6{\widehat{\lambda }}^3+13{\widehat{\lambda }}^2+8\widehat{\lambda }+2\right)}^4}{I}^{-1}\left(\widehat{\lambda }\right).$$

(16)

By combining Equations (15) and (16), the variance of the HT-ZTPISDL can be estimated as

$$\widehat{Var}\left(\widehat{N}\right)=\frac{n{\widehat{\lambda }}^3{\left(\widehat{\lambda }+1\right)}^3\left({\widehat{\lambda }}^2+2\widehat{\lambda }+2\right)\left({\widehat{\lambda }}^2+4\widehat{\lambda }+5\right)}{{\left({\widehat{\lambda }}^4+6{\widehat{\lambda }}^3+13{\widehat{\lambda }}^2+8\widehat{\lambda }+2\right)}^2}+\frac{n{\widehat{\lambda }}^4{\left(\widehat{\lambda }+1\right)}^4{\left({\widehat{\lambda }}^4+10{\widehat{\lambda }}^3+37{\widehat{\lambda }}^2+52\widehat{\lambda }+30\right)}^2}{{\left({\widehat{\lambda }}^4+6{\widehat{\lambda }}^3+13{\widehat{\lambda }}^2+8\widehat{\lambda }+2\right)}^4}{I}^{-1}\left(\widehat{\lambda }\right).$$

(17)

Therefore, the 95% confidence interval for the estimator can be written as $\widehat{N}\mp z_{0.025}\sqrt{\widehat{Var}\left(\widehat{N}\right)}$, where $z_{0.025}=1.96$.

4.3. Simulation Study

A simulation study is conducted to assess the performance of the HT-ZTPISDL estimator in estimating the population size when the data are generated from the ZTPISDL distribution. The algorithm for the simulation study is as follows:

Step 1:

Generate $N=1000, 2000, \dots , 10,000$random data, which follow the ZTPISDL distribution with $\lambda =0.5$.

Step 2:

Obtain $\widehat{\lambda }$ using the MLE and use $\widehat{\lambda }$ to obtain $\widehat{N}$.

Step 3:

Repeat Steps 1–2 for a total of 2000 iterations and obtain the estimates.

Step 4:

Calculate the relative absolute error, RAB values, and the relative standard deviation, RSd values, given, respectively, as $RAB=\left|\overline{\widehat{N}}-N\right|/N$ and $RSd=\frac{1}{N}\sqrt{{\sum }_{i=1}^{2000}{\left(\widehat{N}-\overline{\widehat{N}}\right)}^2/2000}$, where $\overline{\widehat{N}}={\sum }_{i=1}^{2000}\widehat{N}/2000$.

Step 5:

Repeat Steps 3–4 for $\lambda =2.0, 5.0$.

R software is used to obtain $\widehat{\lambda }$ using the optim command, and $\widehat{N}$ is obtained by plugging in $\widehat{\lambda }$ in Equation (10). The results of the simulation study are presented in Figure 8. In this figure, for any value of $\lambda$, as N increases, the MAD and the MSE of $\widehat{N}$ decreases, suggesting that $\widehat{N}$ is asymptotically unbiased and consistent.

5. Medical Data Applications

The applications in the medical datasets are segregated into three subsections with respect to the PISDL distribution, the ZTPISDL distribution, and the estimation of the population size. For a comparison of the model fittings, Akaike’s information criterion, AIC [38], and Bayesian information criterion, BIC [39], are used.

5.1. Model Fittings Using the PISDL Distribution

Two datasets on the number of dicentric chromosomes after being exposed to different doses of radiation (0.405 and 0.600) that were studied by Puig and Barquinero [40] are considered in fitting using the PISDL distribution. The two datasets are overdispersed with dispersion values of 1.2704 and 1.2178, respectively. Since the PISDL distribution is closely related to the Poisson and the PL distributions, the model fittings from the PISDL distribution are compared with those from the Poisson and the PL distributions. The results of the model fittings for the two datasets using the three distributions are summarized in

Table 1. Model fittings of the number of dicentric chromosomes after being exposed to a 0.405 radiation dose using Poisson, PL, and PISDL distributions.

$\mathit{x}$	${\mathit{n}}_{\mathit{x}}$	Distributions
		Poisson	$\mathit{P}\mathit{L}$	$\mathit{P}\mathit{I}\mathit{S}\mathit{D}\mathit{L}$ (MLE)	$\mathit{P}\mathit{I}\mathit{S}\mathit{D}\mathit{L}$ (Moment)
0	437	426.55	433.87	433.73	433.64
1	66	84.50	72.04	72.30	72.36
2	15	8.37	11.81	11.75	11.77
3	1	0.55	1.91	1.87	1.88
4	1	0.03	0.37	0.35	0.37
Total	520	520.00	520.00	520.00	520.00
Parameter
$\widehat{\lambda }$		0.1981	-	6.5464	6.5392
$\widehat{\theta }$		-	5.7953	-	-
Max log-likelihood		−285.14	−279.40	−279.45
$AIC$		572.27	560.80	560.90	-
$BIC$		576.52	565.05	565.15	-
${\chi }^2$		11.55	1.13	1.23	1.22
df		1	1	1	1
p-value		0.0007	0.2878	0.2674	0.2694

and

Table 2. Model fittings of the number of dicentric chromosomes after being exposed to a 0.600 radiation dose using Poisson, PL, and PISDL distributions.

$\mathit{x}$	${\mathit{n}}_{\mathit{x}}$	Distributions
		Poisson	$\mathit{P}\mathit{L}$	$\mathit{P}\mathit{I}\mathit{S}\mathit{D}\mathit{L}$ (MLE)	$\mathit{P}\mathit{I}\mathit{S}\mathit{D}\mathit{L}$ (Moment)
0	473	456.69	475.79	474.91	475.07
1	119	147.65	117.81	118.96	118.88
2	34	23.87	28.53	28.55	28.50
3	3	2.57	6.79	6.64	6.62
4	2	0.22	2.08	1.94	1.93
Total	631	631.00	631.00	631.00	631.00
Parameter
$\widehat{\lambda }$		0.3233	-	4.4001	4.4046
$\widehat{\theta }$		-	3.7420	-	-
Max log-likelihood		−469.65	−464.09	−463.96
$AIC$		941.30	930.17	929.91	-
$BIC$		945.75	934.61	934.36	-
${\chi }^2$		11.85	2.77	2.54	2.54
df		1	2	2	2
p-value		0.0006	0.2503	0.2808	0.2808

.

In both tables, the Poisson distribution does not fit the data based on the p-value. On the other hand, the model fittings based on both PL and PISDL distributions gave similar AIC and BIC values, as well as non-significant p-values, indicating that both distributions can be used for describing the number of dicentric chromosomes after being exposed to radiation of different doses. However, the first dataset was fitted better by the PL distribution, whereas the second dataset was fitted better by the PISDL distribution. Therefore, it is reasonable to suggest that both PL and PISDL distributions are equally competent and can be selected as the best distributions in explaining the number of dicentric chromosomes after being exposed to two different doses of radiation.

The comparison between the empirical plots of the data and the fitted values based on Equation (1) in Theorem 1 for the considered datasets, (i) the number of dicentric chromosomes after being exposed to a 0.405 radiation dose, and (ii) the number of dicentric chromosomes after being exposed to a 0.600 radiation dose are presented in Figure 9.

5.2. Model Fittings Using ZTPISDL Distribution

A dataset on the number of positive samples of Salmonella data, which was initially given in a survey study by Snow et al. [41] and later summarized by Arnold et al. [42], is considered in model fitting using the zero-truncated Poisson (ZTP), the zero-truncated PL (ZTPL), and the ZTPISDL distributions. The data refer to the number of farms with at least one positive sample of Salmonella. The dataset is overdispersed with a dispersion value of 1.4381. The results of the model fittings are given in

Table 3. Model fittings of the number of positive samples of Salmonella using the ZTP, ZTPL, and ZTPISDL distributions.

$\mathit{y}$	${\mathit{n}}_{\mathit{y}}$	Distributions
		$\mathit{Z}\mathit{T}\mathit{P}$	$\mathit{Z}\mathit{T}\mathit{P}\mathit{L}$	$\mathit{Z}\mathit{T}\mathit{P}\mathit{I}\mathit{S}\mathit{D}\mathit{L}$ (MLE)	$\mathit{Z}\mathit{T}\mathit{P}\mathit{I}\mathit{S}\mathit{D}\mathit{L}$ (Moment)
1	17	7.88	15.06	14.01	14.00
2	9	12.12	11.55	11.62	11.62
3	5	12.44	8.45	8.82	8.82
4	6	9.57	5.99	6.32	6.32
5	5	5.89	4.15	4.34	4.34
6	5	3.02	2.83	2.89	2.89
7	6	2.08	4.97	5.00	5.01
Total	53	53.00	53.00	53.00	53.00
Parameter
$\widehat{\lambda }$		3.0778	-	0.8932	0.8928
$\widehat{\theta }$		-	0.6660	-	-
Max log-likelihood		−110.64	−105.28	−105.10
$AIC$		223.27	212.55	212.19	-
$BIC$		225.24	214.53	214.16	-
${\chi }^2$		24.10	3.61	4.16	4.16
df		4	3	4	4
p-value		<0.0001	0.3068	0.3848	0.3848

.

In Table 3, the ZTP distribution does not provide a good fit to the data based on the p-value. On the contrary, the ZTPISDL and the ZTPL distributions provide a good fit to the data based on the AIC and BIC values, as well as the non-significant p-value. However, the dataset was fitted better by the ZTPISDL distribution. Therefore, this suggests that the ZTPISDL is the best distribution in describing the number of positive samples of Salmonella data. The comparison between the empirical plots of the data and the fitted value based on Equation (10) in Theorem 2 for the above dataset is presented in Figure 10.

5.3. Estimating Population Size

The dataset studied in Section 5.2 on the number of positive samples of Salmonella data is considered to estimate the population size of the sample. The Horvitz–Thompson estimator based on the proposed ZTPISDL distribution is compared with those based on the ZTP and the ZTPL distributions. The Horvitz–Thompson estimators based on the ZTP [43,44] and the ZTPL [32] distributions are, respectively, given as

$${\widehat{N}}_{ZTP}=\frac{n}{1-\exp \left(-{\widehat{\lambda }}_{ZTP}\right)},$$

and

$${\widehat{N}}_{ZTPL}=\frac{n{\left(\widehat{\theta }+1\right)}^3}{{\widehat{\theta }}^2+3\widehat{\theta }+1},$$

where ${\widehat{\lambda }}_{ZTP}$refers to the MLE of $\lambda$ for the ZTP distribution and $\widehat{\theta }$refers to the MLE of $\theta$ for the ZTPL distribution. The estimated population sizes and their corresponding standard deviations, as well as the lower and upper limits for the 95% confidence interval, are presented in

Table 4. The estimated population size, the standard deviation, and the lower and upper limits for the 95% confidence interval of the population size estimator based on the ZTP, ZTPL, and ZTPISDL distributions.

Distributions	$\mathbf{Estimated} \widehat{\mathit{N}}$	$\mathit{S}\mathit{D}$	95% Lower Limit	95% Upper Limit
$ZTP$	55.56	1.761	52.11	59.01
$ZTPL$	71.21	4.045	63.28	79.14
$ZTPISDL$ (MLE)	66.64	4.896	57.04	76.24
$ZTPISDL$ (Moment)	63.10	4.961	53.38	72.82

. Based on Table 4, the estimated population size based on the ZTPISDL distribution using the MLE is 66.64 with a 95% confidence interval between 57.04 and 76.24. Since the ZTPISDL distribution based on the MLE provides the best fit for the data (refer to Table 3), the resulting estimated population size is acceptable.

6. Conclusions, Limitations, and Future Research

A tractable one-parameter Poisson-improved second-degree Lindley (PISDL) distribution has been proposed to address the need for modeling count data exhibiting overdispersion. This distribution is composed of three negative binomial distributions, each with fixed mixing proportions and parameters, allowing for the fitting of datasets originating from three sub-populations whose modes are proximal. However, if the three modes are far from each other and clearly visible from the plots of the datasets, the PISDL distribution may not be a good candidate for model fittings. The hazard rate function for the PISDL distribution showed an increasing shape. Parameters of the PISDL distribution have been estimated using maximum likelihood estimation (MLE) and moment methods, and both were found to be asymptotically unbiased and consistent. It has been observed from model fittings that the PISDL distribution performs on par with the PL distribution and surpasses the standard Poisson distribution in describing the number of dicentric chromosomes post-exposure to various radiation doses.

Given that data may not always present the frequency of unobserved events and exhibit dispersion, methodologies like zero truncation or the size-biased approach are employed. Zero truncation is a commonly favored method for handling datasets lacking frequencies of non-observed events. Hence, a zero-truncated version of the proposed distribution, named the zero-truncated PISDL (ZTPISDL) distribution, has been proposed to accommodate data exhibiting either over- or underdispersion. Parameters of the ZTPISDL distribution estimated by the MLE and moment methods have also been shown to be asymptotically unbiased and consistent. When applied to datasets, the MLE technique for estimating parameters of the ZTPISDL distribution has provided the best fit in comparison with zero-truncated Poisson and zero-truncated PL distributions.

When the population size is unknown due to the absence of frequencies for non-observed events in positive count data, estimation has been conducted using the Horvitz–Thompson estimator based on the ZTPISDL distribution. Since the ZTPISDL distribution has provided the best fit for the dataset considered in this study, the acceptance of the resulting estimated population size for the number of positive Salmonella samples is justified. It is suggested that this population size estimate may serve as a lower bound to the actual population size, especially when the ZTPISDL distribution is extended to linear models with the inclusion of relevant covariates. Moreover, the derivation of the variance and confidence interval for the population size estimator is intended to assist policymakers in revising rules and guidelines pertaining to the population, which will, in turn, benefit the population at large.

Despite its flexibility, the PISDL distribution’s mixing proportion is solely dependent on λ, which consequently limits its flexibility. By introducing a new parameter (denoted as α) to influence the mixing proportion alongside λ, greater flexibility is achieved, resulting in a more adaptable PISDL distribution. This would, in turn, lead to a more versatile ZTPISDL distribution and improved estimates for population size.

Further research is anticipated to explore additional modifications and applications of the PISDL distribution. These include, but are not limited to, actuarial measures, such as value-at-risk and tail value-at-risk, reliability measures, like hazard rate and entropy, various forms of inflated models, including zero-inflated, k-inflated, and zero-one-inflated distributions, and weighted models, like size-biased and area-biased distributions. These enhancements are expected to broaden the utility of the PISDL distribution, making it a competitive model in the realm of the statistical literature.