Two new bivariate zero-inflated generalized Poisson distributions with a flexible correlation structure

To model correlated bivariate count data with extra zero observations, this paper proposes two new bivariate zero-inflated generalized Poisson (ZIGP) distributions by incorporating a multiplicative factor (or dependency parameter) λ, named as Type I and Type II bivariate ZIGPλ distributions, respectively. The proposed distributions possess a flexible correlation structure and can be used to fit either positively or negatively correlated and either overor under-dispersed count data, comparing to the existing models that can only fit positively correlated count data with over-dispersion. The two marginal distributions of Type I bivariate ZIGPλ share a common parameter of zero inflation while the two marginal distributions of Type II bivariate ZIGPλ have their own parameters of zero inflation, resulting in a much wider range of applications. The important distributional properties are explored and some useful statistical inference methods including maximum likelihood estimations of parameters, standard errors estimation, bootstrap confidence intervals and related testing hypotheses are developed for the two distributions. A real data are thoroughly analyzed by using the proposed distributions and statistical methods. Several simulation studies are conducted to evaluate the performance of the proposed methods.


Introduction
As the simplest distribution for modeling count data, the Poisson distribution possesses an exceptional property of equi-dispersion, i.e., its mean and variance are identical.When the observed variance is larger than the theoretical variance, over-dispersion has occurred and it is a very common feature in applications because in practice populations are frequently heterogeneous.In such situation, a useful alternative such as the negative-binomial model with an additional free parameter may provide a better fit.Conversely, under-dispersion means that there is less variation in the data than predicted.To model count data with both over-dispersion and/or under-dispersion, the generalized Poisson (GP) distribution can be used and its probability mass function (pmf) with two parameters is defined by ( [7], p.5) x = 0, 1, 2, . . .,
To model two-dimensional correlated count data, the bivariate Poisson, bivariate negative-binomial and bivariate generalized Poisson models have been developed.By using the method of trivariate reduction, Holgate [15] first proposed the bivariate Poisson distribution and Famoye & Consul [9] then extended it to the more general bivariate generalized Poisson distribution.A major disadvantage for the method of trivariate reduction is that so-generated bivariate distributions can be used to model count data with positive correlation only.
However, two-dimensional count data with negative correlation are often encountered in practice.For example, Table 1 below gives a paired count data of twelve patients who experienced frequent premature ventricular contractions (PVCs) and were administered a drug with antiarrhythmic properties.One-minute EKG recordings were taken before and after the drug administration as shown in the second and third rows of Table 1.To describe the data, we calculated the sample means and variances as x = 13.4167,s 2 x = 166.9924,ȳ = 1.9167, s 2 y = 14.4470 and the sample correlation coefficient as r = −0.2305847.Both X and Y are obviously over-dispersed and they are negatively correlated.Thus, the aforementioned distributions are not appropriate to fit this data set.

Table 1
The PVC counts for twelve patients one minute after administering a drug with antiarrhythmic properties ( Lakshminarayana et al. [20] first proposed a new type of bivariate Poisson distribution, whose joint pmf was expressed as a product of two Poisson margins with a multiplicative factor.Subsequently, a new bivariate negativebinomial regression model and a new bivariate GP distribution are extended in a similar way by Famoye [12,13].This kind of expression not only ensures the marginal distribution as the traditional distribution (e.g., Poisson, generalized Poisson, negative-binomial) but also admits positive or negative correlation.The general form of the joint bivariate pmf of (X, Y ) is p(x, y) = p 1 (x)p 2 (y) where g 1 (x) and g 2 (y) are two bounded functions, and λ is called the multiplicative factor or dependency parameter.
It is easy to verify that p 1 (x) and p 2 (y) are two marginal pmfs.When there is a larger frequency of (0,0) observations in the bivariate count data, the zero-inflation phenomenon arises.For example, there are seven zero values of Y out of a total of 12 observations in Table 1.The multivariate zero-inflated Poisson models were developed by using the mixture technique [21] and by using the stochastic representation technique [22].Their models can only fit over-dispersed data with positive correlation.Krishna & Tukaram [19] provided a bivariate power series distribution that utilized the multiplicative factor form and incorporated inflation at zero.Bivariate zero-inflated Poisson distribution were then discussed as a special case of power series distributions and the associated properties were explored.The univariate marginally follows zeroinflated Poisson distribution which is still over-dispersed, and the correlation can have either sign.All these zeroinflated models cannot comprehensively address the correlation and dispersion issues.
The aim of this paper is to develop two new bivariate zero-inflated generalized Poisson (ZIGP) distributions by incorporating a multiplicative factor (or dependency parameter) λ, named as Type I and Type II bivariate 107 ZIGP λ distributions, respectively.The two distributions have the following features: (i) Each marginal follows a univariate ZIGP distribution, which is a mixture of the degenerate distribution with all mass at (0, 0) point and a GP distribution.A ZIGP distribution can model zero-inflated count data with both over-dispersion and/or underdispersion.(ii) They possess a flexible correlation structure and can be used to fit either positively or negatively correlated zero-inflated count data, comparing to the existing models that can only fit positively correlated zeroinflated count data.(iii) The two marginal distributions of Type I bivariate ZIGP λ share a common parameter of zero inflation while the two marginal distributions of Type II bivariate ZIGP λ have their own parameters of zero inflation, resulting in a much wider range of applications.
The rest of the paper is organized as follows.In Section 2, we propose a new bivariate ZIGP distribution indexed by a multiplicative factor via a mixture of a Bernoulli variable with the bivariate GP distribution of Famoye [13], provide some important distributional properties, maximum likelihood estimations of parameters, standard errors estimation, bootstrap confidence intervals and related testing hypotheses.In Section 3, the Type II bivariate zeroinflated generalized Poisson distribution with a multiplicative factor is also developed.The proposed distributions and statistical methods are applied to analyze the Australian health care utilization data in Section 4. In Section 5, several simulation studies are conducted to evaluate the performance of the proposed methods.Section 6 presents a discussion.Some technical details are put in the Appendix.

Type I bivariate zero-inflated generalized Poisson distribution with a multiplicative factor
Famoye [13] introduced a new bivariate GP distribution indexed by an unknown real number λ called the multiplicative factor.Its joint pmf Pr( for x 1 , x 2 = 0, 1, 2, . . ., where We denote this new bivariate GP distribution by It is easy to verify that X i ∼ GP (θ i , α i ) for i = 1, 2, which are irrelevant to the multiplicative factor λ. When α 1 = α 2 = 0, this bivariate GP distribution reduces to the bivariate Poisson distribution with the multiplicative factor λ ∈ R developed by Lakshminarayana et al. [20].
When there is a larger frequency of (0,0) observations in the bivariate count data, the issue of the over-dispersion may arise.To model such extra zero points in the data, we propose a new bivariate ZIGP distribution indexed by a multiplicative factor λ ∈ R via a mixture of a Bernoulli variable with the bivariate GP distribution (1).
is said to follow the Type I bivariate ZIGP distribution indexed by a multiplicative factor λ, denoted by , where ϕ ∈ [0, 1) and c i (i = 1, 2) are defined in (1), the range of λ is to be discussed in Section 2.1.It is easy to derive the joint pmf , we discuss the following three special cases:

(c)
When α 1 = α 2 = λ = 0, it reduces to the Type I bivariate ZIP distribution, denoted by , which is a special case of the Type I multivariate ZIP distribution recently developed by Liu & Tian [22].

Marginal distributions, moments and correlation
From (2), we have Y i = ZX i following an univariate zero-inflated generalized Poisson distribution, denoted by The pmf, expectation and variance of Y i are given by respectively.It is shown that Y 1 and Y 2 can be either over-or under-dispersed.The covariance between Y 1 and Y 2 is where From (2), we can see that the dependency structure between Y 1 and Y 2 hinges on the common Bernoulli variable Z and the correlation between X 1 and X 2 .In fact, the correlation coefficient of Y 1 and Y 2 is which could be positive, zero or negative depending on the values of parameters ϕ and λ. ) respectively, where and .

Maximum likelihood estimations of parameters
λ (ϕ; θ 1 , θ 2 , α 1 , α 2 ) for j = 1, . . ., n, and Y obs = {(y 1j , y 2j ) ⊤ : j = 1, . . ., n} denote the observation data.Let J = {j: (y 1j , y 2j ) ⊤ = 0 0, j = 1, . . ., n} and m 0 denote the number of elements in J. Let θ = (θ 1 , θ 2 , α 1 , α 2 , λ) ⊤ , then the observed-data likelihood function is To obtain the MLEs of the parameters, we employ the EM algorithm.A latent variable Z is introduced to split m 0 into Z and m 0 − Z, so that the conditional predictive distribution of Z given Y obs and (ϕ, θ) is On the other hand, the complete-data likelihood function is proportional to where L 1 (ϕ|Y com ) = ϕ z (1 − ϕ) n−z only involves ϕ and } only involves θ.Obviously, the complete-data MLE of ϕ is φ = z/n.The logarithm of L 2 (θ|Y com ) is given by where Since the complete-data MLEs of θ are not available in closed form, we adopt the Newton-Raphson algorithm to calculate the complete-data MLEs of θ.The score vector and the Hessian matrix are derived in Appendix A. Thus, the M-step is to calculate The E-step is to replace z (t) in ( 11) by the conditional expectation where c

Estimation of standard errors
In this subsection, we are interested in deriving the standard deviations of the MLEs φ and θ of the parameters.Let ℓ(ϕ, θ|Y com ) = log L(ϕ, θ|Y com ) denote the logarithm of the complete-data likelihood function (8).Then, its first and second derivatives are given by ∇ℓ 2 (θ|Y com ) and ∇ 2 ℓ 2 (θ|Y com ) are given by (10).
According to Louis [23], the observed information I( φ, θ|Y obs ) can be expressed as where The computation of the conditional expectations in (13) involves the conditional predictive distribution (7).We have .
The estimated standard errors are the square roots of the diagonal elements of the inverse observed information matrix I −1 ( φ, θ|Y obs ).Finally, we can use these estimated standard errors to construct the (1 − α)100% asymptotic Wald confidence interval (CIs) for the parameters (ϕ, θ).
Alternatively, when the observed information matrix is too complicated, we may use the square roots of the diagonal elements of the inverse complete information matrix to approximate the estimated standard errors.The complete information matrix is defined by

Bootstrap confidence intervals
The Wald CIs are reliable only for large sample sizes and the lower/upper bound may exceed the natural bound of parameters.For example, the Wald CI of ϕ may fall outside the unit interval (0, 1]. The bootstrap method is useful in calculating the bootstrap CIs, which are reliable for small to moderate sample sizes.For an arbitrary function of (ϕ, θ), say ϑ = h(ϕ, θ), the bootstrap method can be used to obtain the bootstrap CI of ϑ.Let ( φ, θ) denote the MLEs of (ϕ, θ) obtained by the EM algorithm embedded with the Newton-Raphson algorithm specified by ( 11)-( 12), then θ = h( φ, θ) is the MLE of ϑ.Based on the obtained MLEs ( φ, θ) = ( φ, θ1 , θ2 , α1 , α2 , λ), we can generate for j = 1, . . ., n using the conditional sampling technique presented in Section 2.2.Based on the bootstrap sample , we first compute the MLEs ( φ * , θ * ) and then obtain a bootstrap replication θ * = h( φ * , θ * ).Independently repeating the above process G times, we obtain G bootstrap samples . The standard error of θ, se( θ), can be estimated by the sample standard deviation of the G replications, i.e., se where θL and θU are the 100(α/2) and 100(1 , respectively.Note that the CI in ( 16) is based on the normality assumption, while (17) provide accurate confidence intervals without making the normality assumption.However, (17) does not adjust the confidence interval to account for skewness in the underlying population or other errors that can result in where θ is not the sample mean.Thus, we consider the bootstrap percentile-t CI, which can adjust such errors.Based on the generated G bootstrap samples where tα/2 and t1−α/2 are the 100(α/2) and 100(1 , respectively.

Testing hypotheses
2.6.1.Likelihood ratio test for zero inflation.Suppose we want to test the null hypothesis Under H 0 , the likelihood ratio test (LRT) statistic where (0, θH0 ) denote the constrained MLEs of (ϕ, θ) under H 0 and ( φ, θ) denote the unconstrained MLEs of (ϕ, θ), which are obtained by the algorithm specified in ( 11)- (12).Since the null hypothesis in (19) corresponds to ϕ being on the boundary of the parameter space and the appropriate null distribution is a 50:50 mixture of χ 2 (0) (i.e., the degenerate distribution with all mass at zero) and χ 2 (1), see Self & Liang [25] and Feng & McCulloch [14].Hence, the p-value ([16], p.78; [17], p.225) is 2.6.2.LRT for independency parameter.If λ = 0, the Type I bivariate ZIGP random vector reduces to the bivariate random vector that is the product of a common zero-inflated Bernoulli variable Z and a random vector with two independent GP random components.Correlation between Y 1 and Y 2 in (4) becomes which must be non-negative.Therefore, we would like to test The corresponding LRT statistic is given by Under H 0 , T 2 approximately follows the χ 2 (1) distribution.The corresponding p-value is 2.6.3.LRT of Type I bivariate ZIGP λ against Type I bivariate ZIP λ .
We have mentioned that when the parameters α 1 = α 2 = 0, the Type I bivariate ZIGP distribution indexed by the multiplicative factor λ reduces to the Type I bivariate ZIP distribution indexed by the multiplicative factor λ. We test the following null hypothesis To assess the adequacy of the Type I bivariate ZIGP over the Type I bivariate ZIP.The LRT statistic is Note that under H 0 , the MLEs of parameters in ZIP (I) λ (ϕ; θ 1 , θ 2 ) can be obtained through an EM algorithm embedded with a Newton-Raphson algorithm in a similar way as that for the Type I bivariate ZIGP distribution.Under H 0 , the test statistic T 3 follows χ 2 (2), thus the corresponding p-value is

Type II bivariate zero-inflated generalized Poisson distribution with a multiplicative factor
According to (3), two marginal distributions of the Type I bivariate ZIGP distribution with a multiplicative factor are zero-inflated GP and are forced into sharing a common parameter of zero inflation.This limitation evokes us to extend the Type I bivariate ZIGP to a Type II bivariate ZIGP distribution with a multiplicative factor in which two different parameters of zero inflation are allowed.
is said to follow the Type II bivariate ZIGP distribution indexed by a multiplicative factor λ, denoted by if In particular, ϕ 1 = 0 indicates that only the Y 2 marginally follows a zero-inflated GP distribution.Similarly, ϕ 2 = 0 indicates that only the Y 1 marginally follows a zero-inflated GP distribution.Similar to Section 2.2, we can generate i.i.d.random vectors from (28) by using the conditional sampling technique.

MLEs of parameters
Let θ = (θ 1 , θ 2 , α 1 , α 2 , λ) ⊤ , the observed-data likelihood function is We adopt the EM algorithm to obtain the MLEs of parameters.First we augment Y obs with latent variables The conditional predictive distributions of these latent variables are given by where The complete data are denoted by Y com = {Y obs , u 0 , u 1 , u 2 , u 3 , {w j } j∈J1 , {z j } j∈J2 }.Thus, the complete-data likelihood function is proportional to where only involving θ.Since the complete-data MLEs of ϕ 1 and ϕ 2 are given by where u + = ∑ 3 i=0 u i .On the other hand, where ∑ n j=1 y 2j .Since the complete-data MLEs of θ are not available in closed form, we adopt the Newton-Raphson algorithm to calculate the complete-data MLEs of θ.The score vector and the Hessian matrix are derived in Appendix B. The M-step is to update = u where and u i .The E-step is to replace (u 0 , u 1 , u 2 , u 3 ), {w j } j∈J1 and {z j } j∈J2 in (34) by their conditional expectations Note that if the initial value of ϕ 1 is set to be ϕ (0) 1 = 0, the above algorithm yields the estimates of parameters, where that only Y 2 is marginally inflated; if the initial value of ϕ 2 is set to be ϕ (0) 2 = 0, it yields the estimates of parameters, where only Y 1 is marginally inflated.
The estimation of standard errors can be obtained via Louis's method, which is very similar to that stated in Section 2.4.The bootstrap method presented in Section 2.5 can be used to construct bootstrap CIs of parameters for the current situation.

LRT for zero inflation.
First, we want to test whether there exists zero inflation in both marginal distributions.Thus, we consider the following null and alternative hypotheses: Stat., Optim.Inf.Comput.Vol.

Applications in Australian health care utilization data
Cameron & Trivedi [5] reported data concerning the demand for health care in Australia which refers to the Australian Health Survey for 1977-1978.Let Y 1 denote the number of consultations with a doctor or a specialist and Y 2 denote the total number of prescribed medications used in past two days.The data are given in Table 2.
Before the data analysis, we first examine the descriptive statistics for the data: ȳ1 = 0.3017341, ȳ2 = 0.8626204; s 2 1 = 0.6370176, s 2 2 = 2.0032855; the sample correlation coefficient r = 0.3077787.We have observed that both marginal count data exhibit over-dispersion and very high frequencies of zero observations since most of the observed frequencies fall in the (0,0) category.The positive sample correlation coefficient r indicates that there is a positive correlation between Y 1 and Y 2 .Based on these observations, we first try to fit the bivariate count data by the Type I bivariate ZIGP model with a multiplicative factor in Section 2.
Besides, we adopt the bootstrap method to compute the bootstrap CIs for these parameters.Using the obtained MLEs, we generate G = 6,000 bootstrap samples and the three types of bootstrap CIs ( 16)-( 18) are summarized in Table 4. std: Calculated as the square roots of the diagonal elements of I −1 ( φ, θ|Y obs ), cf.(13).
Suppose that we want to test the null hypothesis H 0 : ϕ = 0 against the alternative hypothesis H 1 : ϕ > 0. According to (20), we calculate the value of the LRT statistic which is given by t 1 = 20.2776,then from (21), we have the corresponding p v1 = 3.3491 × 10 −6 < α = 0.05.Thus, the null hypothesis should be rejected at 5% significance level.
If we want to test the null hypothesis H 0 : λ = 0 against the alternative hypothesis H 1 : λ ̸ = 0.According to (23), the value of the LRT statistic is t 2 = 88.9215 and from (24) the p-value is p v2 < 0.0001.Thus, we should reject the null hypothesis at the significance level of 5%.
If we want to test the null hypothesis H 0 : α 1 = α 2 = 0.According to (26), the value of the LRT statistic is t 3 = 1615.4960and from (27), we have p v3 ≪ 0.0001.Thus, the null hypothesis is rejected at the significance level of 5%.

Marginal analysis
The sample correlation coefficient r = 0.3077787, indicating that there is a positive correlation between Y 1 and Y 2 .Thus, any two independent univariate distributions cannot be used to model the health care utilization data.Besides, since both marginal count data are over-dispersed according to their marginal sample means and variances, the bivariate Poisson distribution constructed by the trivariate reduction method is also not appropriate due to the property of equi-dispersion.Note that the over-dispersion may result from the excess (0,0) points in the observations, the Type I multivariate ZIP distribution proposed by Liu & Tian [22] could be considered.If we fit the data by the Type I bivariate ZIP distribution, denoted by (Y 1 , Y 2 ) ⊤ ∼ ZIP (I) (ϕ; θ 1 , θ 2 ), the MLEs of three parameters are φ = 0.483007, θ1 = 0.583633, θ2 = 1.668533.The estimated correlation coefficient is which is very close to the sample correlation coefficient r = 0.3077787.From the definition of (Y 1 , Y 2 ) ⊤ ∼ ZIP (I) (ϕ; θ 1 , θ 2 ), we have Y i ∼ ZIP (ϕ; θ i ) for i = 1, 2, indicating that the marginal distributions must share a common zero inflation parameter ϕ.On the other hand, based on the marginal counts, the univariate ZIP distribution of Y 1 and Y 2 are estimated to be ZIP ( φM 1 = 0.650393; θM 1 = 0.863068) and ZIP ( φM 2 = 0.510188; θM 2 = 1.761464), respectively.It is clear that the difference between φM 1 and φM 2 is slight.An alternative to the above ZIP (I) (ϕ; θ 1 , θ 2 ) is the Type I bivariate ZIP distribution indexed by the multiplicative factor λ, i.e., ZIP (I) λ (ϕ; θ 1 , θ 2 ), which is reduced from the proposed Type I bivariate ZIGP distribution with the same multiplicative factor λ by setting α 1 = α 2 = 0.The corresponding MLEs of the four parameters are given by φ = 0.489964, θ1 = 0.588041, θ2 = 1.685382, λ = −0.477044.From (5), the estimated correlation coefficient is Corr(Y 1 , Y 2 ) = 0.2886614, which is close to the sample correlation coefficient r = 0.3077787.Note that the distribution ZIP (I) λ (ϕ; θ 1 , θ 2 ) has the same marginal distributions ZIP (ϕ; θ i ) for i = 1, 2 as ZIP (I) (ϕ; θ 1 , θ 2 ).For the proposed Type I bivariate ZIGP distribution with the multiplicative factor λ, according to (3), we have That is, the two marginal distributions should share a common zero inflation parameter ϕ.However, based on the marginal counts, the univariate ZIGP distribution of Y 1 and Y 2 are estimated to be and Since there is a larger difference between ϕ M 1 and ϕ M 2 , from the viewpoint of marginal analysis, Type I bivariate ZIGP distribution with the multiplicative factor λ is not appropriate to fit the data.
If we want to test the null hypothesis H 0 : λ = 0 against the alternative hypothesis H 1 : λ ̸ = 0.According to (43), the value of the LRT statistic is t 7 = 585.2057and from (44) we have the corresponding p-value is less than 0.0001.Thus, we should reject the null hypothesis.
The marginal analysis based on the Type II bivariate ZIGP λ is the same as that based on the Type I bivariate ZIGP λ , which is still given by ( 45) and (46).

Model comparison
We adopt the Akaike information criterion (AIC; [2]) and Bayesian information criterion (BIC; [24]) to compare six different models: (a) Karlis & Ntzoufras [18] considered two regression models to fit the data; i.e., bivariate Poisson (BP) constructed by the method of trivariate reduction based on three independent Poisson(θ i ) variables for i = 0, 1, 2, and (b) diagonal inflated bivariate Poisson (DIBP), which led to a model with only (0, 0) inflation but without (1,1) inflation in the current case.For the purpose of comparison, we concentrate our attention upon the situation of without covariates.
(c) Liu & Tian [22]  The estimates of parameters for first four models are reported in the second column of Table 7.The values of AIC and BIC based on the Total Likelihood Function are summarized in the last two columns of Table 7. From Table 7, we can see that the Type II bivariate ZIGP λ is the best model in terms of the log-likelihood, AIC and BIC.It also coincides with the result from the marginal analysis.On the other hand, we can compare the performances of these models via the estimated mean, variance and correlation coefficient which are reported in Table 8.
From Table 8, we can see that both the Type I and Type II bivariate ZIGP λ distributions perform better at adjusting the over-dispersion of the data, while the Type I bivariate ZIP and Type I bivariate ZIP λ have better estimates on 124 BIVARIATE ZIGP DISTRIBUTION WITH A FLEXIBLE CORRELATION STRUCTURE the correlation coefficient.The key point lies on that the GP distribution in the ZIGP λ has one more parameter than the Poisson distribution, thus it can better accommodate the degree of over-dispersion and zero-inflation in the data.Moreover, the GP distribution can also model under-dispersed data depending on the values of the two parameters.Besides, the addition of the multiplicative factor λ in the pmf (1) can further improve the model fitting so that it can fit either positively or negatively correlated count data.

Simulation studies
To evaluate the performance of the proposed methods for the two types of bivariate ZIGP λ distributions, we first investigate the accuracy of point estimators and interval estimators.Next, we focus on the type I error rates and powers of the LRT in the Type I and Type II bivariate ZIGP λ distributions, respectively.

Accuracy of point estimators and interval estimators
In this subsection, we assess the accuracy of the point estimators and interval estimators in the two types of bivariate ZIGP λ distributions.We consider one case for each type.For the Type I bivariate ZIGP λ distribution, the true values of the parameters (ϕ; θ 1 , θ 2 , α 1 , α 2 , λ) are set to be (0.4; 0.5, 0.3, 0.6, 0.2, 2).Based on (4), the estimated correlation coefficient between the two random variables Y 1 and Y 2 is 0.276215; i.e., they are positively correlated.
For the Type II bivariate ZIGP λ distribution, we consider the situation with zero inflation only in one margin, and the true values of the parameters (ϕ 1 , ϕ 2 ; θ 1 , θ 2 , α 1 , α 2 , λ) are set to be (0, 0.4; 1, 0.6, 0.3, 0.4, −2).Based on the first formula of (31), the estimated correlation coefficient between Y 1 and Y 2 is −0.1513673, indicating negative correlation.Sample size is set to be n = 1, 000.In the first case, we generate for j = 1, . . ., n, then calculate the MLEs via the EM algorithm embedded with the Newton-Raphson algorithm specified by ( 11)-( 12) and 95% bootstrap CIs with G = 5, 000 bootstrap replications.In the second case, we generate for j = 1, . . ., n based on the given parameter configuration and calculate the MLEs via the EM algorithm embedded with the Newton-Raphson algorithm specified by (34)-(35) and 95% bootstrap CIs with G = 5, 000 bootstrap replications.Next, we independently repeat each process for 1000 times, the resultant mean of the MLEs (denoted by MLE), the mean squared errors of the estimates (denoted by MSE), the average width of confidence intervals for parameters (denoted by Width) and the average coverage proportion of bootstrap CIs (denoted by CP) are reported in Table 9, respectively.

Performance of the LRT in Type I bivariate ZIGP λ
In this subsection, we investigate the performance of the LRT for testing hypotheses (19) and (22) for the Type I bivariate ZIGP λ distribution.We conduct extensive simulations to observe the changes of type I error rates and powers against the sample sizes which are set to be n = 500(50)950.

LRT for testing H
To estimate the type I error rates (with H 0 : ϕ = 0) and powers (with H 1 : ϕ > 0) for various sample sizes, where the values of ϕ in H 1 are chosen to be 0.01 and 0.05.Given a combination of parameters (n, ϕ, θ 1 = 1, θ 2 = 0.6, α 1 = 0.3, α 2 = 0.4, λ = 2), we generate for l = 1, . . ., L (L = 5, 000).For each group of samples , we conduct the testing hypothesis.Let r 1 denote the number of rejecting the null hypothesis H 0 : ϕ = 0 by the test statistic T 1 given by (20).Then the actual significance level can be estimated by r 1 /L with ϕ = 0 and the power of the test statistic T 1 can be estimated by r 1 /L with ϕ > 0.
The empirical levels/powers of the LRT statistic T 1 are summarized in Table 10. Figure 3 plots the type I error rates and powers of the LRT in testing H 0 : ϕ = 0 against H 1 : ϕ > 0 with two different values of ϕ > 0 for various sample sizes.12) and ( 34)-(35), respectively; MSE is equal to the sum of the variance and the squared bias of the estimator; Width and CP are the average width and coverage proportion of 1000 bootstrap confidence intervals.for l = 1, . . ., L (L = 5, 000).For each group of samples , we conduct the testing hypothesis.Let r 2 denote the number of rejecting the null hypothesis H 0 : λ = 0 by the test statistic T 2 given by (23).Then the actual significance level can be estimated by r 2 /L with λ = 0 and the power of the test statistic T 2 can be estimated by r 2 /L with λ ̸ = 0.
The empirical levels/powers of the LRT statistic T 2 are summarized in Table 11 and Figure 4 shows the changes of the type I error rates and powers of the LRT in testing H 0 : λ = 0 against H 1 : λ ̸ = 0 with three different values of λ ̸ = 0 for various sample sizes.

Performance of the LRT in Type II bivariate ZIGP λ
In this subsection, we investigate the performance of the LRT for the testing hypotheses (36) and (39) for the Type II bivariate ZIGP λ distribution.We conduct extensive simulations to observe the changes of type I error rates and powers of the LRT.The sample sizes are set to be n = 500(50)950.
The empirical levels of the LRT statistic T 4 are summarized in Table 12 and Figure 5 plots the type I error rates of the LRT in testing H 0 : (ϕ 1 , ϕ 2 ) = (0, 0) against H 1 : (ϕ 1 , ϕ 2 ) ̸ = (0, 0) for different sample sizes.The empirical levels/powers of the LRT statistic T 6 are summarized in Table 13 and Figure 6 portrays the trend of the type I error rates and powers of the LRT in testing H 0 : ϕ 2 = 0 against H 1 : ϕ 2 > 0 with three different values of ϕ 2 > 0 for various sample sizes.
From Figures 3-6, we can see that the lines for empirical levels of these tests all fluctuate near the line of α = 0.05, indicating that they perform well in controlling the type I error rates around the pre-chosen nominal level.The curves for empirical powers under different situations are non-decreasing in general although including some downs, which still shows that the tests tend to be more powerful as the sample size n becomes larger.

Discussion
This paper extended the new bivariate generalized Poisson distribution [13] to zero-inflated situation that can model bivariate count data with zero-inflated marginal distributions.An important characteristic is that it can model negatively correlated data while most existing models cannot.It retains the flexibility of GP distribution in allowing either over-or under-dispersion and accepts either positive or negative correlation from data.A common parameter or two different parameters on zero inflation of the margins are both considered, including only one marginal zero inflation in the model as a special case.However, due to the limitation that it can only apply to the bivariate case, zero-inflated multivariate GP model are desired on the basis of the new multivariate GP distribution [1] in which the correlation can also be positive, zero or negative.Another problem on the function in the multiplicative factor are all chosen to be g(y) = e −y when setting distribution, see [20,12,13,1].While no evidence is provided to show it outperformed others.Actually, the constructed structure in the expression of probability mass function represents a group of distributions.In the future, we may extend the proposed distribution to parametric model by considering covariates and semi-parametric random effects model ( [3]).
Appendix A: Derivation of (10) Since c i = exp[θ i (s i − 1)] and log s i − α i θ i (s i − 1) + 1 = 0 for i = 1, 2, the related partial derivatives are given by First we set The first and second partial derivatives in (10) are given by the following expressions: and where ℓ 2 = ℓ 2 (θ|Y com ) which is given in (9).

Figure 3 .
Figure 3. (i) The type I error rates for testing H 0 : ϕ = 0 against H 1 : ϕ > 0 in the Type I bivariate ZIGP λ distribution and the dashed line is set as the predetermined significance level of α = 0.05; (ii) the powers when ϕ = 0.01 in H 1 ; (iii) the powers when ϕ = 0.05 in H 1 .

Figure 4 .
Figure 4. (i) The type I error rates for testing H 0 : λ = 0 against H 1 : λ ̸ = 0 in the Type I bivariate ZIGP λ distribution and the dashed line is set as the predetermined significance level of α = 0.05; (ii) the powers when λ = −1 in H 1 ; (iii) the powers when λ = 1 in H 1 ; (iv) the powers when λ = 2 in H 1 .

Table 2
Cross tabulation of the health care utilization data in the Australian Health Survey for 1977-1978(Cameron and

Table 4
Three 95% boostrap CIs for the Australian Health Survey data (Type I) (15)The sample standard deviation based on the bootstrap samples, cf.(15).

Table 5
MLEs and CIs of parameters for the Australian Health Survey data (Type II) (14) Calculated as the square roots of the diagonal elements of I −1 ( φ2 , θ|Ycom), cf.(14).

Table 6
Two 95% bootstrap CIs for the Australian Health Survey data (Type II) (15)B : The sample standard deviation based on the bootstrap samples, cf.(15).
recently proposed the Type I multivariate ZIP distribution to fit the data in which Y 1 and Y 2 are correlated just through a common zero inflation factor., θ 2 ), see the third paragraph of Section 4.2.(e) Type I bivariate ZIGP λ , see Section 4.1.(f) Type II bivariate ZIGP λ , see Section 4.3.

Table 8
Comparison of estimated mean, variance and correlation coefficient for four different models

Table 9
MLEs and bootstrap confidence intervals of parameters NOTE: MLE is the mean of the 1000 point estimates via the EM algorithm embedded with the Newton-Raphson algorithm (11)-(

Table 13
Empirical levels/powers of the LRT statistic T 6 based on L = 1,000 replications