On bivariate Poisson regression models

1 Introduction

Joint modeling of two or more counts data has received a great deal of attention in recent years. Bivariate count models are used in cases where two count variables are correlated and need to be jointly estimated. For example, variables can include the number of emergency and non-emergency visits by a person to the hospital, the number of insurance claims with and without bodily injuries, or the number of voluntary and involuntary job changes.

The bivariate Poisson is the most widely used model for bivariate counts. It was proposed by Holgate (1964) and presented by Johnson and Kotz (1969). The definition of the bivariate Poisson distribution is not unique. Several approaches have been discussed by Kocherlakota and Kocherlakota (1992). Here, we adopt the trivariate reduction method to construct the distribution (Johnson et al., 1997). Consider three independent Poisson random variables, $X_k\text{,}k=1\text{,}2\text{,}3$ , with parameters (means) ${\lambda }_k>0\text{,}k=1\text{,}2\text{,}3$ . Then, the random variables $Y_1=X_1+X_3$ and $Y_2=X_2+X_3$ are set to follow joint bivariate Poisson distribution JBP¹(λ₁, λ₂, λ₃). The joint probability mass function (p.m.f) is given by

It can be easily shown that $Y_1$ and $Y_2$ are marginally distributed as Poisson with means ${\lambda }_1+{\lambda }_3$ and ${\lambda }_2+{\lambda }_3$ , respectively. The covariance of $Y_1$ and $Y_2$ is ${\lambda }_3$ , and hence, ${\lambda }_3$ is a measure of dependence between the two random variables. If ${\lambda }_3=0$ , then the two variables are independent and the bivariate Poisson distribution reduces to the product of two independent Poisson distributions (referred to as double Poisson distribution “DP”²).

By means of conditional probability theory, the joint density (1) can be written as the product of a marginal and a conditional distribution. Hence:

In general, every decomposition of $f(y_1\text{,}{y}_2)$ will lead to different marginal and conditional distributions.

Returning to Eq. (2), $Y_1$ is a Poisson distribution with parameter ${\mu }_1={\lambda }_1+{\lambda }_3$ and the conditional distribution of $Y_2$ given $Y_1$ is expressed as

The expression in (4) is a convolution of a Poisson variable with parameter ${\lambda }_2$ and a binomial with parameters $(y_1\text{,}{p}_1)$ . The conditional mean and variance are given by: $$\begin{array}{cc}\hfill E(Y_2\left|Y_1\right.)& ={\lambda }_2+p_1{y}_1\text{,}\hfill \\ \hfill \mathit{Var}(Y_2\left|Y_1\right.)& ={\lambda }_2+p_1(1-p_1)y_1\text{.}\hfill \end{array}$$

Hence, the joint p.m.f of conditional model 1 “CM1”³ given in (2):

For this model $\mathit{Cov}(y_1\text{,}{y}_2)=p_1{\mathrm{\mu }}_1$ . Hence, the covariance is zero if and only if $p_1=0$ (note that ${\mathrm{\mu }}_1=0$ correspond to the univariate case).

In the same way, Eq. (3) can be written as:

Also, for this model $\mathit{Cov}(y_1\text{,}{y}_2)=p_2{\mathrm{\mu }}_2$ . Hence, the covariance is zero if and only if $p_2=0$ .

Note that for CM1 and CM2, $p_1$ and $p_2$ play the same role as ${\lambda }_3$ in JBP.

The bivariate Poisson model was applied by King (1989) to the annual number of presidential vetoes of social welfare bills and defense bills, by Jung and Winkelmann (1993) to the number of voluntary and involuntary changes, and by Ozuna and Gomez (1994) to the number of trips to different recreational sites. They model the marginal expectation of $Y_1$ and $Y_2$ as a log linear function of explanatory variables.

The disadvantage of this particular model is that it does not allow for over/under dispersion (the marginal distributions are Poisson) or negative correlation, and thus lacks generality. For the case of over-dispersed count data, the mixed Poisson models are potentially useful (Bermúdez and Karlis, 2012; Ghitany et al., 2012; Gurmu and Elder, 2000; Munkin and Trivedi, 1999). There are also some other models that allow for negative correlation (Berkhout and Plug, 2004; Chib and Winkelmann, 2001; Gurmu and Elder, 2008; Karlis and Meligkotsidou, 2007; Van Ophem, 1999).

It is clear that there are one-to-one transformations between the parameters of the three representations of the bivariate Poisson distribution. Hence, according to the maximum likelihood invariance principle, we obtain identical estimates of the parameters no matter which representation has been used. It will be of interest to see whether such an invariance property holds when we have explanatory variables.

In this paper, we compare fitting a bivariate Poisson regression model using both the joint and conditional arguments described above. We also consider the inflated versions of these models to allow for over-dispersion and discuss inferences related to the parameters involved in the models and the execution times, using SAS and R software.

The outline of the paper is as follows. In Section 2, we present the bivariate Poisson regression model and discuss the estimation procedure using both the joint and conditional procedures. In Section 3, we discuss the zero-inflation version of the joint bivariate Poisson regression model and the proposed conditional models. In Section 4, a simulation study is conducted to compare the three models. In Section 5, applications on two data sets from Australia and Saudi Arabia are illustrated using the considered bivariate regression models. Finally, some concluding remarks are provided in Section 6.

2 Bivariate Poisson regression model

Here, we assume that the parameters of the models depend on explanatory variables. In the joint bivariate Poisson regression model “JBPM”⁴, ${\lambda }_k>0$ with k = 1, 2 and 3 can be related to various explanatory variables using the classical exponential link functions. Therefore, the joint bivariate Poisson regression model can take the following form:

In the case of the explanatory variables, two aspects should be stressed. First, different covariates can be used to model each parameter (or same covariates) and second, covariates can be introduced to model ${\lambda }_3$ in order to learn more about the influence of the covariates on each pair of variables (or, to facilitate interpretation, no covariates are used to model ${\lambda }_3$ ).

Gourieroux et al. (1984) derived pseudo maximum likelihood estimation methods. Jung and Winkelmann (1993) and Kocherlakota and Kocherlakota (2001) considered the joint bivariate Poisson regression model using a Newton–Raphson procedure. Karlis and Ntzoufras (2005) constructed an EM algorithm to remedy convergence problems encountered with the Newton Raphson procedure. Ho and Singer (2001) proposed a generalized least squares method for maximizing the log likelihood. Karlis and Tsiamyrtzis (2008) implemented a Bayesian inference that does not rely on the MCMC scheme. Tsionas (2001), Ma and Kockelman (2006) and Choe et al. (2012) considered a Bayesian approach based on the MCMC technique to execute some computations.

Here we use the maximum likelihood method for estimating the parameters. Consider independent observations $(y_{1i}\text{,}{y}_{2i})$ with the i-th vector having the join bivariate Poisson distribution $$f(y_{1i}\text{,}{y}_{2i})=e^{-({\lambda }_{1i}+{\lambda }_{2i}+{\lambda }_{3i})}\phi (y_{1i}\text{,}{y}_{2i})$$ where $$\phi (y_{1i}\text{,}{y}_{2i})={\displaystyle \sum _{r=0}^{\mathit{min}(y_{1i}\text{,}{y}_{2i})}}\frac{{\lambda }_{3i}^r{\lambda }_{1i}^{(y_{1i}-r)}{\lambda }_{2i}^{(y_{2i}-r)}}{r!(y_{1i}-r)!(y_{2i}-r)!}\text{.}$$

Then, the corresponding score functions $U_{\mathit{kj}}=\frac{\partial \log L}{\partial {\alpha }_{\mathit{kj}}}\text{,}k=1\text{,}2\text{,}3\text{;}j=0\text{,}1\text{,}\dots \text{,}l$ are $$\begin{array}{cc}\hfill U_{1j}& ={\displaystyle \sum _{i=1}^n}{\lambda }_{1i}{w}_{\mathit{ji}}\left[\frac{\phi (y_{1i}-1\text{,}{y}_{2i})}{\phi (y_{1i}\text{,}{y}_{2i})}-1\right]\text{,}\hfill \\ \hfill U_{2j}& ={\displaystyle \sum _{i=1}^n}{\lambda }_{2i}{w}_{\mathit{ji}}\left[\frac{\phi (y_{1i}\text{,}{y}_{2i}-1)}{\phi (y_{1i}\text{,}{y}_{2i})}-1\right]\text{,}\hfill \\ \hfill U_{3j}& ={\displaystyle \sum _{i=1}^n}{\lambda }_{3i}{w}_{\mathit{ji}}\left[\frac{\phi (y_{1i}-1\text{,}{y}_{2i}-1)}{\phi (y_{1i}\text{,}{y}_{2i})}-1\right]\text{.}\hfill \end{array}$$

For the conditional Poisson regression models, we will only consider the first conditional model (CM1), as the argument for the second conditional model (CM2) is the same. The likelihood of CM1 can be written as the product of two likelihoods; the first one is that of the conditional model and the second is that of the marginal model. Hence, assuming that the regression parameters are different, one can maximize each likelihood separately with respect to the corresponding parameters. For the CM1 we denote the parameters by ${\mathit{\beta }}_k\text{,}k=1\text{,}2\text{,}3$ .

For the marginal Poisson regression model, let ${\mu }_{1i}$ be the mean of $Y_{1i}\text{,}i=1\text{,}2\text{,}\dots \text{,}n$ . Then, using the standard Poisson regression model, we set:

For the conditional regression model, the conditional density (4) of $\left.Y_{2i}\right|Y_{1i}\text{,}i=1\text{,}2\text{,}\dots \text{,}n$ involves the parameters $p_{1i}$ and ${\lambda }_{2i}$ . We use the following link functions:

The score functions $V_{\mathit{kj}}=\frac{\partial \log L}{\partial {\beta }_{\mathit{kj}}}\text{,}k=1\text{,}2\text{,}3\text{;}j=0\text{,}1\text{,}\dots \text{,}l$ for each regression parameters are summarized as follows:

By equating the score functions in (11) to zero, one can obtain the parameter estimates directly. We note that the likelihood equations are non-linear and do not have a closed-form solution, so they can be solved numerically.

The following theorem gives the relations between the score functions of JBPM and CM1.

Theorem 1

Note that when we set $U_{\mathit{kj}}=0$ for all $k=1\text{,}2\text{,}3$ and $j=0\text{,}1\text{,}\dots \text{,}l$ we get $V_{\mathit{kj}}=0$ , for all $k=1\text{,}2\text{,}3$ and $j=0\text{,}1\text{,}\dots \text{,}l$ and vice versa. This shows that the likelihood invariance property holds when we have explanatory variables.

3 A bivariate zero-inflated model

The zero-inflated model is used when a count data set shows a large proportion of zeros. A bivariate zero-inflated model can be constructed by increasing the probability of the event $(y_1=0\text{,}{y}_2=0)$ and decreasing the other joint probabilities.

The zero-inflated joint bivariate Poisson model “ZIJBPM”⁵ is specified by the probability function:

This model was proposed by Karlis and Ntzoufras (2003, 2005) in order to allow for over-dispersion of the corresponding marginal distributions. It was applied by Wang et al. (2003) to analyze two types of occupational injuries and by Bermúdez (2009) in the automobile insurance context for a bivariate case.

One can easily see that the marginal distributions are no longer simple Poisson distributions, but are now zero-inflated versions. Note that one may define more complicated models by assuming other kinds of inflations (Karlis and Ntzoufras, 2005). Moreover, one may add covariates to $\pi$ , implying that inflation depends on external factors.

For estimation of the parameters of ZIJBPM, we use the maximum likelihood procedure as follows. Let $(y_{1i}\text{,}{y}_{2i})\text{,}i=1\text{,}2\text{,}\dots \text{,}n$ be a random sample observed from ZIJBPM. The likelihood function for the observed random sample is given by: $$L({\pi }_i\text{,}{\lambda }_{1i}\text{,}{\lambda }_{2i}\text{,}{\lambda }_{3i}\text{;}{y}_{1i}\text{,}{y}_{2i})={\displaystyle \prod _{i=1}^n}\{{({\pi }_i+(1-{\pi }_i)f_{\mathit{JBP}}(0\text{,}0\text{;}{\lambda }_{1i}\text{,}{\lambda }_{2i}\text{,}{\lambda }_{3i}))}^{1-a_i}{((1-{\pi }_i)f_{\mathit{JBP}}(y_{1i}\text{,}{y}_{2i}\text{;}{\lambda }_{1i}\text{,}{\lambda }_{2i}\text{,}{\lambda }_{3i}))}^{a_i}\}\text{,}$$ where $a_i=1$ if $(y_{1i}\text{,}{y}_{2i})\ne (0\text{,}0)$ and $a_i=0$ otherwise. In order to fit this model, we used the link functions in (7) for ${\lambda }_{\mathit{ki}}$ , and for modeling the mixing proportion ${\pi }_i$ , we used a logit link function: $$\text{logit}({\pi }_i)={\gamma }_0+{\displaystyle \sum _{j=1}^l}{\gamma }_j{w}_{\mathit{ji}}\text{.}$$

The maximum likelihood estimates (MLEs) of the parameters can be obtained by solving equations $\frac{\partial \log L}{\partial {\alpha }_{1j}}=0$ , $\frac{\partial \log L}{\partial {\alpha }_{2j}}=0$ , $\frac{\partial \log L}{\partial {\alpha }_{3j}}=0$ and $\frac{\partial \log L}{\partial {\gamma }_j}=0$ for $j=0\text{,}1\text{,}\dots \text{,}l$ , simultaneously.

The zero-inflated conditional model “ZICM”⁶ is given by:

For parameter estimations of ZICM1, let $(y_{1i}\text{,}{y}_{2i})\text{,}i=1\text{,}2\text{,}\dots \text{,}n$ be a random sample observed from ZICM1. The likelihood function for the observed random sample is given by $$L({\pi }_i\text{,}{\lambda }_{1i}\text{,}{\lambda }_{2i}\text{,}{\lambda }_{3i}\text{;}{y}_{1i}\text{,}{y}_{2i})={\displaystyle \prod _{i=1}^n}\{{({\pi }_i+(1-{\pi }_i)f_{\text{CM1}}(0\text{,}0\text{;}{p}_{1i}\text{,}{\lambda }_{2i}\text{,}{\mathrm{\mu }}_{1i}))}^{1-a_i}{((1-{\pi }_i)f_{\text{CM1}}(y_{1i}\text{,}{y}_{2i}\text{;}{p}_{1i}\text{,}{\lambda }_{2i}\text{,}{\mathrm{\mu }}_{1i}))}^{a_i}\}\text{,}$$ where $a_i=1$ if $(y_{1i}\text{,}{y}_{2i})\ne (0\text{,}0)$ and $a_i=0$ otherwise. In order to fit this model, we used the link functions, as in (8)–(10), for ${\mu }_{1i}\text{,}{\lambda }_{2i}$ and $p_{1i}$ , respectively, and to model the mixing proportion ${\pi }_i$ , we used a logit link function: $$\text{logit}({\pi }_i)={\gamma }_0+{\displaystyle \sum _{j=1}^l}{\gamma }_j{w}_{\mathit{ji}}\text{.}$$

The MLEs of the parameters can be obtained by solving equations $\frac{\partial \log L}{\partial {\beta }_{1j}}=0$ , $\frac{\partial \log L}{\partial {\beta }_{2j}}=0$ , $\frac{\partial \mathit{LogL}}{\partial {\beta }_{3j}}=0$ and $\frac{\partial \mathit{LogL}}{\partial {\gamma }_j}=0$ for $j=0\text{,}1\text{,}\dots \text{,}l$ , simultaneously.

4 Experiment results

We conducted a simulation study to examine the performance of the proposed conditional model compared with that of the joint bivariate Poisson model and the double Poisson model. We used three variables, namely gender, age and income from Health Care Australian data (Cameron et al., 1988) as covariates. For each of the following four cases we generated 1000 samples each of size 5190.

Case 1. Simulation from joint bivariate Poisson model.

We have simulated the data points $(y_{1i}\text{,}{y}_{2i})$ from the joint bivariate Poisson regression model with ${\lambda }_{1i}\text{,}{\lambda }_{2i}\text{,}{\lambda }_{3i}$ given by: $$\begin{array}{cc}\hfill {\lambda }_{1i}& =\exp ({\alpha }_{10}+{\alpha }_{11}{w}_{1i}+{\alpha }_{12}{w}_{2i}+{\alpha }_{13}{w}_{3i})\text{,}\hfill \\ \hfill {\lambda }_{2i}& =\exp ({\alpha }_{20}+{\alpha }_{21}{w}_{1i}+{\alpha }_{22}{w}_{2i}+{\alpha }_{23}{w}_{3i})\text{,}\hfill \\ \hfill {\lambda }_{3i}& =\exp ({\alpha }_{30}+{\alpha }_{31}{w}_{1i}+{\alpha }_{32}{w}_{2i}+{\alpha }_{33}{w}_{3i})\text{.}\hfill \end{array}$$

We have used the following parameter values: $${\alpha }_{10}=-1.9983\text{,}{\alpha }_{11}=0.2268\text{,}{\alpha }_{12}=0.6902\text{,}{\alpha }_{13}=-0.2004\text{,}{\alpha }_{20}=-2.1718\text{,}{\alpha }_{21}=0.6436\text{,}{\alpha }_{22}=3.1092\text{,}{\alpha }_{23}=-0.07716\text{,}{\alpha }_{30}=-3.0328\text{,}{\alpha }_{31}=0.1966\text{,}{\alpha }_{32}=2.1272\text{,}{\alpha }_{33}=-0.3984\text{.}$$

Case 2. Simulation from conditional model.

We generated the data points $(y_{1i}\text{,}{y}_{2i})$ from the conditional regression model 1, assuming the following: $$\begin{array}{cc}\hfill {\mu }_{1i}& =\exp ({\beta }_{10}+{\beta }_{11}{w}_{1i}+{\beta }_{12}{w}_{2i}+{\beta }_{13}{w}_{3i})\text{.}\hfill \\ \hfill {\lambda }_{2i}& =\exp ({\beta }_{20}+{\beta }_{21}{w}_{1i}+{\beta }_{22}{w}_{2i}+{\beta }_{23}{w}_{3i})\text{,}\hfill \\ \hfill \text{logit}(p_{1i})& ={\beta }_{30}+{\beta }_{31}{w}_{1i}+{\beta }_{32}{w}_{2i}+{\beta }_{33}{w}_{3i}\text{.}\hfill \end{array}$$

Then, $$\begin{array}{cc}\hfill {\lambda }_{3i}& =p_{1i}\ast {\mu }_{1i}\text{,}\hfill \\ \hfill {\lambda }_{1i}& ={\mu }_{1i}-{\lambda }_{3i}\text{.}\hfill \end{array}$$

We have used the following parameter values: $${\beta }_{10}=-1.71473\text{,}{\beta }_{11}=0.21565\text{,}{\beta }_{12}=1.23798\text{,}{\beta }_{13}=-0.27726\text{,}{\beta }_{20}=-2.1658\text{,}{\beta }_{21}=0.6316\text{,}{\beta }_{22}=3.1162\text{,}{\beta }_{23}=-0.08287\text{,}{\beta }_{30}=-1.2005\text{,}{\beta }_{31}=0.1642\text{,}{\beta }_{32}=1.4778\text{,}{\beta }_{33}=-0.08393\text{.}$$

Case 3. Simulation from zero-inflated joint bivariate Poisson model.

We have simulated the data points $(y_{1i}\text{,}{y}_{2i})$ from the zero-inflated joint bivariate Poisson regression model with ${\lambda }_{1i}\text{,}{\lambda }_{2i}\text{,}{\lambda }_{3i}$ and ${\pi }_i$ given by: $$\begin{array}{cc}\hfill {\lambda }_{1i}& =\exp ({\alpha }_{10}+{\alpha }_{11}{w}_{1i}+{\alpha }_{12}{w}_{2i}+{\alpha }_{13}{w}_{3i})\text{,}\hfill \\ \hfill {\lambda }_{2i}& =\exp ({\alpha }_{20}+{\alpha }_{21}{w}_{1i}+{\alpha }_{22}{w}_{2i}+{\alpha }_{23}{w}_{3i})\text{,}\hfill \\ \hfill {\lambda }_{3i}& =\exp ({\alpha }_{30}+{\alpha }_{31}{w}_{1i}+{\alpha }_{32}{w}_{2i}+{\alpha }_{33}{w}_{3i})\text{,}\hfill \\ \hfill \text{logit}({\pi }_i)& =\gamma \text{.}\hfill \end{array}$$

We have used the following parameter values: $${\alpha }_{10}=-1.4601\text{,}{\alpha }_{11}=0.2223\text{,}{\alpha }_{12}=0.7534\text{,}{\alpha }_{13}=-0.1644\text{,}{\alpha }_{20}=-1.6012\text{,}{\alpha }_{21}=0.5709\text{,}{\alpha }_{22}=2.8359\text{,}{\alpha }_{23}=-0.03434\text{,}{\alpha }_{30}=-2.7627\text{,}{\alpha }_{31}=-0.2254\text{,}{\alpha }_{32}=1.9464\text{,}{\alpha }_{33}=-0.5794\text{,}\gamma =-0.7853\text{.}$$

Case 4. Simulation from zero-inflated conditional model.

We generated the data points $(y_{1i}\text{,}{y}_{2i})$ from the zero-inflated conditional regression model 1, assuming the following: $$\begin{array}{cc}\hfill {\mu }_{1i}& =\exp ({\beta }_{10}+{\beta }_{11}{w}_{1i}+{\beta }_{12}{w}_{2i}+{\beta }_{13}{w}_{3i})\text{.}\hfill \\ \hfill {\lambda }_{2i}& =\exp ({\beta }_{20}+{\beta }_{21}{w}_{1i}+{\beta }_{22}{w}_{2i}+{\beta }_{23}{w}_{3i})\text{,}\hfill \\ \hfill \text{logit}(p_{1i})& ={\beta }_{30}+{\beta }_{31}{w}_{1i}+{\beta }_{32}{w}_{2i}+{\beta }_{33}{w}_{3i}\text{,logit}({\pi }_i)=\gamma \text{.}\hfill \end{array}$$

Then, $${\lambda }_{3i}=p_{1i}*{\mu }_{1i}\text{,}$$ $${\lambda }_{1i}={\mu }_{1i}-{\lambda }_{3i}\text{.}$$

We have used the following parameter values: $${\beta }_{10}=-1.2134\text{,}{\beta }_{11}=0.1149\text{,}{\beta }_{12}=1.034\text{,}{\beta }_{13}=-0.257\text{,}{\beta }_{20}=-1.6096\text{,}{\beta }_{21}=0.5437\text{,}{\beta }_{22}=2.8886\text{,}{\beta }_{23}=-0.04059\text{,}{\beta }_{30}=-1.2839\text{,}{\beta }_{31}=-0.1031\text{,}{\beta }_{32}=0.6529\text{,}{\beta }_{33}=-0.2774\text{,}\gamma =-0.789\text{.}$$

Table 1 illustrates the means and standard deviations of the AIC for the fitted models under the four cases. Bias and mean square error (MSE) of the estimated parameters are presented in Table 2. Table 1 shows that the average values of the AIC for the joint bivariate Poisson model are almost identical with those obtained from conditional model 1 irrespective of the mechanism of data generations. On the other hand, when the data are generated from zero-inflated models, the fitted models ZIDPM, ZIJBP and ZICM1 perform better than their non-inflated counter parts, as expected.

To perform official test for the hypothesis that there is no significant difference of the AIC’s, we preformed Vuong test. Vuong (1989) set the information criterion in a testing framework in which the null hypothesis is that the two competing models are equally close to the true model. Under the null hypothesis that the models are indistinguishable, the test statistic is asymptotically distributed standard normal. We applied the Vuong test to investigate if there is statistically significant difference between the joint bivariate Poisson model and the conditional model 1 in terms of AICs. We also performed the test to compare the zero inflated models. In all cases the Vuong test confirmed that there is no significant difference between the two models. The smallest p-value obtained was 0.529. Both models always perform better than the double Poisson model.

In general, the simulation study confirmed the expected similarity of the results obtained under joint and conditional bivariate Poisson models. This concludes the invariance property holds when we have explanatory variables as expected from the results of Theorem 1.

5 Applications

6 Conclusion

In this paper, we used conditional argument to introduce a two-stage procedure for estimating the bivariate Poisson regression model. Our study showed that this method has the same performance as that of using the joint density. The conditional method is simple as it partitions the joint likelihood into two univariate likelihoods. We also applied the two methods to two sets of real data. The conclusion is that the two ways of modeling have the same AIC but the execution times for the two-stage models are shorter.