A Family of Finite Mixture Distributions for Modelling Dispersion in Count Data

Abstract

This paper considers the construction of a family of discrete distributions with the flexibility to cater for under-, equi- and over-dispersion in count data using a finite mixture model based on standard distributions. We are motivated to introduce this family because its simple finite mixture structure adds flexibility and facilitates application and use in analysis. The family of distributions is exemplified using a mixture of negative binomial and shifted negative binomial distributions. Some basic and probabilistic properties are derived. We perform hypothesis testing for equi-dispersion and simulation studies of their power and consider parameter estimation via maximum likelihood and probability-generating-function-based methods. The utility of the distributions is illustrated via their application to real biological data sets exhibiting under-, equi- and over-dispersion. It is shown that the distribution fits better than the well-known generalized Poisson and COM–Poisson distributions for handling under-, equi- and over-dispersion in count data.

1. Introduction

In the statistical analysis of count data, an important feature is the under-, equi- or over- dispersion of the data. There are many models that cater for a particular type of dispersion. In contrast, there are not many distributions with the flexibility to model under-, equi- and over-dispersion. Examples of these distributions are the hyper-Poisson [1], double Poisson [2] and weighted Poisson [3] distributions. Among these distributions, two popular distributions are the generalized Poisson distribution [4] and the COM–Poisson distribution [5,6,7,8].

The generalized Poisson distribution (GPD) has a probability mass function (pmf)

$$P(X=x)=\left\{\begin{array}{cc}\hfill \frac{\theta (\theta +x\lambda {)}^{x-1}{e}^{-\theta -x\lambda }}{x!}\hfill & \mathrm{for} x=0, 1, 2, \dots \hfill \\ \hfill 0\hfill & \mathrm{for} x>m, \mathrm{when} \lambda <0\hfill \end{array}\right.$$

where $\theta >0$, max $(-1,-\theta /m)<\lambda \le 1$ and $m (\ge 4)$ are the largest positive integers satisfying $\theta +m\lambda >0$ when $\lambda$ is negative. For $\lambda <0$, the probabilities do not sum to 1. The GPD is under-, equi- or over-dispersed when λ is negative, zero or positive, respectively. If $\lambda =0$, the Poisson distribution is obtained. However, for under-dispersion, the parameters must be restricted for the GPD to be a proper distribution [9]; see also [10] for further discussion.

The COM–Poisson distribution arises from a single queue–single server system with state-dependent service rates [5]. Recently, this distribution has attracted a lot of interest and has been applied in many areas [6]. Refs. [11,12] considered the COM–Poisson distribution in the context of modelling time series of counts. The COM–Poisson distribution is a member of the exponential family of distributions and confers several desirable properties for statistical inference. The COM–Poisson distribution has a pmf given by

$$P(X=x)=\frac{{\lambda }^x}{(x!{)}^{v}Z(\lambda ,v)}$$

where $Z(\lambda ,v)={\sum }_{j=0}^{\infty }\frac{{\lambda }^j}{(j!{)}^v}$for $\lambda >0$ and $v\ge 0$. The COM–Poisson distribution is under-dispersed or over-dispersed when $\nu >1$ or $\nu <1,$ respectively. The case $\nu =1$ gives the Poisson distribution. The COM–Poisson distribution cannot be directly expressed in terms of its mean. To facilitate regression modelling, ref. [13] considered a mean-parametrized COM–Poisson model based on an iterative solution of a mean equation. Compared with over-dispersion, there is relatively less discussion in the literature about under-dispersion in count data. Ref. [14] reviewed under-dispersed count models and various sources of data under-dispersion. Ref. [15] discussed count models that are arbitrarily under-dispersed and showed that the mean-parametrized COM–Poisson distribution can handle arbitrarily small under-dispersion.

It is to be noted that the distributions mentioned above have a restricted parameter range (GPD) or normalizing constants in the form of infinite series (double Poisson, COM–Poisson) that may add to the complexity of a statistical analysis. Recently, motivated by these issues, a particular case of the generalized shifted inverse trinomial (GIT) distribution [16], designated $\mathrm{GI}{\mathrm{T}}_{\mathrm{3,1}}$ [17], has been proposed in order to cater for under-, equi- or over-dispersion in count data. It is noted that when the GIT is equi-dispersed, the distribution is non-Poisson. The GIT has a simple probabilistic structure as the convolution of negative binomial and binomial distributions. This is apparent from its probability-generating function (pgf), given by

$$G(t)={\left(\frac{1-p_3}{1-p_{3}t}\right)}^n{\left(\frac{p_1+p_{2}t}{p_1+p_2}\right)}^n,1-p_3=p_1+p_2,n>0.$$

It is a member of a general family of distributions formed via the convolution of binomial and pseudo-binomial variables [18]. This family has a pgf in the form

$${\left(\frac{1-Q_{1}t}{1-Q_1}\right)}^{U_1}{\left(\frac{1-Q_{2}t}{1-Q_2}\right)}^{U_2},$$

and Kemp gave conditions for this to be a valid pgf. Some members of this family of distributions have been considered in the literature.

The pmf of the GIT distribution is given by

$$f_n(x)={\sum }_{i=0}^{\mathit{min}(n,x)}\frac{n}{n+x-i}\left(\begin{array}{c}n+x-i\\ n-i,i,x-i\end{array}\right)p_1^{n-i}{p}_2^i{p}_3^{x-i}$$

for $x=\mathrm{0,1},2,\dots$, where n is a positive integer and $p_i\ge 0 \mathrm{for} i=\mathrm{1,2},3$; ${\sum }_{i=1}^3{p}_i=1$.

The GIT pmf may also be expressed as

$$f_n(x)=\left(\begin{array}{c}n+x-1\\ x\end{array}\right)p_1^n{p}_3^x{}_{2}F_1\left(-n,-x;-n-x+1;-\frac{p_2}{p_1{p}_3}\right)$$

in terms of the Gauss hypergeometric function ${}_{2}F_1$. Other recent efforts to develop distributions that can handle under-, equi- or over-dispersion have been performed through COM–Poisson-type extensions on binomial and negative binomial distributions [19,20,21,22,23,24]. Ref [25] examined the flexible weighted Poisson that extends the weighted Poisson of [3] and contains the COM–Poisson and hyper-Poisson distributions as special cases.

Even though the GIT has a simple probabilistic structure as a convolution of binomial and negative binomial random variables and is free from some of the limitations of the other popular distributions, the positive integer index parameter n can be an impediment in parameter estimation. Section 2 proposes a finite mixture model for the construction of a family of distributions that can cater for under-, equi- or over-dispersion. This is exemplified by a finite mixture of negative binomial and shifted negative binomial distributions. The negative binomial distribution is over-dispersed, but the finite mixture exhibits under-, equi- and over-dispersion. Section 2 also presents some basic properties of the finite mixture of negative binomial distributions, while Section 3 discusses parameter estimation and test of equi-dispersion. Empirical modelling of biological count data is considered in Section 4, with concluding remarks presented in Section 5.

2. Finite Mixture of Distributions

Finite mixture models have been widely employed due to their flexibility and computationally convenient representation of complex distributions. Mixture formulations lead to simple techniques in statistical analysis; for example, ref [26] gave a simple method for generating bivariate binomial samples with given marginals and correlation. A good review of this topic is found in [27]. A general family of two-component finite mixture distributions is first established and then a particular case of interest is examined.

2.1. Finite Mixture of Distribution and Its Shifted Distribution

Let $P\left(k;\mathsf{\Theta }\right)$ be a probability mass function (pmf), where $\mathsf{\Theta }$ is the vector of parameters, $k=0, 1, 2, 3,\dots .$ and $P\left(k-1\right)$ is the shifted pmf. Define the pmf $P\left(k;\mathsf{\Theta },p\right)$ of a random variable $X$ as

$$P\left(k;\mathsf{\Theta },p\right)=\left(1-p\right)P\left(k;\mathsf{\Theta }\right)+pP\left(k-1;\mathsf{\Theta }\right),k=0,1,2,3,\dots ,$$

(1)

where $P\left(-1;\mathsf{\Theta }\right)=0$ and $0<p<1$. It is easy to check that

$$\sum _{k=0}^{\mathrm{\infty }}P\left(k;\mathsf{\Theta },p\right)=1.$$

If $G\left(z\right)$ is the probability-generating function (pgf) of $P\left(k;\mathsf{\Theta }\right)$, the pgf of $X$ is

$$G_X\left(z\right)=\left(1-p\right)G\left(z\right)+pzG\left(z\right).$$

Moments of the finite mixture may be obtained from this pgf. Let $\mu$ and ${\sigma }^2$ be the mean and variance of the component distribution $P\left(k;\mathsf{\Theta }\right)$. The mean and variance of $X$ with pmf given by Equation (1) are

$${\mu }_X=\mu +p,{\sigma }_X^2={\sigma }^2+p\left(1-p\right).$$

(2)

The index of dispersion $ID$ is given by

$$ID=\frac{{\sigma }_X^2}{{\mu }_X}=\frac{{\sigma }^2+p\left(1-p\right)}{\mu +p}.$$

The distribution of finite mixture (1) is under-, equi- or over-dispersed according to

$$ID<,=,>1,$$

that is,

$${\sigma }^2-\mu <,=,>p^2.$$

The quantity ${\sigma }^2-\mu <,=,>0$ depends upon whether the component distribution is under-, equi- or over-dispersed. This shows that the choice of either an equi-dispersed or over-dispersed component distribution will make the finite mixture (1) an under, equi- or over-dispersed distribution. This is illustrated in the next section by using a mixture of negative binomial (NB) distributions.

2.2. Finite Mixture of Negative Binomial Distributions

The negative binomial distribution is a popular model that only caters for over-dispersion. However, when used as components in the finite mixture model (1), the resulting mixture model can handle under-, equi- or over-dispersion. The mixture of NB distributions has pmf

$$\begin{gathered} P\left(0\right)=\left(1-p\right){\left(1-\theta \right)}^v,k=0 \\ P(k)=\left(1-p\right)\frac{{\left(v\right)}_k}{k!}{\theta }^k{\left(1-\theta \right)}^v+p\frac{{\left(v\right)}_{k-1}}{\left(k-1\right)!}{\theta }^{k-1}{\left(1-\theta \right)}^v,k\ge 1. \end{gathered}$$

(3)

where $0<p<\mathrm{1,0}<\theta <1,\nu >0 \mathrm{a}\mathrm{n}\mathrm{d}{\left(v\right)}_k=\frac{\mathsf{\Gamma }\left(v+k\right)}{\mathsf{\Gamma }\left(v\right)}$. Let $NB\left(\theta ,\nu \right)$ denote a NB random variable with parameters $\left(\theta ,\nu \right)$. Thus (3) may be interpreted as follows: an $NB\left(\theta ,\nu \right)$ chance mechanism contributes a proportion $1-p$ of the counts $k\ge 0$. However, for non-zero counts where $k\ge 1$, a shifted $NB\left(\theta ,\nu \right)$ also contributes a proportion $p$.

The pmf (3) can be written as

$$\begin{gathered} P(k)=\left(1-p\right)\frac{{\left(v\right)}_k}{k!}{\theta }^k{\left(1-\theta \right)}^v+p\frac{k}{\theta \left(\nu +k-1\right)}\frac{{\left(v\right)}_k}{k!}{\theta }^k{\left(1-\theta \right)}^v \\ =\left[\left(1-p\right)+p\frac{k}{\theta \left(\nu +k-1\right)}\right]\frac{{\left(v\right)}_k}{k!}{\theta }^k{\left(1-\theta \right)}^v,k\ge 0 \end{gathered}$$

(4)

It is easy to compute the probability $P\left(k\right)$ from (4) by using the one-term recurrence formula for the NB pmf

$$P\left(k\right)=\frac{{\left(v\right)}_k}{k!}{\theta }^k{\left(1-\theta \right)}^v,$$

that is, $P(k+1)=\left(v+k\right)\theta P\left(k\right)/\left(k+1\right)$.

The pgf of pmf (3) is given by

$$G(z)=(1-p){\left(\frac{1-\theta }{1-\theta z}\right)}^v+pz{\left(\frac{1-\theta }{1-\theta z}\right)}^v$$

(5)

The form of the pgf given by (5) shows that this finite mixture of NB distributions with pmf (3) can be called an NB-shifted NB distribution. The pgf (5) may also be written as

$$G\left(z\right)=\left[\left(1-p\right)+pz\right]{\left(\frac{1-\theta }{1-\theta z}\right)}^v,\left|z\right|<\frac{1}{\left(1-\theta \right)}.$$

This is seen to be the convolution of Bernoulli and NB distributions. The mean and variance of this special case are, respectively,

$$\mu =p+\frac{\nu \theta }{1-\theta }$$

and

$${\sigma }^2=p\left(1-p\right)+\frac{\nu \theta }{{\left(1-\theta \right)}^2}.$$

2.3. Weighted Negative Binomial Distribution

Let $X$ be a random variable with pmf p(k) and assume that the probability of ascertaining the event $X=k$ has a weighting factor w(k). The weighted distribution [28,29] with weight w(x) has a pmf

$$P\left(k\right)=w\left(k\right)p\left(k\right)/\sum w\left(x\right)p\left(x\right).$$

(6)

If w(k) = k, the distribution with pmf (6) is known as the size(length)-biased distribution.

Let $Y$ be the size-biased version of a random variable $X$ with weight $w\left(k\right)=k.$ Consider another size-biased random variable $X^w$ with weight $w\left(k\right)=k+\gamma$. The pmf of $X^w$ is given by

$$P(X^w=k)=\frac{\gamma }{\mu +\gamma }P(X=k)+\frac{\mu }{\mu +\gamma }P(Y=k)$$

(7)

where $\mu =\sum xp\left(x\right)$. If $\gamma \to \mathrm{\infty }$, (7) reduces to $P(X=k).$

By comparing Equations (4) and (7), it is seen that the finite mixture of NB can be regarded as a weighted NB distribution with weight

$$w\left(k\right)=\left(1-p\right)+p\frac{k}{\theta \left(\nu +k-1\right)}.$$

2.4. Conditions for Under-, Equi- and Over-Dispersion

The index of dispersion is

$$ID=\frac{\mathrm{Variance}}{\mathrm{Mean}}=\frac{p\left(1-p\right)+\frac{\nu \theta }{{\left(1-\theta \right)}^2}}{p+\frac{\nu \theta }{1-\theta }}.$$

The distribution is under, equi- and over-dispersed according to whether $ID$ is less than, equal to or more than 1. Thus, the conditions for under, equi- and over-dispersion are as follows:

If $\sqrt{\nu }\theta /\left(1-\theta \right)<p$, the distribution is under-dispersed. Conversely, it is equi-dispersed if $\sqrt{\nu }\theta /\left(1-\theta \right)=p$ and over-dispersed if $\sqrt{\nu }\theta /\left(1-\theta \right)>p$.

Remark 1.

It should be noted that for the special case of equi-dispersion, the distribution does not reduce to the Poisson distribution. With the substitution$p=\sqrt{\nu }\theta /\left(1-\theta \right)$into Equation (4), a weighted negative binomial distribution is obtained.

2.5. Shapes of the Distribution

To examine the shapes of NB-shifted NB distribution, the probabilities in

Table 1. Different combination of parameters under NB-shifted NB distribution.

	$\mathit{\theta }$$\mathbf{and} \mathit{p} \mathbf{Are} \mathbf{Fixed} \mathbf{and} \mathit{\nu }$ Is Increased			$\mathit{\nu }$$\mathbf{Is} \mathbf{Fixed} \mathbf{and} \mathit{\theta }$ and p Are Varied
	Case (a)	Case (b)	Case (c)	Case (d)	Case (e)	Case (f)
$\nu$	1	10	50	10	10	10
$\theta$	0.1	0.1	0.1	0.1	0.4	0.8
p	0.1	0.1	0.1	0.8	0.4	0.1
ID	1.01	1.09	1.11	0.73	1.61	4.99

are computed for a different combination of parameters.

Figure 1i shows that increasing the parameter $\nu$ decreases the probability at zero count dramatically but increases the index of dispersion mildly. The computed pmf drops drastically from 0.81 (case (a)) to almost zero (case (c)) when $\nu$ increases from 1 to 50. For case (a), the computed index of dispersion is close to 1 and can be considered to be an equi-dispersed distribution. When $\theta$ and p are fixed and $\nu$ (>1) is increased as shown in Table 1, the computed index of dispersion is always larger than 1, meaning that it is an over-dispersed distribution. Besides that, the modes are shifted to the right when $\nu$ is increased. An increase to the length of the distribution is found and it is always skewed where the right tails are heavier.

On the other hand, Table 1 shows that, when $\nu$ is fixed and $\theta$ and p are varied, we obtain either over- or under- dispersion. The length of the distribution is longer when $\theta$ is larger than p, as presented in Figure 1ii. Meanwhile, case (f) has the highest computed index of dispersion. However, the distribution seems to be flattened when compared with other cases.

2.6. Log-Concavity, Strong Unimodality and Reliability Properties

A distribution is said to be log-concave if its pmf $\left\{f_k\right\}$, $f_k>0,k=0,\mathrm{1,2},3,\dots$ satisfies $f_k^2\ge f_{k+1}{f}_{k-1}$.

From the results of [30] (p. 388), we can assert that the NB-shifted NB distribution is log-concave since both the binomial (Bernoulli) and negative binomial distributions are log-concave and as log-concavity is preserved under convolution.

From Theorem 3 of [30] (p. 386), which states that a necessary and sufficient condition of pmf $\left\{f_k\right\}$ being strongly unimodal is that $f_k$ be log-concave for all k, it follows that the finite mixture of NB distributions is strongly unimodal.

The failure rate of a distribution with pmf $\left\{f_k\right\}$, $f_k>0,k=0, \mathrm{1,2},3,\dots$ is defined by

$$r\left(k\right)=f_k/{\displaystyle {\sum }_{i\ge k}^{\infty }}{f}_i.$$

From the log-concavity property, the NB-shifted NB distribution has an increasing failure rate (IFR) [31]. Furthermore, the following implications hold

$$IFR\Rightarrow IFRA\Rightarrow NBU\Rightarrow NBUE\Rightarrow HNBUE$$

where IFRA is ‘increasing failure rate average’, NBU is ‘new better than used’, NBUE is ‘new better than used in expectation’ and HNBUE is ‘harmonic new better than used in expectation’. Thus, the finite mixture of NB-shifted NB distribution is IFR, IFRA, NBU, NBUE and HNBUE.

3. Statistical Inferences

3.1. Parameter Estimation

In this section, maximum likelihood (ML) estimation and probability-generating function-based estimation (pgf-estimator) are employed.

ML estimation is an efficient method for the estimation of unknown parameters, but the score equations involved can be complicated and difficult to solve. To overcome numerical complexity, numerical optimization via the simulated annealing algorithm [32] is used.

An alternative estimation method based on the pgf is suggested because the NB-shifted NB pgf has a simple form. This method has been shown to provide quick and consistent estimates for discrete distributions and is robust to outliers [33]. The pgf-based statistic considered here is

$$T={\int }_0^1\left(F_N\right(t)-F(t){)}^{2}dt$$

(8)

where $F_N(t)=\frac{1}{N}{\sum }_{i=1}^N{t}^{x_i}$ and $F(t)=E\left[t^X\right]$ are, respectively, the empirical pgf and the theoretical pgf of the distribution. Ref [33] demonstrated that the pgf estimator outperformed the ML estimator in terms of achieving low values of mean-square error and bias.

Remark 2.

In the statistic (8), the integral over (0, 1) may be interpreted as averaging over the auxiliary variable t by a uniform distribution. A non-uniform distribution can be used.

3.2. Test for Equi-Dispersion

The NB-shifted NB distribution is equi-dispersed when $p=\sqrt{\nu }\theta /\left(1-\theta \right)$. Hence, the following set of hypotheses is considered. Let $h\left(p, \nu , \theta \right)=p-\sqrt{\nu }\theta /\left(1-\theta \right)$

$$H_0:h\left(p,\nu ,\theta \right)=0 \mathrm{vs}.H_1:h\left(p,\nu ,\theta \right)\ne 0.$$

The study of the power of this statistical hypothesis test is developed in the next subsection and the results will be provided in

Table 2. Simulated power of Score and GLR tests for NB-shifted NB distribution. (Power = number of rejections divided by number of repetitions).

			Equi-Dispersion	Over-Dispersion		Under-Dispersion
$\nu$			10	10	40	25	3
$\theta$			0.1	0.1	0.1	0.1	0.1
$p$			0.35	0.1	0.1	0.8	0.8
Index of dispersion			1	1.09	1.11	0.91	0.47
Effect size			0	0.25	0.60	0.24	0.61
N	$\alpha$	Method
100	0.05	score	0.030	0.086	0.114	0.010	0.048
		GLR	0.036	0.077	0.105	0.072	0.911
	0.1	score	0.056	0.148	0.160	0.030	0.362
		GLR	0.082	0.133	0.175	0.144	0.960
500	0.05	score	0.035	0.280	0.391	0.059	1.000
		GLR	0.050	0.274	0.389	0.279	1.000
	0.1	score	0.069	0.394	0.493	0.124	1.000
		GLR	0.099	0.387	0.509	0.407	1.000
1000	0.05	score	0.032	0.494	0.637	0.152	1.000
		GLR	0.046	0.483	0.654	0.538	1.000
	0.1	score	0.080	0.626	0.742	0.286	1.000
		GLR	0.088	0.621	0.758	0.655	1.000

.

3.2.1. Rao’s Score Test [34]

The log-likelihood function can be written as $\ln L={\sum }_i^N{I}_{\left\{Y_i=y_i\right\}}lnP(Y_i=y_i),$ where $Y_1, Y_2, {\cdots , Y}_N$ is a random sample. The Rao’s score test statistic is given by T = $V{I}^{-1}{V}^T$, where the score vector V and information matrix I are evaluated at the restricted ML estimates.

The score vector $V$ is given as $\left(\frac{\partial \ln L}{\partial p},\frac{\partial \ln L}{\partial \nu },\frac{\partial \ln L}{\partial \theta }\right)$, where

$$\begin{gathered} \frac{\partial \ln L}{\partial p}={\sum }_i^N{I}_{\left\{y_i\ge 0\right\}}\frac{\partial lnP\left(Y_i=y_i\right)}{\partial p} \\ ={\sum }_i^N{I}_{\left\{y_i\ge 0\right\}}\left\{\frac{y_i\left(\theta -1\right)+(v-1)\theta }{\alpha \left(i\right)}\right\}, \end{gathered}$$

where $\alpha \left(i\right)=y_i\left(p\left(\theta -1\right)-\theta \right)+\left(p-1\right)\left(v-1\right)\theta .$

$$\begin{gathered} \frac{\partial \ln L}{\partial \nu }={\sum }_i^N{I}_{\left\{y_i\ge 0\right\}}\frac{\partial lnP(Y_i=y_i)}{\partial \nu } \\ ={\sum }_i^N{I}_{\left\{y_i\ge 0\right\}}\left\{\frac{y_{i}p}{\left(y_i+v-1\right)\alpha \left(i\right)}-{\psi }_0\left(v\right)+{\psi }_0\left(v+y_i\right)+\ln \left(1-\theta \right)\right\} \end{gathered}$$

$$\begin{gathered} \frac{\partial \mathit{ln}L}{\partial \theta }={\sum }_i^N{I}_{\left\{y_i\ge 0\right\}}\frac{\partial lnP(Y_i=y_i)}{\partial \theta } \\ ={\sum }_i^N{I}_{\left\{y_i\ge 0\right\}}\left\{\left(\frac{\nu }{\theta -1}\right)+\frac{y_i}{\theta }\left[1+\frac{p}{\alpha \left(i\right)}\right]\right\}. \end{gathered}$$

The elements of the information matrix are

$$\begin{array}{l}\mathrm{E}\left[\frac{{\partial }^2\mathit{ln}L}{\partial p^2}\right]=-{\sum }_i^{N}E\left(I_{\left\{y_i\ge 0\right\}}\right){\left(\frac{y_i\left(\theta -1\right)+(v-1)\theta }{\alpha \left(i\right)}\right)}^2,\\ \mathrm{E}\left[\frac{{\partial }^2\mathit{ln}L}{\partial p\partial \nu }\right]=-{\sum }_i^{N}E\left(I_{\left\{y_i\ge 0\right\}}\right)\frac{y_i\theta }{{\alpha \left(i\right)}^2},\\ \mathrm{E}\left[\frac{{\partial }^2\mathit{ln}L}{\partial p\partial \theta }\right]=-{\sum }_i^{N}E\left(I_{\left\{y_i\ge 0\right\}}\right)\frac{y_i(y_i+v-1)}{{\alpha \left(i\right)}^2},\\ \mathrm{E}\left[\frac{{\partial }^2\mathit{ln}L}{\partial \nu \partial \theta }\right]={\sum }_i^{N}E\left(I_{\left\{y_i\ge 0\right\}}\right)\left\{\frac{1}{\theta -1}-\frac{{y_i}^2{p}^2}{{\left(y_i+v-1\right)}^3{\left[1+p\left(\frac{y_i}{y_i+v-1}-1\right)\right]}^2{\theta }^3}+\frac{y_{i}p}{{\left(y_i+v-1\right)}^2\left[1+p\left(\frac{y_i}{y_i+v-1}-1\right)\right]{\theta }^2}\right\},\\ \mathrm{E}\left[\frac{{\partial }^2\mathit{ln}L}{\partial v^2}\right]={\sum }_i^{N}E\left(I_{\left\{y_i\ge 0\right\}}\right)\left\{-\frac{y_{i}p\left[2\left(p-1\right)\left(v-1\right)\theta +y_i\left(p(2\theta -1)-2\theta \right)\right]}{{\left(y_i+v-1\right)}^2{\alpha \left(i\right)}^2}-{\psi }_1\left(v\right)+{\psi }_1\left(v+y_i\right)\right\},\\ \mathrm{E}\left[\frac{{\partial }^2\mathit{ln}L}{{\partial \theta }^2}\right]={\sum }_i^{N}E\left(I_{\left\{y_i\ge 0\right\}}\right)\left\{-\frac{v}{{\left(\theta -1\right)}^2}-\frac{y_i}{{\theta }^2}\left[1+\frac{y_i{p}^2}{{\alpha \left(i\right)}^2}+\frac{2p}{\alpha \left(i\right)}\right]\right\},\end{array}$$

where ${\psi }_0$ is digamma function and ${\psi }_1$ is polygama function of order 1.

3.2.2. Generalized Likelihood Ratio Test

The generalized likelihood ratio (GLR) test statistic is given by

$$\mathrm{G}\mathrm{LR}=-2\ln \left(\frac{L(\widehat{\beta }*;x)}{L(\widehat{\beta };x)}\right)$$

where $\widehat{\beta }*$ is the vector of restricted ML estimators under $H_0$ and $\widehat{\beta }$ is the unrestricted ML estimators under $H_1$. The distribution of the test is ${\chi }^2$ with 1 degree of freedom.

3.3. Statistical Power Analysis of the Rao’s Score and Generalized Likelihood Ratio Tests

In this section, Table 2 displays the power of score and GLR tests from a Monte Carlo simulation study that is conducted with 1000 repetitions. In the simulation study, the significance level $\alpha$ is set at 5% and 10%. The sample sizes considered are N =100 (small), 500 (moderate) and 1000 (large), while different levels of under-, equi- and over-dispersion are incorporated. Rao’s score test and the GLR test are asymptotically equivalent.

Table 2 demonstrates that, in the scenario of equi-dispersion $\left(p=\sqrt{\nu }\theta /\left(1-\theta \right)\right)$, the estimated empirical levels of Rao’s score test and GLR test exhibit proximity to each other as the sample size N increases. In the case of under-dispersion, the GLR test consistently exhibits a higher power than Rao’s score test, even with small sample sizes. Meanwhile, a 100% detection rate is achieved for an effect size of 0.61 when the sample size $N\ge 500$. For over-dispersion, when the effect size is 0.25, Rao’s score test outperforms the GLR test marginally. However, as the effect size increases to 0.60, a slightly stronger detection is observed with the GLR test, particularly when performed with larger sample sizes ($N=1000$). It is evident that, as the deviation from $p=\sqrt{\nu }\theta /\left(1-\theta \right)$ increases, the power of the test also increases.

4. Modeling of Biological Count Data

The NB-shifted NB distribution has been fitted to count data for different indexes of dispersion. The three real-life data sets are as follows:

(1)

Frequency and distribution of chromosome aberrations (dicentrics plus ring) in peripheral blood lymphocytes irradiated in vitro with γ-rays (dose 10 in

Table 3. Frequency and distribution of chromosome aberrations (dicentrics plus ring) in peripheral blood lymphocytes irradiated in vitro with γ-rays [35].

No. of Cells	Observed Frequency	Expected Frequency
	Dose 10	GPD	GIT	COM– Poisson	NB-Shifted NB
		MLE	MLE	MLE	MLE	pgf-Estimator
0	0	1.62	1.49	1.26	0.07	0.12
1	9	8.50	8.24	8.01	7.57	8.44
2	26	21.54	21.52	21.73	23.56	24.39
3	33	35.03	35.49	35.92	37.97	37.53
4	39	41.08	41.60	41.77	41.94	40.48
5	36	37.02	37.10	36.97	35.63	34.30
6	26	26.66	26.35	26.18	24.81	24.29
7	23	15.77	15.42	15.36	14.74	14.95
8	3	7.81	7.65	7.65	7.68	8.20
9	2	3.28	3.29	3.30	3.58	4.09
10	2	1.18	1.25	1.25	1.52	1.88
11	1	0.50	0.62	0.60	0.91	1.33
Total	200
${\chi }^2$		10.67	10.60	10.54	9.85	9.86
p value		0.30	0.30	0.31	0.28	0.28
Parameter Estimates		$\widehat{\theta }$ = 4.82	${\widehat{p}}_2$ = 0.26	$\widehat{v}$ = 1.22	$\widehat{\theta }$ = 0.08	$\widehat{\theta }$ = 0.15
		$\widehat{\lambda }$ = −0.09	${\widehat{p}}_3$ = 0.13	$\widehat{\lambda }$ = 6.34	$\widehat{p}$ = 0.99	$\widehat{p}$ = 0.99
			n = 10		$\widehat{v}$ = 38.58	$\widehat{v}$ = 19.75
ID		0.85

Test for equi-dispersion: Rao’s score test statistic = 0.06 and GLR test statistic = 2.17 (Do not reject the null hypothesis).

and dose 6 in

Table 4. Frequency and distribution of chromosome aberrations (dicentrics plus ring) in peripheral blood lymphocytes irradiated in vitro with γ-rays [35].

No. of Cells	Observed Frequency		Expected Frequency
	Dose 6	Poisson	GPD	GIT	COM– Poisson	NB-Shifted NB
		MLE	MLE	MLE	MLE	MLE	pgf-Estimator
0	19	23.77	22.87	19.18	21.91	19.08	18.95
1	56	50.62	50.50	55.43	50.88	56.17	56.93
2	60	53.92	54.82	59.54	55.70	56.33	55.94
3	31	38.28	39.01	35.04	39.28	36.73	36.22
4	18	20.38	20.46	17.28	20.27	18.74	18.61
5	11	8.68	8.44	7.82	8.21	8.14	8.24
6	5	4.35	3.91	5.71	3.74	4.81	5.12
Total	200
${\chi }^2$		4.60	4.77	1.88	4.61	2.17	2.01
p value		0.47	0.31	0.76	0.33	0.54	0.57
Parameter Estimates		$\widehat{\lambda }$ = 2.13	$\widehat{\theta }$ = 2.17	${\widehat{p}}_2$ = 0.34	$\widehat{v}$ = 1.09	$\widehat{\theta }$ = 0.19	$\widehat{\theta }$ = 0.21
			$\widehat{\lambda }$ = −0.02	${\widehat{p}}_3$ = 0.35	$\widehat{\lambda }$ = 2.32	$\widehat{p}$ = 0.63	$\widehat{p}$ = 0.65
				n = 2		$\widehat{v}$ = 6.44	$\widehat{v}$ = 5.48
ID		0.97

Test for equi-dispersion: Rao’s score test statistic = 0.03 and GLR test statistic = 0.05 (Do not reject the null hypothesis).

) to represent under- and equi-dispersed data [35].

(2)

Fetal movement (in

Table 5. Fetal movement [36].

Number of Movements	Observed Frequency	Expected Frequency
	Number of Interval	GPD	GIT	COM– Poisson	NB-Shifted NB
		MLE	MLE	MLE	MLE	pgf-Estimator
0	182	182.50	176.61	176.67	182.02	181.96
1	41	39.49	46.71	46.63	41.22	41.60
2	12	11.62	12.29	12.30	10.30	10.14
3	2	3.95	3.23	3.24	3.78	3.69
4	2	1.46	0.85	0.85	1.52	1.49
5	0	0.57	0.22	0.23	0.65	0.63
6	0	0.23	0.06	0.06	0.28	0.27
7	1	0.17	0.02	0.02	0.23	0.23
Total	240
${\chi }^2$		6.09	48.75	48.32	4.73	4.87
p value		0.30	0.0	0.0	0.32	0.30
Parameter estimates		$\widehat{\theta }$ = 0.66	${\widehat{p}}_2$ = 0.001	$\widehat{v}$ = 0.001	$\widehat{\theta }$ = 0.50	$\widehat{\theta }$ = 0.50
		$\widehat{\lambda }$ = 0.22	${\widehat{p}}_3$ = 0.26	$\widehat{\lambda }$ = 0.26	$\widehat{p}$ = 0.08	$\widehat{p}$ = 0.09
			n = 1		$\widehat{v}$ = 0.28	$\widehat{v}$ = 0.27
ID		1.84

Test for equi-dispersion: Rao’s score test statistic = 26.85 and GLR test statistic = 19.67 (Reject the null hypothesis).

) to represent the over-dispersed data [36].

The performance of the NB-shifted NB distribution is compared with GPD, GIT and COM–Poisson distributions. The parameter n under GIT distribution is chosen based on the lowest chi-square values. For the data with ID of 0.97 (a value very close to 1), Poisson distribution fits well as a comparison study for equi-dispersed data. Under NB-shifted NB distribution conditions, both ML estimation (MLE) and pgf-based estimation are used.

For Table 3, the fits obtained by the NB-shifted NB distribution are comparable with the GPD, GIT and COM–Poisson distributions derived based on the p value. Both estimation methods work well, with similar chi-square values achieved. Meanwhile, Table 4 shows that, based upon the p value and chi-square value, the fits by the NB-shifted NB and GIT distributions are significantly better than the others. For this data set, pgf-based estimation is slightly better than MLE methods, with a lower chi-square value achieved. For over-dispersed data, the NB-shifted NB distribution achieves the highest p value (for MLE) and the lowest chi-square values if compared with others. The MLE method slightly outperforms pgf-based estimation, as presented in Table 5. In addition, the tests of equi-dispersion (Rao’s score and GLR test) indicate that the null hypothesis should be rejected in Table 5, which agrees with the computed index of dispersion. For the two data sets in Table 3 and Table 4, the simulation results presented in

Table 6. Summary statistics for MLE and pgf-estimator based on bootstrap simulation.

	Table 3
	MLE			pgf-Estimator
	$\widehat{\mathit{\theta }}$	$\widehat{\mathit{p}}$	$\widehat{\mathit{v}}$	$\widehat{\mathit{\theta }}$	$\widehat{\mathit{p}}$	$\widehat{\mathit{v}}$
Estimated parameters	0.08	0.99	38.58	0.15	0.99	19.75
Mean	0.09	0.98	57.32	0.14	0.98	36.97
Standard error	0.07	0.06	36.64	0.08	0.01	32.35
Confidence interval	(0.03, 0.24)	(0.83, 0.99)	(10.51, 100)	(0.03, 0.32)	(0.94, 0.99)	(7.41, 100)
Bias	0.01	−0.02	18.74	−0.01	−0.01	17.22
	Table 4
	MLE			pgf-Estimator
	$\widehat{\mathit{\theta }}$	$\widehat{\mathit{p}}$	$\widehat{\mathit{v}}$	$\widehat{\mathit{\theta }}$	$\widehat{\mathit{p}}$	$\widehat{\mathit{v}}$
Estimated parameters	0.19	0.63	6.44	0.21	0.65	5.48
Mean	0.17	0.59	20.50	0.20	0.61	18.35
Standard error	0.10	0.19	30.70	0.11	0.16	28.73
Confidence interval	(0.02, 0.36)	(0.01, 0.81)	(2.41, 100)	(0.02, 0.40)	(0.08, 0.81)	(2.16, 99.99)
Bias	−0.02	−0.05	14.06	−0.02	−0.04	12.88
	Table 5
	MLE			pgf-Estimator
	$\widehat{\mathit{\theta }}$	$\widehat{\mathit{p}}$	$\widehat{\mathit{v}}$	$\widehat{\mathit{\theta }}$	$\widehat{\mathit{p}}$	$\widehat{\mathit{v}}$
Estimated parameters	0.5	0.08	0.28	0.5	0.09	0.27
Mean	0.48	0.07	0.40	0.49	0.08	0.37
Standard error	0.14	0.05	0.36	0.14	0.06	0.28
Confidence interval	(0.20, 0.72)	(0.01,0.17)	(0.08, 1.28)	(0.22, 0.74)	(0.01, 0.18)	(0.07, 1.10)
Bias	−0.02	−0.01	0.12	−0.01	−0.01	0.10

show that the bias for the estimated parameter $\widehat{v}$ is relatively large.

To evaluate the performance of the pgf-estimator in fitting the NB-shifted NB distribution to real-life data sets, a non-parametric bootstrap simulation [37] was conducted with 1000 repetitions. The simulation involved re-sampling, with replacement, from the original data set to generate a multitude of bootstrap samples. For each of these samples, the maximum likelihood estimate was computed. Non-parametric bootstrap simulation ensures that these insights are obtained without making specific assumptions about the underlying population’s distribution.

The results of this simulation are summarized in Table 6, which provides key statistics the mean, standard error, and confidence interval from the bootstrap analysis for both the MLE and the pgf-estimator. Additionally, the bias of the MLE and pgf-estimator was assessed by comparing the average of the bootstrap estimates to their respective estimated parameters. The findings presented in Table 6 indicate that pgf-estimator has a lower bias than MLE. The mean, standard error, confidence intervals, and bias computed for both estimators are in close agreement. For the two data sets in Table 3 and Table 4, the standard errors for the estimated parameter $\widehat{v}$ are relatively large, indicating low precision in the estimates.

5. Concluding Remarks

A two-component mixture model of a distribution with its shifted counterpart has been proposed in order to generate models that can cater for under-, equi- and over-dispersion. As an example, the paper considers the finite mixture of NB and shifted NB distributions. The distribution does not, like the generalized Poisson distribution, have a problem with the range of the parameters. Even though the COM–Poisson distribution has a simple pmf, the computation of the normalizing constant could encounter problems for extreme values of the parameters. The GIT distribution does not face these issues. However, in comparison, the finite mixture of NB distributions (NB-shifted NB) is much simpler. The proposed distribution has been fitted to a variety of biological data sets. The good fits, shown by the distribution relative to the GPD, COM–Poisson distribution and GIT distribution, proves that it is a viable and simple alternate model for under-, equi- and over-dispersion. The use of components other than NB and shifted NB for the finite mixture will be considered elsewhere.