Pseudo-Poisson Distributions with Concomitant Variables

Abstract

It has been argued in Arnold and Manjunath (2021) that the bivariate pseudo-Poisson distribution will be the model of choice for bivariate data with one equidispersed marginal and the other marginal over-dispersed. This is due to its simple structure, straightforward parameter estimation and fast computation. In the current note, we introduce the effects of concomitant variables on the bivariate pseudo-Poisson parameters and explore the distributional and inferential aspects of the augmented models. We also include a small simulation study and an example of application to real-life data.

1. Introduction

The classical “one-dimensional” Poisson distribution has historically been found to be useful in modeling a wide variety of “integer-valued” phenomena, such as the number of accidents and associated fatalities, disease advances, rate of rare event occurrences and so on. With regard to the Poisson model with concomitants, i.e., Poisson regression or count regression, its best known applications are in $(i)$ modeling counts of bacteria exposed to various environmental conditions and dilutions; $(ii)$ modeling counts of infant mortality among groups with different demographics. All these examples are typically modeled under the assumption of equi-dispersion. However, count data also exhibits over and under dispersion. In this context, the one-dimensional Conway–Maxwell–Poisson model or its regression version fills the bill precisely by allowing us to model over, equi- and under-dispersion data.

In general, bivariate count data, along with having marginal over-, under- or equi-dispersion, will also exhibit a variety of dependence structures. In particular, for linear dependence, the possible relations are positive or negative correlation. The classical bivariate Poisson model is appropriate for data having equi-dispersed marginals with positive correlation. Here again, the bivariate Conway–Maxwell–Poisson is more flexible in that it can adapt to both under and over dispersed data, see Sellers et al. [1]. Concerning bivariate Poisson regression models, the first version involving explanatory variables acting on the marginal means was introduced in Kocherlakota and Kocherlakota [2] based on the classical bivariate Poisson model. In addition, the derivation of Wald, score and likelihood ratio test statistics for testing a single coefficient parameter vector are discussed in Riggs et al. [3]. Zamani et al. [4] proposed a bivariate Poisson model which can be fitted to both positive and negatively correlated data. Recently, Chowdhury et al. [5] considered the Poisson–Poisson regression model (which is the particular case of the bivariate pseudo-Poisson model) to analyze the impact of covariates on the daily new cases and fatalities associated with the COVID-19 pandemic. Finally, we refer to Karlis and Ntzoufras [6] and the R package $bivpois$ for maximum likelihood estimation, using an Expectation-Maximization (EM) algorithm, for diagonally inflated bivariate Poisson regression models.

In recent work, Arnold and Manjunath [7] recommended the bivariate pseudo-Poisson model to fit data which have one marginal and other conditional of the Poisson form. Due to its straightforward structure with no restrictions on the conditional mean function, it allows us to include a variety of dependence structures, including positive and negative correlation. In the following, we introduce explanatory variables acting on the pseudo-Poisson parameters. Thanks to the simple structure, the concomitant effects can be introduced into each of the parameters to generate a family of models with a variety of dependence structures. We refer to Arnold et al. [8] and Veeranna et al. [9] on Bayesian and goodness-of-fit tests for the bivariate pseudo-Poisson model, respectively, which can also be adapted to accommodate the presence of concomitant variables. We refer to Arnold et al. [10] and Filus et al. [11] for further reading on conditional specified models and the triangular transformations, respectively. Finally, we refer Ghosh et al. [12] on the recent results on bivariate count model which has both conditionals with a Poisson structure.

We next review the concept of multivariate pseudo-Poisson distributions, as discussed in Arnold and Manjunath [7].

Definition 1.

A k-dimensional random variable $\underline{X}=(X_1,X_2,\cdots ,X_k)$ is said to have a k-dimensional pseudo-Poisson distribution if there exists a positive constant ${\lambda }_1$ such that

$$X_1\sim \mathcal{P}({\lambda }_1)$$

and $k-1$ functions $\{{\lambda }_{\ell }:\ell =2,3,\cdots ,k\}$ where, for each ℓ, ${\lambda }_{\ell }:{\{0,1,2,\cdots \}}^{(\ell -1)}\to (0,\infty )$ such that

$$X_{\ell }|{\underline{X}}_{(\ell -1)}={\underline{x}}_{(\ell -1)}\sim \mathcal{P}({\lambda }_{\ell }({\underline{x}}_{(\ell -1)})),$$

where${\underline{X}}_{(\ell -1)}={(X_1,\cdots ,X_{l-1})}^{\top }$. Note that there are no constraints on the forms of the functions ${\lambda }_{\ell },\ell =2,3,\cdots ,k$ that appear in the definition, save for measurability. In applications, it would typically be the case that the ${\lambda }_{\ell }$’s would be chosen to be relatively simple functions depending on a limited number of parameters.

Definition 2.

A random pair of variables $(X_1,X_2)$ is said to have a bivariate pseudo-Poisson distribution if there exists a positive constant ${\lambda }_1$ such that

$$X_1\sim \mathcal{P}({\lambda }_1)$$

and a function ${\lambda }_2:\{0,1,2,\cdots \}\to (0,\infty )$ such that, for every non-negative integer $x_1$,

$$X_2|X_1=x_1\sim \mathcal{P}({\lambda }_2(x_1)).$$

The fact that there are no constraints on the ${\lambda }_2(x_1)$ allows us to adapt to a variety of dependence structures including positive or negative correlation.

Example 1.

A judicious choice of a parametric family for ${\lambda }_2(x_1)$ will admit positive and negative correlation between $X_1$ and $X_2$. For example, if we consider

$${\lambda }_2(x_1;\gamma ,\delta )=1+(2\gamma -1)(1-e^{-\delta x_1}).$$

(1)

For $\delta >0$, the above function will be increasing if $\gamma >1/2$, decreasing if $\gamma <1/2$ and constant if $\gamma =1/2$. Consequently, $X_1$ and $X_2$ will have a positive correlation if $\gamma >1/2$, negative correlation if $\gamma <1/2$ and will be uncorrelated if $\gamma =1/2.$ A more general model with the same properties can be obtained by replacing $1-e^{-\delta x_1}$ by $F(x_1;\underline{\theta })$, a parameterized family of distribution functions with support $(0,\infty )$.

2. Incorporating Concomitant Variables

In many (perhaps, most) applications, in addition to the observed values of the $(X_{1,i},X_{2,i})$’s pairs, there will be available values of arrays of concomitant variables which are expected to influence the stochastic behavior of the observed data points. A straightforward manner in which to incorporate vectors of concomitant variables ${\underline{u}}_i={(u_{1i},\cdots ,u_{di})}^{\top }$ into the model is as follows:

$$X_1\sim \mathcal{P}\left({\lambda }_{1}exp\left({\underline{\alpha }}^{\top }\underline{u}\right)\right)$$

(2)

and

$$X_2|X_1=x_1\sim \mathcal{P}\left({\lambda }_{2}exp\left({\underline{\beta }}^{\top }\underline{u}\right)+{\lambda }_{3}exp\left({\underline{\gamma }}^{\top }\underline{u}\right)x_1\right)$$

(3)

where ${\lambda }_1>0$, ${\lambda }_2\ge 0$, ${\lambda }_3>0$, $\underline{\alpha }={({\alpha }_1,\cdots ,{\alpha }_d)}^{\top }$, $\underline{\beta }={({\beta }_1,\cdots ,{\beta }_d)}^{\top }$ and $\underline{\gamma }={({\gamma }_1,\cdots ,{\gamma }_d)}^{\top }$ are d-dimensional unknown parameters.

There are certainly many other manners in which one can model the influence of concomitant variables. If there is scientific justification for alternative models that do not introduce the concomitants via log-linear adjustments of the form specified in (2) and (3), then one should certainly utilize the scientifically appropriate link functions.

Just as in classical multiple regression, it is worthwhile to determine whether a simple linear dependence assumption for the effect of concomitants will be adequate to fit the data. In the remainder of this paper, we will focus on the simple model (2) and (3).

3. Moments

In the following, we derive some population moments for the model specified in (2) and (3).

$$\begin{array}{ccc}\hfill E(X_1)=Var(X_1)& =& \hfill {\lambda }_{1}exp\left({\underline{\alpha }}^{\top }\underline{u}\right)\hfill \\ \hfill E(X_2)& =& \hfill {\lambda }_{2}exp\left({\underline{\beta }}^{\top }\underline{u}\right)+{\lambda }_1{\lambda }_{3}exp\left({(\underline{\gamma }+\underline{\alpha })}^{\top }\underline{u}\right)\hfill \\ \hfill V(X_2)& =& \hfill {\lambda }_{2}exp\left({\underline{\beta }}^{\top }\underline{u}\right)+{\lambda }_1{\lambda }_{3}exp\left({(\underline{\gamma }+\underline{\alpha })}^{\top }\underline{u}\right)\hfill \\ & & \hfill +{\lambda }_1{\lambda }_3^{2}exp\left({(2\underline{\gamma }+\underline{\alpha })}^{\top }\underline{u}\right)\hfill \\ \hfill Cov(X_1,X_2)& =& \hfill {\lambda }_1{\lambda }_{3}exp\left({(\underline{\gamma }+\underline{\alpha })}^{\top }\underline{u}\right).\hfill \end{array}$$

The marginal dispersion indices are

$$\begin{array}{ccc}\hfill DI(X_1)& =& \hfill \frac{Var(X_1)}{E(X_1)}=1.\hfill \\ \hfill DI(X_2)& =& \hfill \frac{Var(X_2)}{E(X_2)}=1+\frac{{\lambda }_1{\lambda }_3^{2}exp\left({(2\underline{\gamma }+\underline{\alpha })}^{\top }\underline{u}\right)}{{\lambda }_{2}exp\left({\underline{\beta }}^{\top }\underline{u}\right)+{\lambda }_1{\lambda }_{3}exp\left({(\underline{\gamma }+\underline{\alpha })}^{\top }\underline{u}\right)}.\hfill \end{array}$$

For ${\lambda }_2=0$

$$\begin{array}{c}\hfill DI(X_2)=1+{\lambda }_{3}exp\left({\underline{\gamma }}^{\top }\underline{u}\right).\end{array}$$

Define

$$\begin{array}{ccc}\hfill E(\underline{X})& =& \hfill {(E(X_1),E(X_2))}^{\top }\hfill \\ & =& \hfill {\left({\lambda }_{1}exp\left({\underline{\alpha }}^{\top }\underline{u}\right),{\lambda }_{2}exp\left({\underline{\beta }}^{\top }\underline{u}\right)+{\lambda }_1{\lambda }_{3}exp\left({(\underline{\gamma }+\underline{\alpha })}^{\top }\underline{u}\right)\right)}^{\top }\hfill \end{array}$$

and

$$\begin{array}{c}\hfill Cov(\underline{X})=\left[\begin{array}{cc}{\lambda }_{1}exp\left({\underline{\alpha }}^{\top }\underline{u}\right)& {\lambda }_1{\lambda }_{3}exp\left({(\underline{\gamma }+\underline{\alpha })}^{\top }\underline{u}\right)\\ {\lambda }_1{\lambda }_{3}exp\left({(\underline{\gamma }+\underline{\alpha })}^{\top }\underline{u}\right)& {\lambda }_{2}exp\left({\underline{\beta }}^{\top }\underline{u}\right)+{\lambda }_1{\lambda }_{3}exp\left({(\underline{\gamma }+\underline{\alpha })}^{\top }\underline{u}\right)+\\ & +{\lambda }_1{\lambda }_3^{2}exp\left({(2\underline{\gamma }+\underline{\alpha })}^{\top }\underline{u}\right)\end{array}\right]\end{array}$$

By using the definition given in the paper by Kokonendji and Puig [13] page 183, the bivariate Fisher index of dispersion is given by

$$\begin{array}{c}\hfill GDI(\underline{X})=\frac{\sqrt{E{(\underline{X})}^{\top }}Cov(\underline{X})\sqrt{E(\underline{X})}}{E{(\underline{X})}^{\top }E(\underline{X})}\end{array}$$

which is a case of over-disperson, cf. Arnold and Manjunath [7] page 2311 for the dispersion index proof for the bivariate pseudo-Poisson distribution.

4. Statistical Inference

In this section, we obtain maximum likelihood estimators (m.l.e.) of parameters ${\lambda }_1$, ${\lambda }_2$, ${\lambda }_3$, $\underline{\alpha }$, $\underline{\beta }$ and $\underline{\gamma }$. In addition, we construct the likelihood ratio test for the possible parallelism, coincidence and significance of each of the regression coefficients.

4.1. Estimation

Let ${(X_{1i},X_{2i})}^{\top }$, $i=1,2,\cdots ,n$ be a bivariate count sample from the pseudo-Poisson distribution (in Section 2) and let ${\underline{u}}_1,\cdots ,{\underline{u}}_n$ be d-dimensional known covariates. Then, the log-likelihood function is

$$\begin{array}{ccc}\hfill logL& =& \hfill -{\lambda }_1\sum _{i=1}^{n}exp\left({\underline{\alpha }}^{\top }{\underline{u}}_i\right)+\sum _{i=1}^n{x}_{1i}log\left({\lambda }_{1}exp\left({\underline{\alpha }}^{\top }{\underline{u}}_i\right)\right)\hfill \\ & & \hfill -\sum _{i=1}^n\left({\lambda }_{2}exp\left({\underline{\beta }}^{\top }{\underline{u}}_i\right)+{\lambda }_{3}exp\left({\underline{\gamma }}^{\top }{\underline{u}}_i\right)x_{1i}\right)\hfill \\ & & \hfill +\sum _{i=1}^n{x}_{2i}log\left({\lambda }_{2}exp\left({\underline{\beta }}^{\top }{\underline{u}}_i\right)+{\lambda }_{3}exp\left({\underline{\gamma }}^{\top }{\underline{u}}_i\right)x_{1i}\right)\hfill \\ & & \hfill -\sum _{i=1}^{n}log(x_{1i}!x_{2i}!).\hfill \end{array}$$

(4)

Partial differentiation with respect to each parameters ${\lambda }_1$, ${\lambda }_2$ and ${\lambda }_3$ and equating to zero gives

$$\begin{array}{c}\hfill -\sum _{i=1}^{n}exp\left({\underline{\alpha }}^{\top }{\underline{u}}_i\right)+\sum _{i=1}^n{x}_{1i}\frac{exp\left({\underline{\alpha }}^{\top }{\underline{u}}_i\right)}{{\lambda }_{1}exp\left({\underline{\alpha }}^{\top }{\underline{u}}_i\right)}=0\end{array}$$

(5)

$$\begin{array}{c}\hfill -\sum _{i=1}^{n}exp\left({\underline{\beta }}^{\top }{\underline{u}}_i\right)+\sum _{i=1}^n{x}_{2i}\frac{exp\left({\underline{\beta }}^{\top }{\underline{u}}_i\right)}{{\lambda }_{2}exp\left({\underline{\beta }}^{\top }{\underline{u}}_i\right)+{\lambda }_{3}exp\left({\underline{\gamma }}^{\top }{\underline{u}}_i\right)x_{1i}}=0\end{array}$$

(6)

$$\begin{array}{c}\hfill -\sum _{i=1}^n{x}_{1i}exp\left({\underline{\gamma }}^{\top }{\underline{u}}_i\right)+\sum _{i=1}^n{x}_{1i}{x}_{2i}\frac{exp\left({\underline{\gamma }}^{\top }{\underline{u}}_i\right)}{{\lambda }_{2}exp\left({\underline{\beta }}^{\top }{\underline{u}}_i\right)+{\lambda }_{3}exp\left({\underline{\gamma }}^{\top }{\underline{u}}_i\right)x_{1i}}=0.\end{array}$$

(7)

Now, taking partial derivatives of $logL$ with respect to ${\alpha }_j$, ${\beta }_j$ and ${\gamma }_j$ for $j\in \{1,\cdots ,d\}$ and equating to zero yields

$$\begin{array}{c}\hfill -{\lambda }_1\sum _{i=1}^n{u}_{ji}exp\left({\underline{\alpha }}^{\top }{\underline{u}}_i\right)+\sum _{i=1}^n{x}_{1i}{u}_{ji}=0\end{array}$$

(8)

$$\begin{array}{c}\hfill -\sum _{i=1}^n{u}_{ji}exp\left({\underline{\beta }}^{\top }{\underline{u}}_i\right)+\sum _{i=1}^n{x}_{2i}{u}_{ji}\frac{exp\left({\underline{\beta }}^{\top }{\underline{u}}_i\right)}{{\lambda }_{2}exp\left({\underline{\beta }}^{\top }{\underline{u}}_i\right)+{\lambda }_{3}exp\left({\underline{\gamma }}^{\top }{\underline{u}}_i\right)x_{1i}}=0\end{array}$$

(9)

$$\begin{array}{c}\hfill -\sum _{i=1}^n{x}_{1i}{u}_{ji}exp\left({\underline{\gamma }}^{\top }{\underline{u}}_i\right)+\sum _{i=1}^n{x}_{1i}{x}_{2i}{u}_{ji}\frac{exp\left({\underline{\gamma }}^{\top }{\underline{u}}_i\right)}{{\lambda }_{2}exp\left({\underline{\beta }}^{\top }{\underline{u}}_i\right)+{\lambda }_{3}exp\left({\underline{\gamma }}^{\top }{\underline{u}}_i\right)x_{1i}}=0.\end{array}$$

(10)

In particular, consider $d=1$ and let $u_1,\cdots ,u_n$ be the observed covariates. The likelihood equations from (5) to (10) simplify to become (with notation ${\alpha }_1=\alpha$, ${\beta }_1=\beta$ and ${\gamma }_1=\gamma$)

$$\begin{array}{ccc}\hfill {\lambda }_1\sum _{i=1}^{n}exp(\alpha u_i)& =& \hfill \sum _{i=1}^n{x}_{1i}\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill \sum _{i=1}^{n}exp(u_i\beta )& =& \hfill \sum _{i=1}^n{x}_{2i}\frac{1}{{\lambda }_2+{\lambda }_{3}exp(u_i(\gamma -\beta ))x_{1i}}\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill \sum _{i=1}^n{x}_{1i}exp(u_i\gamma )& =& \hfill \sum _{i=1}^n{x}_{1i}{x}_{2i}\frac{1}{{\lambda }_{2}exp(u_i(\beta -\gamma ))+{\lambda }_3{x}_{1i}}.\hfill \end{array}$$

(13)

In the same way,

$$\begin{array}{ccc}\hfill {\lambda }_1\sum _{i=1}^n{u}_{i}exp\left(\alpha u_i\right)& =& \hfill \sum _{i=1}^n{x}_{1i}{u}_i\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill \sum _{i=1}^n{u}_{i}exp\left(\beta u_i\right)& =& \hfill \sum _{i=1}^n{x}_{2i}{u}_i\frac{1}{{\lambda }_2+{\lambda }_3{x}_{1i}exp(u_i(\gamma -\beta ))}\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill \sum _{i=1}^n{x}_{1i}{u}_{i}exp\left(\gamma u_i\right)& =& \hfill \sum _{i=1}^n{x}_{1i}{x}_{2i}{u}_i\frac{1}{{\lambda }_{2}exp(u_i(\beta -\gamma ))+{\lambda }_3{x}_{1i}}.\hfill \end{array}$$

(16)

Note that the equations from (11) to (16) do not yield explicit expressions for the maximum likelihood estimates. However, one can use numerical methods to solve the system of six equations with six unknown parameters.

4.2. Likelihood Ratio Test

The general form of a generalized likelihood ratio test statistic is of the form

$$\begin{array}{c}\hfill \Lambda =\frac{{sup}_{\theta \in {\Theta }_0}L(\theta )}{{sup}_{\theta \in \Theta }L(\theta )}\end{array}$$

(17)

Here, ${\Theta }_0$ is a subset of $\Theta$ and we envision testing $H_0:\theta \in {\Theta }_0$. We reject the null hypothesis for small values of $\Lambda$.

Now, for the bivariate pseudo-Poisson model, the natural parameter space under the full model is $\Theta =\{{({\lambda }_1,{\lambda }_2,{\lambda }_3,\underline{\alpha },\underline{\beta },\underline{\gamma })}^{\top }:{\lambda }_1\ge 0,{\lambda }_2\ge 0,{\lambda }_3\ge 0,\underline{\alpha }\in {\mathbb{R}}^d,\underline{\beta }\in {\mathbb{R}}^d,\underline{\gamma }\in {\mathbb{R}}^d\}$. The m.l.e.’s under the complete parameter space are obtained by taking partial differentiation of Equation (4) with respect to ${\lambda }_1,{\lambda }_2,{\lambda }_3,\underline{\alpha },\underline{\beta },\underline{\gamma }$ and equating to zero. We denote the obtained numerical solution m.l.e.’s by $\widehat{{\lambda }_1},\widehat{{\lambda }_2},\widehat{{\lambda }_3}$ and $\widehat{\underline{\alpha }}$, $\widehat{\underline{\beta }}$, $\widehat{\underline{\gamma }}$.

Remark 1.

We used the “maxLik” optimization function in R (in the package “maxLik”) to obtain the m.l.e.’s by numerical solution. This function also allows us to use a different methods of optimization using algorithms such as Newton–Raphson, Broyden–Fletcher–Goldfarb–Shanno, Berndt–Hall–Hall–Hausman, Berndt–Hall–Hall–Hausman, Simulated Annealing, Conjugate Gradients and Nelder–Mead methods. In the current paper, we use the Newton–Raphson method to estimate parameters and to compute their standard errors.

4.2.1. Testing $H_0:\underline{\alpha }=\underline{\beta }=\underline{\gamma }=\underline{0}$

In the following, we will construct a likelihood ratio test for testing whether the observed concomitant does not affect the distribution of $(X_1,X_2)$. Under the null hypothesis, the natural parameter space is ${\Theta }_0=\{{({\lambda }_1,{\lambda }_2,{\lambda }_3,\underline{\alpha },\underline{\beta },\underline{\gamma })}^{\top }:{\lambda }_1>0,{\lambda }_2\ge 0,{\lambda }_3\ge 0,\underline{\alpha }=\underline{0},\underline{\beta }=\underline{0},\underline{\gamma }=\underline{0}\}$. Now, taking partial derivatives of Equation (4) with respect to each parameters ${\lambda }_1$, ${\lambda }_2$, ${\lambda }_3$ and equating to zero yields

$$\begin{array}{ccc}\hfill -n+\frac{1}{{\lambda }_1}\sum _{i=1}^n{X}_{1i}& =& \hfill 0\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill -n+\sum _{i=1}^n\frac{X_{2i}}{{\lambda }_2+{\lambda }_3{X}_{1i}}& =& \hfill 0\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill -\sum _{i=1}^n{X}_{1i}+\sum _{i=1}^n\frac{X_{1i}{X}_{2i}}{{\lambda }_2+{\lambda }_3{X}_{1i}}& =& \hfill 0.\hfill \end{array}$$

(20)

Equation (18) is readily solved, to obtain the m.l.e. for ${\lambda }_1$, namely, $\widehat{{\lambda }_1^{*}}={\overline{X}}_1$. The remaining two Equations (19) and (20) must be solved numerically to obtain $\widehat{{\lambda }_2^{*}}$, $\widehat{{\lambda }_3^{*}}$.

Now let $\widehat{{\lambda }_1},\widehat{{\lambda }_2},\widehat{{\lambda }_3}$ and $\widehat{\underline{\alpha }}$, $\widehat{\underline{\beta }}$, $\widehat{\underline{\gamma }}$ be the m.l.e. estimates on unrestricted space. Then, the likelihood (as defined in Equation (4)) ratio test statistic is

$$\begin{array}{c}\hfill {\Lambda }_1=\frac{L(\widehat{{\lambda }_1^{*}},\widehat{{\lambda }_2^{*}},\widehat{{\lambda }_3^{*}},\underline{0},\underline{0},\underline{0})}{L(\widehat{{\lambda }_1},\widehat{{\lambda }_2},\widehat{{\lambda }_3},\widehat{\underline{\alpha }},\widehat{\underline{\beta }},\widehat{\underline{\gamma }})}.\end{array}$$

(21)

If n is large, then $-2log({\Lambda }_1)$ may be compared with a suitable ${\chi }_{3d}^2$ percentile in order to decide whether $H_0$ should be rejected or not.

4.2.2. Testing $H_0:\underline{\alpha }=\underline{0}$

Here, we are testing that the observed concomitant does not affect the marginal distribution of $X_1$. Note that under the null hypothesis, the natural parameter space is ${\Theta }_0=\{{({\lambda }_1,{\lambda }_2,{\lambda }_3,\underline{\alpha },\underline{\beta },\underline{\gamma })}^{\top }:{\lambda }_1>0,{\lambda }_2\ge 0,{\lambda }_3\ge 0,\underline{\alpha }=\underline{0},\underline{\beta }\in {\mathbb{R}}^d,\underline{\gamma }\in {\mathbb{R}}^d\}$. Now, again taking partial derivatives of Equation (4) with respect to parameters ${\lambda }_1$,${\lambda }_2$, ${\lambda }_3$ & $\underline{\beta }$, $\underline{\gamma }$ and equating to zero gives m.l.e.’s, denoted by $\widehat{{\lambda }_1^{*}}$, $\widehat{{\lambda }_2^{*}}$, $\widehat{{\lambda }_3^{*}}$, ${\widehat{\underline{\beta }}}^{*}$ and ${\widehat{\underline{\gamma }}}^{*}$, respectively. The likelihood ratio test statistic is

$$\begin{array}{c}\hfill {\Lambda }_2=\frac{L(\widehat{{\lambda }_1^{*}},\widehat{{\lambda }_2^{*}},\widehat{{\lambda }_3^{*}},\underline{0},{\widehat{\underline{\beta }}}^{*}{\widehat{\underline{\gamma }}}^{*})}{L(\widehat{{\lambda }_1},\widehat{{\lambda }_2},\widehat{{\lambda }_3},\widehat{\underline{\alpha }},\widehat{\underline{\beta }},\widehat{\underline{\gamma }})}\end{array}$$

(22)

If n is large, then $-2log({\Lambda }_2)$ may be compared with a suitable ${\chi }_d^2$ percentile in order to decide whether $H_0$ should be rejected or not.

4.2.3. Testing $H_0:\underline{\beta }=\underline{\gamma }=\underline{0}$

In this case, we are testing whether the observed concomitant does not affect the conditional distribution of $X_2$ given $X_1$. Under the null hypothesis, the natural parameter space is ${\Theta }_0=\{{({\lambda }_1,{\lambda }_2,{\lambda }_3,\underline{\alpha },\underline{\beta },\underline{\gamma })}^{\top }:{\lambda }_1>0,{\lambda }_2\ge 0,{\lambda }_3\ge 0,\underline{\alpha }\in {\mathbb{R}}^d,\underline{\beta }=\underline{0},\underline{\gamma }=\underline{0}\}$. Again, taking partial derivatives of Equation (4) with respect to each parameters ${\lambda }_1$, ${\lambda }_2$, ${\lambda }_3$ & $\underline{\alpha }$ and equating to zero gives to m.l.e.’s denoted by $\widehat{{\lambda }_1^{*}}$, $\widehat{{\lambda }_2^{*}}$, $\widehat{{\lambda }_3^{*}}$, ${\widehat{\underline{\alpha }}}^{*}$. The likelihood ratio test statistic is

$$\begin{array}{c}\hfill {\Lambda }_3=\frac{L(\widehat{{\lambda }_1^{*}},\widehat{{\lambda }_2^{*}},\widehat{{\lambda }_3^{*}},{\widehat{\underline{\alpha }}}^{*},\underline{0},\underline{0})}{L(\widehat{{\lambda }_1},\widehat{{\lambda }_2},\widehat{{\lambda }_3},\widehat{\underline{\alpha }},\widehat{\underline{\beta }},\widehat{\underline{\gamma }})}\end{array}$$

(23)

If n is large, then $-2log({\Lambda }_3)$ may be compared with a suitable ${\chi }_{2d}^2$ percentile in order to decide whether $H_0$ should be rejected or not.

4.2.4. Testing $H_0:\underline{\beta }=\underline{0}$

Here, we are interested in testing whether the observed concomitant does not affect the intercept term of the pseudo-Poisson model. Now, under the null hypothesis, the natural parameter space is ${\Theta }_0=\{{({\lambda }_1,{\lambda }_2,{\lambda }_3,\underline{\alpha },\underline{\beta },\underline{\gamma })}^{\top }:{\lambda }_1>0,{\lambda }_2\ge 0,{\lambda }_3\ge 0,\underline{\alpha }\in {\mathbb{R}}^d,\underline{\beta }=\underline{0},\underline{\gamma }\in {\mathbb{R}}^d\}$. Again, taking partial derivatives of Equation (4) with respect to each parameters ${\lambda }_1$,${\lambda }_2$, ${\lambda }_3$ & $\underline{\alpha }$, $\underline{\gamma }$ and equating to zero gives to m.l.e.’s denoted by $\widehat{{\lambda }_1^{*}}$, $\widehat{{\lambda }_2^{*}}$, $\widehat{{\lambda }_3^{*}}$, ${\widehat{\underline{\alpha }}}^{*}$, ${\widehat{\underline{\gamma }}}^{*}$. The likelihood ratio test statistic is

$$\begin{array}{c}\hfill {\Lambda }_4=\frac{L(\widehat{{\lambda }_1^{*}},\widehat{{\lambda }_2^{*}},\widehat{{\lambda }_3^{*}},{\widehat{\underline{\alpha }}}^{*},\underline{0},{\widehat{\underline{\gamma }}}^{*})}{L(\widehat{{\lambda }_1},\widehat{{\lambda }_2},\widehat{{\lambda }_3},\widehat{\underline{\alpha }},\widehat{\underline{\beta }},\widehat{\underline{\gamma }})}\end{array}$$

(24)

If n is large, then $-2log({\Lambda }_4)$ may be compared with a suitable ${\chi }_d^2$ percentile in order to decide whether $H_0$ should be rejected or not.

4.2.5. Testing $H_0:\underline{\gamma }=\underline{0}$

In this case, we wish to determine whether the concomitant does not affect the dependence structure of the pseudo-Poisson model. Thus, under the null hypothesis, parameter space is ${\Theta }_0=\{{({\lambda }_1,{\lambda }_2,{\lambda }_3,\underline{\alpha },\underline{\beta },\underline{\gamma })}^{\top }:{\lambda }_1>0,{\lambda }_2\ge 0,{\lambda }_3\ge 0,\underline{\alpha }\in {\mathbb{R}}^d,\underline{\beta }\in {\mathbb{R}}^d,\underline{\gamma }=\underline{0}\}$. Now, taking partial derivatives of Equation (4) with respect to the parameters ${\lambda }_1$,${\lambda }_2$, ${\lambda }_3$ & $\underline{\alpha }$, $\underline{\beta }$ and equating to zero gives, m.l.e.’s denoted by $\widehat{{\lambda }_1^{*}}$, $\widehat{{\lambda }_2^{*}}$, $\widehat{{\lambda }_3^{*}}$, ${\widehat{\underline{\alpha }}}^{*}$, ${\widehat{\underline{\beta }}}^{*}$. The likelihood ratio test statistic is

$$\begin{array}{c}\hfill {\Lambda }_5=\frac{L(\widehat{{\lambda }_1^{*}},\widehat{{\lambda }_2^{*}},\widehat{{\lambda }_3^{*}},{\widehat{\underline{\alpha }}}^{*},{\widehat{\underline{\beta }}}^{*},\underline{0})}{L(\widehat{{\lambda }_1},\widehat{{\lambda }_2},\widehat{{\lambda }_3},\widehat{\underline{\alpha }},\widehat{\underline{\beta }},\widehat{\underline{\gamma }})}\end{array}$$

(25)

If n is large, then $-2log({\Lambda }_5)$ may be compared with a suitable ${\chi }_d^2$ percentile in order to decide whether $H_0$ should be rejected or not.

In the next examples, we are interested in testing some hypotheses concerning the relationship between the explanatory and response variables. In particular, we are interested in testing whether the regression planes are parallel or if they are coincident. We illustrate the testing procedure using the simple sub-model given by

$$X_1\sim \mathcal{P}\left(exp\left(\sum _{j=1}^d{u}_{ij}{\alpha }_j\right)\right)$$

(26)

and

$$X_2|X_1=x_1\sim \mathcal{P}\left(exp\left(\sum _{j=1}^d{u}_{ij}{\gamma }_j\right)x_1\right).$$

(27)

4.2.6. Testing for Parallelism

In the following, we are interested in testing whether the planes on which the means lie are parallel. If we set $u_{1i}=1$ for $i\in \{1,\cdots ,n\}$ then the two marginal means are

$$\begin{array}{ccc}\hfill log(E(X_1))& =& \hfill {\alpha }_1+\sum _{j=2}^d{u}_{ij}{\alpha }_j\hfill \end{array}$$

(28)

$$\begin{array}{ccc}\hfill log(E(X_2))& =& \hfill {\alpha }_1+{\gamma }_1+\sum _{j=2}^d{u}_{ij}({\alpha }_j+{\gamma }_j).\hfill \end{array}$$

(29)

For the bivariate pseudo-Poisson regression model specified in (26) and (27), now it is interesting to examine the hypothesis that the planes on which the mean lies are parallel. This is equivalent to testing for the hypothesis $H_0:{\gamma }_j=0$, for $j\in \{1,\cdots ,d\}$. Under the null hypothesis, the pseudo-Poisson regression model will be

$$X_1\sim \mathcal{P}\left(exp\left({\alpha }_1+\sum _{j=2}^d{u}_{ij}{\alpha }_j\right)\right)$$

(30)

and

$$X_2|X_1=x_1\sim \mathcal{P}\left(exp({\gamma }_1)x_1\right).$$

(31)

The log-likelihood is

$$\begin{array}{ccc}\hfill log(L_{PH})& =& \hfill -\sum _{i=1}^{n}exp\left({\alpha }_1+\sum _{j=2}^d{u}_{ij}{\alpha }_j\right)+\sum _{i=1}^n{x}_{1i}log\left(exp\left({\alpha }_1+\sum _{j=2}^d{u}_{ij}{\alpha }_j\right)\right)\hfill \\ & & \hfill -\sum _{i=1}^{n}exp({\gamma }_1)x_{1i}+\sum _{i=1}^n{x}_{2i}log(exp({\gamma }_1)x_{1i})-\sum _{i=1}^{n}log(x_{1i}!x_{2i}!).\hfill \end{array}$$

Note that testing for parallelism for the model specified in (30) and (31) is equivalent to testing for the observed concomitant and has no effect on the conditional distribution of $X_2$ given $X_1$. Now, partial differentiation with respect to ${\gamma }_1$ and ${\alpha }_j$, $j\in \{1,\cdots ,d\}$ and equating to zero gives us

$$\begin{array}{ccc}\hfill log({\overline{X}}_1)& =& \hfill {\alpha }_1+\sum _{j=2}^d{u}_{ij}{\alpha }_j\hfill \\ \hfill \sum _{i=1}^n{x}_{1i}{u}_{ij}& =& \hfill \sum _{i=1}^{n}exp\left({\alpha }_1+\sum _{j=2}^d{\alpha }_j{u}_{ij}\right),j\in \{2,\cdots ,d\}\hfill \end{array}$$

(32)

Solving the above d equations leads us to the m.l.e. of ${\alpha }_i$ denoted by $\widehat{{\alpha }_{Pj}^{*}}$, $j\in \{1,\cdots ,d\}$ and the m.l.e. of ${\gamma }_1$ is

$$\begin{array}{c}\hfill \widehat{{\gamma }_{P1}^{*}}=log\left(\sum _{i=1}^n{x}_{2i}log(x_{1i})-\sum _{i=1}^n{x}_{1i}\right).\end{array}$$

(33)

Now, we denote the obtained m.l.e.’s under the complete parameter space by ${\widehat{\alpha }}_{Pj}$ and ${\widehat{\gamma }}_{Pj}$, $j\in \{1,\cdots ,d\}$. The likelihood ratio test statistic is

$$\begin{array}{c}\hfill {\Lambda }_P=\frac{L_{PH}(\underline{\widehat{{\alpha }_P^{*}}},\underline{0})}{L_P(\underline{\widehat{{\alpha }_P}},\underline{\widehat{{\gamma }_P}})},\end{array}$$

where $L_P(.,.)$ is the likelihood of the model in (26) & (27) and (30) & (31). If n is large, then $-2log({\Lambda }_P)$ may be compared with a suitable ${\chi }_{d-1}^2$ percentile in order to decide whether $H_0$ should be rejected or not.

4.2.7. Testing for Coincidence

Here, we assume that the regression relationship does not change from time 1 to time 2 which will occur if the planes on which means lies are coincident. Now, for the given model in $(26)$ and $(27)$, the two marginal means are

$$\begin{array}{ccc}\hfill log(E(X_1))& =& \hfill \sum _{j=1}^d{u}_{ij}{\alpha }_j\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill log(E(X_2))& =& \hfill \sum _{j=1}^d{u}_{ij}({\alpha }_j+{\gamma }_j).\hfill \end{array}$$

(35)

The assumption of coincidence leads us to test $H_0:{\gamma }_j=0$, for $j\in \{1,\cdots ,d\}$. Denote by ${\widehat{\alpha }}_{Cj}^{*}$ for $j\in \{1,\cdots ,d\}$ are m.l.e.’s under the null hypothesis and by ${\widehat{\alpha }}_{Cj}$ and ${\widehat{\gamma }}_{Cj}$ for $j\in \{1,\cdots ,d\}$ are m.l.e.’s under complete parameter space, for . Now, the likelihood ratio test statistic is

$$\begin{array}{c}\hfill {\Lambda }_C=\frac{L_{CH}(\underline{\widehat{{\alpha }_P^{*}}},\underline{0})}{L_C(\underline{\widehat{{\alpha }_P}},\underline{\widehat{{\gamma }_P}})},\end{array}$$

where $L_{CH}(.,)$ and $L_C(.,.)$ are likelihood under null and complete parameter space, respectively. If n is large, then $-2log({\Lambda }_C)$ may be compared with a suitable ${\chi }_d^2$ percentile to decide whether $H_0$ should be rejected or not.

5. Applications

In the following two subsections, we illustrate a simulation study and give examples of real-life applications of the bivariate pseudo-Poisson regression model.

5.1. Simulation

We have simulated 2000 data sets of sample size $n=20,30,50,100,200,500,$ 1000 for the parameter values ${\lambda }_1=1$, ${\lambda }_2=1$, ${\lambda }_3=4$, ${\alpha }_1=1$, ${\alpha }_2=0$, ${\alpha }_3=-1$, ${\beta }_1=0$, ${\beta }_2=1$, ${\beta }_3=1$, ${\gamma }_1=0$, ${\gamma }_2=0$ and ${\gamma }_3=1$ from the pseudo-Poisson regression model. We refer to Figure 1, Figure 2, Figure 3 and Figure 4 for the bootstrapped distribution of each of the parameters. The numerical evidence suggests that as sample size increases, m.l.e.’s approach the true parameter values with standard errors that are decreasing as the sample size increases.

5.2. Real-Life Data

5.2.1. Australian Health Service Utilization Data: 1977–1978

We consider a data set which is mentioned in Islam and Chowdhury [14] that is part of the Health and Retirement Study (HRS). The data represent the number of conditions ever had $(X_1)$ as mentioned by the doctors and utilization of healthcare services (say, hospital, nursing home, doctor and home care) $(X_2)$. The concomitant variables are Gender, Age, Hispanic, and Veteran.

The marginal estimated dispersion indices are $0.779$ and $1.029$. The sample Pearson correlation coefficient between $X_1$ and $X_2$ is $0.063$. We can conclude that marginal $X_1$ is approximately equi-dispersed and marginal $X_2$ is slightly over-dispersed. Further, the data were also tested for independence and it was concluded that the assumption was rejected, cf. Arnold and Manjunath [7] pages 2321–2322.

We refer to the

Table 1. Models for the Australian Health Service Utilization Data.

Models	No. Parameters	Log-Likelihood
Full Model	15	−52,654.13
Mirrored, Full Model	15	−16,229.24
Sub-Model I (${\lambda }_2=0$)	11	−16,586.07
Sub-Model II (${\lambda }_2={\lambda }_3$)	11	−16,371.95
Mirrored Sub-Model II (${\lambda }_2={\lambda }_3$)	11	−17,585.37

for the log-likelihood values for the following models:

Full Model:

The parameters are ${\lambda }_1$, ${\lambda }_2$, ${\lambda }_3$, $\underline{\alpha }={({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4)}^{\top }$,

$\underline{\beta }={({\beta }_1,{\beta }_2,{\beta }_3,{\beta }_4)}^{\top }$, $\underline{\gamma }={({\gamma }_1,{\gamma }_2,{\gamma }_3,{\gamma }_4)}^{\top }$

Mirrored, Model

(in which $X_1$ and $X_2$ are interchanged): The parameters are ${\lambda }_1$, ${\lambda }_2$, ${\lambda }_3$,

$\underline{\alpha }={({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4)}^{\top }$, $\underline{\beta }={({\beta }_1,{\beta }_2,{\beta }_3,{\beta }_4)}^{\top }$, $\underline{\gamma }={({\gamma }_1,{\gamma }_2,{\gamma }_3,{\gamma }_4)}^{\top }$

Sub-Model I

(${\lambda }_2=0$): The parameters are ${\lambda }_1$,${\lambda }_3$, $\underline{\alpha }={({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4)}^{\top }$,

$\underline{\gamma }={({\gamma }_1,{\gamma }_2,{\gamma }_3,{\gamma }_4)}^{\top }$

Sub-Model II

(${\lambda }_2={\lambda }_3$): The parameters are ${\lambda }_1$,${\lambda }_3$, $\underline{\alpha }={({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4)}^{\top }$,

$\underline{\gamma }={({\gamma }_1,{\gamma }_2,{\gamma }_3,{\gamma }_4)}^{\top }$

Sub-Model II

(Mirrored): The parameters are ${\lambda }_1$,${\lambda }_3$, $\underline{\alpha }={({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4)}^{\top }$,

$\underline{\gamma }={({\gamma }_1,{\gamma }_2,{\gamma }_3,{\gamma }_4)}^{\top }$.

The mirrored Full Model fits the data best. For the detailed discussion on the mirrored model, see Arnold and Manjunath [7] page 2323. In Islam and Chowdhury [14], page 122, the authors fitted the Poisson–Poisson regression model for the same data set. Note that the Poisson–Poisson regression model is a sub-model of the pseudo-Poisson regression model when ${\lambda }_2=0$. Hence, we conclude that our generalized pseudo-Poisson mirrored model fits the data better than the Poisson–Poisson regression model. The parameter estimates for the pseudo-Poisson mirrored model and their standard errors are displayed in

Table 2. Final model estimates and its standard error (s.e.) for the Australian Health Service Utilization Data.

Parameter	m.l.e.	s.e.
${\alpha }_1$	$0.292$	$0.039$
${\alpha }_2$	$-0.008$	$0.004$
${\alpha }_3$	$-0.186$	$0.058$
${\alpha }_4$	$0.140$	$0.042$
${\beta }_1$	$-0.132$	$0.0273$
${\beta }_2$	$0.016$	$0.0036$
${\beta }_3$	$0.038$	$0.0277$
${\beta }_4$	$0.053$	$0.035$
${\gamma }_1$	$1.636$	$0.656$
${\gamma }_2$	$-0.025$	$0.039$
${\gamma }_3$	$-0.996$	−
${\gamma }_4$	$-0.148$	$0.273$
${\lambda }_1$	$1.172$	$0.385$
${\lambda }_2$	$0.824$	$0.224$
${\lambda }_3$	$0.313$	$1.037$

.

Further, we tested for the significance of the regression coefficients. With reference to

Table 3. Hypothesis testing for the Australian Health Service Utilization Data.

Hypothesis	$log{\Lambda }^{*}-log\Lambda$	d.f.
$\underline{\alpha }=\underline{\beta }=\underline{\gamma }=\underline{0}$	$-113.7227$	12
$\underline{\alpha }=\underline{0}$	$-103.483$	4
$\underline{\beta }=\underline{\gamma }=\underline{0}$	$-28.26604$	8
$\underline{\beta }=\underline{0}$	$-24.26104$	4
$\underline{\gamma }=\underline{0}$	$-6.857175$	4

, the computed $-2log\lambda$ and compared with ${\chi }^2$ table values with respective degrees of freedom and the size of $0.05$ or $0.10$ and concluded that there is not enough evidence to accept the null hypotheses.

5.2.2. Road Safety Data

The second data set is on road safety, published by the Department for Transport, United Kingdom. The data comprise information about personal injury road accidents in Great Britain and the consequent casualties on public roads. The concomitant variables are Gender of the driver (Male = 1, Female = 0), Area (Urban = 0, Rural = 1), Accident Severity (Fatal Severity = 1 else = 0), Accident Severity ( Serious Severity = 1, else = 0), and Light condition (Daylight = 1, Others = 0).

We refer to

Table 4. Models for the Road safety data.

Models	No. Parameters	Log-Likelihood
Full Model	18	−223,743.3
Mirrored Full Model	18	−243,538.7
Sub-Model I(${\lambda }_2=0$)	11	−251,937.1
Mirrored Sub-Model I(${\lambda }_2=0$)	11	−37,599.63
Sub-Model II(${\lambda }_2={\lambda }_3$)	11	−36,201.52
Mirrored Sub-Model II(${\lambda }_2={\lambda }_3$)	11	−36,516.22

for the log-likelihood values for the following:

Full Model:

parameters are ${\lambda }_1$, ${\lambda }_2$, ${\lambda }_3$, $\underline{\alpha }={({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5)}^{\top }$,

$\underline{\beta }={({\beta }_1,{\beta }_2,{\beta }_3,{\beta }_4,{\beta }_5)}^{\top }$, $\underline{\gamma }={({\gamma }_1,{\gamma }_2,{\gamma }_3,{\gamma }_4,{\gamma }_5)}^{\top }$

Mirrored, Model

($X_1$ and $X_2$ are interchanged): parameters are ${\lambda }_1$, ${\lambda }_2$, ${\lambda }_3$,

$\underline{\alpha }={({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5)}^{\top }$, $\underline{\beta }={({\beta }_1,{\beta }_2,{\beta }_3,{\beta }_4,{\beta }_5)}^{\top }$, $\underline{\gamma }={({\gamma }_1,{\gamma }_2,{\gamma }_3,{\gamma }_4,{\gamma }_5)}^{\top }$

Sub-Model I

(${\lambda }_2=0$): parameters are ${\lambda }_1$,${\lambda }_3$, $\underline{\alpha }={({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5)}^{\top }$,

$\underline{\gamma }={({\gamma }_1,{\gamma }_2,{\gamma }_3,{\gamma }_4,{\gamma }_5)}^{\top }$

Sub-Model I

(Mirrored): parameters are ${\lambda }_1$,${\lambda }_3$, $\underline{\alpha }={({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5)}^{\top }$,

$\underline{\gamma }={({\gamma }_1,{\gamma }_2,{\gamma }_3,{\gamma }_4,{\gamma }_5)}^{\top }$

Sub-Model II

(${\lambda }_2={\lambda }_3$): parameters are ${\lambda }_1$,${\lambda }_3$, $\underline{\alpha }={({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5)}^{\top }$,

$\underline{\gamma }={({\gamma }_1,{\gamma }_2,{\gamma }_3,{\gamma }_4,{\gamma }_5)}^{\top }$

Sub-Model II

(Mirrored): parameters are ${\lambda }_1$,${\lambda }_3$, $\underline{\alpha }={({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5)}^{\top }$,

$\underline{\gamma }={({\gamma }_1,{\gamma }_2,{\gamma }_3,{\gamma }_4,{\gamma }_5)}^{\top }$.

We refer to Table 4 and conclude that the Full Model fits the road safety data and refer to

Table 5. Final model estimates and its standard error (s.e.) Road safety data.

Parameter	m.l.e.	s.e.
${\alpha }_1$	$1.002$	$0.006$
${\alpha }_2$	$0.999$	$0.005$
${\alpha }_3$	$0.999$	$0.017$
${\alpha }_4$	$0.999$	$0.005$
${\beta }_1$	$1.000$	$0.005$
${\beta }_2$	$1.004$	$0.005$
${\beta }_3$	$1.003$	$0.005$
${\beta }_4$	$1.000$	$0.015$
${\gamma }_1$	$0.999$	$0.0036$
${\gamma }_2$	$1.005$	$0.004$
${\gamma }_3$	$1.355$	–
${\gamma }_4$	$1.319$	–
${\lambda }_1$	$1.010$	–
${\lambda }_2$	$1.105$	–
${\lambda }_3$	$-0.078$	$0.007$

for the estimates and their standard errors.

6. Concluding Remarks

The bivariate pseudo-Poisson model with its straightforward structure with no restrictions on the conditional mean function allows us to model a variety of dependence structures, including positive and negative correlation. Introducing explanatory variables in such models will be a useful additional to the toolkit for modelers dealing with bivariate count data which have positive or negative correlation. In the current note, we explored distributional and inferential aspects of such models and also included a simulation and real-life data applications. We emphasize the advantage of considering the current model over other available count regression models in Section 5.2. The bivariate pseudo-Poisson regression model has a simple structure, straightforward parameter estimation and fast computation, and will deserve a place in the analysis of count data sets with concomitants.