Paired count regressions for modeling the number of doctor consultations and non-prescribed drugs intake

Abstract

There are sundry practical situations in which paired count variables are correlated, thus requiring a joint estimation method. In this article, we introduce a flexible class of bivariate mixed Poisson regression models, which settle into an exponential-family (EF) distributed component for unobserved heterogeneity. The proposed bivariate mixed Poisson models deal with the phenomenon of overdispersion, typical of count data, and have flexibility in terms of the correlation structure. Thus, this novel class of models has a distinct advantage over the most widely used models because it captures both positive and negative correlations in the count data. Under the bivariate mixed Poisson model, inference of the parameters is conducted through the maximum likelihood method. Monte Carlo studies on assessing the finite-sample performance of the estimators of the parameters are presented. Furthermore, we employ a likelihood ratio statistic for testing the significance of certain sources of correlation and evaluate its performance via simulation studies. Moreover, model adequacy is addressed by using simulated envelopes for residual analysis, and also a randomized probability integral transformation for calibration model control. The proposed bivariate mixed Poisson model is considered for analyzing a healthcare dataset from the Australian Health Survey, where our aim is to study the association between the number of consultations with a doctor and the number of non-prescribed drug intake.

1. Introduction

Bivariate count data are often dependent, especially for those that come from clinical, healthcare, and social surveys. In this article, we develop a model for analyzing counts of medical visits and frequency of intake of non-prescription medication. This is very important to the health insurance system because it is essential to provide a reliable estimate of the number of claim events in an insurance portfolio because this is the starting point for establishing a fair rate premium (an amount payable for desired benefits in an insurance policy). In particular, in the Brazilian healthcare system, actuaries of health insurance companies need to predict the number of doctor appointments and the number of medical exams, among other groups of medical procedures such as hospitalizations, therapies, other ambulatory visits, and dental care. Since a medical exam is realized when the insured has a prescribed medical indication, these variables are naturally correlated and therefore they should be ideally jointly analyzed. Despite the widespread need to analyze associations in count data, the major limitation of the most widely used models is that they do not offer the flexibility of capturing negative correlation. As a consequence, this leads to misleading conclusions and downstream managerial and policy decisions. In this article, we develop a novel model that overcomes this serious limitation, which is demonstrated in the analysis of the Australian healthcare data.

We also point to the presence of unobserved heterogeneity in such type of data which the current models fail to consider. As discussed above, the unobserved heterogeneity perhaps can be explained by portfolio aging, socioeconomic and demographic factors, available technology, epidemiological conditions, and quality of the services provided, among other elements. Sometimes it is possible to evaluate the unobserved heterogeneity by incorporating explanatory variables in the model. However, many of these covariates that account for the unmeasured heterogeneity are difficult to access, such as measuring the quality of service providers. This problem of heterogeneity is addressed by our proposed model.

The bivariate Poisson (BP) distribution proposed by Campbell, ¹ and further explored by Holgate ² and Johnson and Kotz, ³ is an option for jointy modeling count data. The trivariate reduction method ⁴ considers three independent Poisson random variables, yielding the BP distribution as sums of these independent Poisson with a common element. However, its weakness is that it neither allows overdispersion nor negative correlation. For more details and properties of multivariate discrete distributions, we refer to Marshall and Olkin, ⁵ Marshall and Olkin, ⁶ and Cameron and Trivedi, ⁷ to name a few. For bivariate count distributions that allow for negative correlation, we refer to the class proposed by Sarmanov. ⁸ A multivariate extension was developed by Lee. ⁹ According to Lakshminarayana et al., ¹⁰ a BP distribution that allows negative correlation was introduced, which was derived from the product of two marginal Poisson random variables and a multiplicative factor. This family is a member of the Sarmanov family of bivariate distributions.

Concerning the problem of unobserved heterogeneity, Gurmu and Elder ¹¹ proposed a bivariate regression model based on a first-degree polynomial expansion of the gamma distribution. This model produced a generalized bivariate negative binomial (BNB) model, although it only accommodates positive correlation. Later, Gurmu and Elder ¹² introduced a more flexible class of bivariate count regression models allowing for both positive and negative correlations. Another version of the bivariate negative regression model was proposed by Famoye, ¹³ based on a similar framework of the BP distribution introduced by Lakshminarayana et al. ¹⁰ This model is employed in a healthcare dataset and compared to the BP log-normal modeling by Ho and Singer ¹⁴ and Aitchison Ho. ¹⁵ To account for zero inflation in count data, Faroughi and Ismail ¹⁶ suggest a class of nested bivariate zero-inflated negative binomial regression models, derived from the product of two zero-inflated negative binomial distributions with a multiplicative factor, and also demonstrates their application to healthcare.

According to Consul and Famoye, ¹⁷ the bivariate generalized Poisson distribution based on the trivariate reduction method was introduced and the performance of the method of moments and the method of maximum likelihood were compared. Since this model is limited only to positive correlation, Famoye ¹⁸ introduces a new bivariate generalized Poisson distribution, applying the same framework introduced by Lakshminarayana et al. ¹⁰ This model deals not only with positive or negative correlation but also with dispersion phenomena (both under-dispersion and over-dispersion). In this field, Zamani et al. ¹⁹ introduce several forms of bivariate generalized Poisson regression models, which includes the bivariate generalized Poisson model introduced by Famoye ¹⁸ as a particular case.

A bivariate Conway-Maxwell-Poisson (COM-Poisson) distribution was introduced by Sellers et al. ²⁰ It allows for overdispersion in count data and embraces the BP, the bivariate Bernoulli, and the bivariate geometric distributions, as particular cases. Concerning the correlation range, when the model reduces for the bivariate geometric distribution, it can assume a bounded negative correlation for a specific configuration of the dispersion parameter. However, as pointed out by these authors, it is not simple to establish the correlation support for a general level of the precision parameter. Still, in the context of the COM-Poisson distribution, a recent contribution by Piancastelli et al. ²¹ proposed a multivariate count distribution with COM-Poisson marginals by developing a modification on the Sarmanov method for building multivariate distributions.

The main goal of this article is to develop a novel class of bivariate mixed Poisson (BMP) regression models, which takes into account the unobserved heterogeneity through a random effect belonging to the EF. Our models are tailored to handle paired bivariate counts with heterogeneity and flexible correlation structure. We assume a common random effect for the counts, which makes sense to analyze the count of doctor visits and non-prescribed drug intake (paired counts). The advantage of our proposed BMP model, over existing models, is its ability to simultaneously handle these challenges. Next, we outline the main contributions of thisarticle.

(1)

Our class of bivariate count regression models is a multivariate extension of the univariate mixed Poisson (MP) regression models (with constant dispersion/precision) proposed by Barreto-Souza and Simas. ²²

(2)

The BMP regressions have the advantage over some existing models of being flexible about the correlation structure, which combines two sources of dependence and allows for both positive and negative correlations. This will be important especially when we analyze the count doctor consultations and non-prescribed drug intake.

(3)

Our general bivariate count models consolidate a new framework for modeling bivariate count data. Moreover, they contain, as special cases, novel BNB and Poisson-inverse Gaussian regression models; the latter is an important approach for modeling bivariate counts with heavy tails.

(4)

In addition to examining computational aspects for model inference, we also explore residual analysis and a model diagnostic measure.

The rest of the article is organized as follows. In Section 2, we introduce the BMP distributions and discuss, in detail, the model definition and its properties. In Section 3, the estimation method using the maximum likelihood is developed. Moreover, simulation studies aiming at the finite-sample performance evaluation of the maximum likelihood estimators (MLEs) are presented. In Section 4, we discuss hypothesis testing with respect to one parameter associated with the correlation structure. Furthermore, we propose a simulated envelope technique for residual analysis and randomized probability integral transformations as a goodness-of-fit tool. Section 5 is devoted to the statistical analysis of the paired counts doctor consultations and the number of non-prescribed medications variables from the Australian Health Survey (AHS). Concluding remarks and future research are addressed in Section 6. Proofs of propositions are provided in the Appendix. This article contains Supplemental Material with additional numerical results.

2. Model definition and properties

We begin this section by introducing some important concepts for the construction of our class of bivariate count models. Lakshminarayana et al. ¹⁰ proposed a BP distribution allowing for positive and negative correlations. A bivariate vector $(W_1,W_2)$ follows a BP distribution if its joint probability mass function (pmf) assumes the form

$$\begin{array}{rl}p(w_1,w_2)& \equiv P(W_1=w_1,W_2=w_2)={\displaystyle \frac{e^{-({\mu }_1+{\mu }_2)}{\mu }_1^{w_1}{\mu }_2^{w_2}}{w_1!w_2!}}[1+\delta (e^{-w_1}-e^{-{\mu }_{1}c})(e^{-w_2}-e^{-{\mu }_{2}c})]\end{array}$$

(1)

for $(w_1,w_2)\in \mathbb{N}\times \mathbb{N}$ , with $\mathbb{N}\equiv \{0,1,2,\dots \}$ , where $c=1-e^{-1}$ , $\delta$ is a parameter controlling the dependence with parameter space given in (2), and ${\mu }_1,{\mu }_2>0$ . The model in (1) belongs to the family of bivariate distributions introduced by Sarmanov ⁸ (for instance, see discussion in Subsection 2.2 from Piancastelli et al. ²³ ), and the range of $\delta$ can be assessed using Corollary 2 provided by Lee, ⁹ which is given by the following equation:

$$-{\displaystyle \frac{1}{max\{e^{-({\mu }_1+{\mu }_2)c},(1-e^{-{\mu }_{1}c})(1-e^{-{\mu }_{2}c})\}}}<\delta <{\displaystyle \frac{1}{max\{e^{-{\mu }_{1}c}(1-e^{-{\mu }_{2}c}),e^{-{\mu }_{2}c}(1-e^{-{\mu }_{1}c})\}}}$$

(2)

Along with this article, we consider the notation $(W_1,W_2)\sim \text{BP}({\mu }_1,{\mu }_2,\delta )$ . The marginals of the BP model are Poisson distributed and the independent case is obtained by taking $\delta =0$ . The joint moment generating function of $(W_1,W_2)$ , denoted ${\mathrm{\Phi }}_{W_1,W_2}(s,t)\equiv E(e^{s{W}_1+t{W}_2})$ , and the mean vector and covariance matrix are, respectively, derived to be

$$\begin{array}{r}\begin{array}{rl}& {\mathrm{\Phi }}_{W_1,W_2}(s,t)=e^{{\mu }_1(e^s-1)+{\mu }_2(e^t-1)}+\delta [e^{{\mu }_1(e^{s-1}-1)}-e^{{\mu }_1(e^s-1-c)}][e^{{\mu }_2(e^{t-1}-1)}-e^{{\mu }_2(e^t-1-c)}],s,t\in \mathbb{R}\\ & E\left[\begin{array}{c}{W}_1\\ W_2\end{array}\right]=\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\end{array}\right],\text{and}\mathrm{\Sigma }=\left[\begin{array}{cc}{\mu }_1& \delta {\mu }_1{\mu }_2{c}^2{e}^{-({\mu }_1+{\mu }_2)c}\\ \delta {\mu }_1{\mu }_2{c}^2{e}^{-({\mu }_1+{\mu }_2)c}& {\mu }_2\end{array}\right]\end{array}\end{array}$$

(3)

The correlation of the BP distribution is given by $\rho =\mathrm{cov}(W_1,W_2)/\sqrt{\mathrm{Var}(W_1)\mathrm{Var}(W_2)}=\sqrt{{\mu }_1{\mu }_2}\delta c^2{e}^{-({\mu }_1+{\mu }_2)c}$ , that is negatively correlated if $\delta <0$ and positively correlated if $\delta >0$ .

We now present the class of univariate MP distributions proposed by Barreto-Souza ²⁴ and Barreto-Souza and Simas, ²² which combines Poisson distribution and EF. This class will play an important role in this article. Let $Z$ be a continuous positive random variable belonging to the EF with a density function given by the following equation:

$$f(z)=\exp \{\varphi [z{\xi }_0-b({\xi }_0)]+c(z;\varphi )\},z>0$$

(4)

where $\varphi >0$ , $b(\cdot )$ has derivatives of all orders, ${\xi }_0$ is a constant such that $b^{\prime }({\xi }_0)=1$ , with $b^{\prime }(\cdot )$ denoting the first derivative of $b(\cdot )$ , and $c(\cdot ;\cdot )$ is a function mapping ${\mathbb{R}}^{+}\times {\mathbb{R}}^{+}$ into $\mathbb{R}$ . Thus, we denote the unobserved heterogeneity by $Z\sim \text{EF}(\varphi )$ , which has mean and variance given by $E(Z)=b^{\prime }({\xi }_0)=1$ , and $\text{Var}(Z)={\varphi }^{-1}{b}^{″}({\xi }_0)$ , respectively, with $b^{″}(\cdot )$ denoting the second derivative of the function $b(\cdot )$ . The moments of $Z$ can be derived from its moment generating function

$${\mathrm{\Phi }}_Z(t)\equiv E(e^{tZ})=\exp \{-\varphi [b({\xi }_0)-b(\frac{t}{\varphi }+{\xi }_0)]\}$$

(5)

where $t$ belongs to some interval containing 0.

A count random variable $Y$ belongs to the class of univariate MP distributions by Barreto-Souza ²⁴ and Barreto-Souza and Simas ²² if it satisfies the stochastic representation $Y|Z=z\sim \text{Poisson}(\mu z)$ , with $Z\sim \text{EF}(\varphi )$ , where $\mu >0$ and $\varphi >0$ . We now have all the ingredients to define our bivariate class of count distributions.

Definition 2.1

A random vector $(Y_1,Y_2)$ follows a BMP distribution if it satisfies

$$(Y_1,Y_2)|Z=z\sim \text{BP}({\mu }_{1}z,{\mu }_{2}z,\delta )$$

(6)

where $Z\sim \text{EF}(\varphi )$ , for ${\mu }_1,{\mu }_2>0$ , $\varphi >0$ , and $\delta$ satisfying

$${\int }_0^{\mathrm{\infty }}{e}^{-({\mu }_1+{\mu }_2)z}{z}^{y_1+y_2}[1+\delta (e^{-y_1}-e^{-{\mu }_{1}zc})(e^{-y_2}-e^{-{\mu }_{2}zc})]f(z)\mathrm{d}z\ge 0\mathrm{\forall }(y_1,y_2)\in {\mathbb{N}}^2$$

(7)

with $f(\cdot )$ given in (4). We denote $(Y_1,Y_2)\sim \text{BMP}({\mu }_1,{\mu }_2,\delta ,\varphi )$ .

Remark 2.1

Our bivariate model has two sources of correlation. The first is due to the conditional distribution of $(Y_1,Y_2)|Z$ through the parameter $\delta$ . In this regard, we call attention that the constraint in (7) is necessary to ensure that the model is well-defined. The second source of correlation is driven by the random effect $Z$ common to the paired counts $Y_1$ and $Y_2$ , which acts to capture the heterogeneity caused by the non-observation of significant covariates. Such a second source of correlation is controlled by the parameter $\varphi$ .

Remark 2.2

Note that it makes sense in our model construction, which aims to study the count of doctor consultations and non-prescribed drugs intake (to be presented in Section 5), to assume a common random effect $Z$ in (6) since the bivariate response is from the same individual, and therefore possible important missing covariates would be the same to explain both responses.

We obtain an explicit form for the joint probability function of the BMP distributions in the next proposition. In particular, this expression will enable us to perform the maximum likelihood estimation of parameters.

Proposition 2.3

Let $(Y_1,Y_2)\sim \text{BMP}({\mu }_1,{\mu }_2,\delta ,\varphi )$ . The joint pmf of $(Y_1,Y_2)$ can be expressed by the following equation:

$$\begin{array}{rl}p(y_1,y_2)=& {\displaystyle \frac{{\mu }_1^{y_1}{\mu }_2^{y_2}}{y_1!y_2!}}{e}^{-\varphi b({\xi }_0)}\{\mathrm{\Psi }(y_1+y_2;{\mu }_1+{\mu }_2,\varphi )+\delta [e^{-(y_1+y_2)}\mathrm{\Psi }(y_1+y_2;{\mu }_1+{\mu }_2,\varphi )\\ & -e^{-y_1}\mathrm{\Psi }(y_1+y_2;{\mu }_1+{\mu }_2(1+c),\varphi )-e^{-y_2}\mathrm{\Psi }(y_1+y_2;{\mu }_1(1+c)+{\mu }_2,\varphi )\\ & +\mathrm{\Psi }(y_1+y_2;({\mu }_1+{\mu }_2)(1+c),\varphi )]\},\text{for}(y_1,y_2)\in {\mathbb{N}}^2\end{array}$$

(8)

where

$$\mathrm{\Psi }(y;\tau ,\varphi )={\displaystyle \frac{d^y}{d{t}^y}}\exp \left\{\varphi b({\displaystyle \frac{t-\tau }{\varphi }}+{\xi }_0)\right\}{|}_{t=0},\text{for}y\in \mathbb{N}\text{and}\tau >0$$

Next, we provide the moment generating function, a useful tool to explicit the moment structure of the proposed BMP distributions.

Proposition 2.4

The moment generating function of the BMP distributions, say ${\mathrm{\Phi }}_{Y_1,Y_2}(s,t)\equiv E(e^{s{Y}_1+t{Y}_2})$ , and its moment structure is given by the following equation:

$$\begin{array}{rl}{\mathrm{\Phi }}_{Y_1,Y_2}(s,t)& ={\mathrm{\Phi }}_Z[{\mu }_1(e^s-1)+{\mu }_2(e^t-1)]\\ & +\delta \{{\mathrm{\Phi }}_Z[{\mu }_1(e^{s-1}-1)+{\mu }_2(e^{t-1}-1)]-{\mathrm{\Phi }}_Z[{\mu }_1(e^{s-1}-1)+{\mu }_2(e^t-1-c)]\\ & -{\mathrm{\Phi }}_Z[{\mu }_1(e^s-1-c)+{\mu }_2(e^{t-1}-1)]+{\mathrm{\Phi }}_Z[{\mu }_1(e^s-1-c)+{\mu }_2(e^t-1-c)]\}\end{array}$$

$$E\left[\begin{array}{c}{Y}_1\\ Y_2\end{array}\right]=\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\end{array}\right],\mathrm{Var}\left[\begin{array}{c}{Y}_1\\ Y_2\end{array}\right]=\left[\begin{array}{c}{\mu }_1+{\mu }_1^2{\displaystyle \frac{b^{″}({\xi }_0)}{\varphi }}\\ {\mu }_2+{\mu }_2^2{\displaystyle \frac{b^{″}({\xi }_0)}{\varphi }}\end{array}\right],\text{and}$$

$$\begin{array}{rl}& \mathrm{cov}(Y_1,Y_2)={\mu }_1{\mu }_2\delta c^2\exp \left\{\varphi [b({\xi }_0-{\displaystyle \frac{({\mu }_1+{\mu }_2)c}{\varphi }})-b\left({\xi }_0\right)]\right\}\\ & \times \left\{\begin{array}{c}{\displaystyle \frac{b^{″}({\xi }_0-{\displaystyle \frac{({\mu }_1+{\mu }_2)c}{\varphi }})}{\varphi }}+{\left[b^{\prime }({\xi }_0-{\displaystyle \frac{({\mu }_1+{\mu }_2)c}{\varphi }})\right]}^2\end{array}\right\}+{\mu }_1{\mu }_2{\displaystyle \frac{b^{″}({\xi }_0)}{\varphi }}\end{array}$$

where ${\mathrm{\Phi }}_Z(\cdot )$ is the moment generating function stated in (5).

The covariance between $Y_1$ and $Y_2$ can be derived from the moment generating function presented in Proposition 2.4. Alternatively, it can be obtained by using the total covariance law, that is,

$$\begin{array}{rl}\mathrm{cov}(Y_1,Y_2)& =E\left[\mathrm{cov}(Y_1,Y_2|Z)\right]+\mathrm{cov}[E\left(Y_1|Z\right),E\left(Y_2|Z\right)]\\ & =E[E\left(Y_1{Y}_2|Z\right)-E\left(Y_1|Z\right)E\left(Y_2|Z\right)]+\mathrm{cov}({\mu }_{1}Z,{\mu }_{2}Z)\\ & ={\mu }_1{\mu }_2\delta c^{2}E\left[Z^2{e}^{-({\mu }_1+{\mu }_2)cZ}\right]+{\mu }_1{\mu }_{2}V(Z)\\ & ={\mu }_1{\mu }_2\delta c^2\exp \left\{\varphi [b({\xi }_0-{\displaystyle \frac{({\mu }_1+{\mu }_2)c}{\varphi }})-b\left({\xi }_0\right)]\right\}\\ & \times \left\{\begin{array}{c}{\displaystyle \frac{b^{″}({\xi }_0-{\displaystyle \frac{({\mu }_1+{\mu }_2)c}{\varphi }})}{\varphi }}+{\left[b^{\prime }({\xi }_0-{\displaystyle \frac{({\mu }_1+{\mu }_2)c}{\varphi }})\right]}^2\end{array}\right\}+{\mu }_1{\mu }_2{\displaystyle \frac{b^{″}({\xi }_0)}{\varphi }}\end{array}$$

(9)

where $E[Z^2{e}^{-({\mu }_1+{\mu }_2)cZ}]$ can be analytically computed through integration, by identifying the kernel of an EF density function. The correlation coefficient $\eta =\mathrm{cov}(Y_1,Y_2)/\sqrt{\mathrm{Var}(Y_1)\mathrm{Var}(Y_2)}$ is given by the following equation:

$$\begin{array}{rl}\eta & ={\displaystyle \frac{\sqrt{{\mu }_1{\mu }_2}\{\delta c^2\exp \left\{\varphi [b\left({\theta }_0\right)-b\left({\xi }_0\right)]\right\}\{{\displaystyle \frac{b^{″}\left({\theta }_0\right)}{\varphi }}+{\left[b^{\prime }\left({\theta }_0\right)\right]}^2\}+{\displaystyle \frac{b^{″}({\xi }_0)}{\varphi }}\}}{\sqrt{(1+{\displaystyle \frac{{\mu }_1{b}^{″}\left({\xi }_0\right)}{\varphi }})(1+{\displaystyle \frac{{\mu }_2{b}^{″}\left({\xi }_0\right)}{\varphi }})}}}\end{array}$$

(10)

Through moment structure, we notice that the BMP distribution is overdispersed. In fact, for any positive random variable $Z$ , it turns out that the associated MP distribution is overdispersed, as noted by Barreto-Souza and Simas ²² and Gonçalves and Barreto-Souza. ²⁵

Remark 2.5

A small correlation coefficient study between our BMP and BP models is provided in Section 1 of the Supplemental Material.

The conditional BP distribution is straightforward to obtain, as provided by Cui and Zhu. ²⁶ Hence, in what follows, we present its probability mass function and moment generating function since we will use these attributes to derive some results needed to generate samples from our proposed class of BMP distributions. Thus, the probability mass function of the conditional BP distribution is given by the following equation:

$$p(w_2|w_1)={\displaystyle \frac{e^{-{\mu }_2}{\mu }_2^{w_2}}{w_2!}}[1+\delta (e^{-w_1}-e^{-{\mu }_{1}c})(e^{-w_2}-e^{-{\mu }_{2}c})],w_2\in \mathbb{N}$$

(11)

The distribution of $W_1$ conditional on $W_2$ follows similarly. Therefore, since we have computed the conditional distribution in (11), we are able to generate samples of the BMP distributions. For that, we follow analogously to the conditional sampling method proposed by Cui and Zhu, ²⁶ and the steps for generating samples $(y_1,y_2)$ of $(Y_1,Y_2)\sim \text{BMP}({\mu }_1,{\mu }_2,\varphi ,\delta )$ are, in Algorithm 1, described.

In Table 1, we present some quantities associated with two models belonging to the BMP class, the BNB and BP-inverse Gaussian (BPIG) distributions, which arise when the random effect $Z$ is gamma or inverse-Gaussian distributed, respectively.

Table 1.

Quantities associated with gamma and the inverse-Gaussian heterogeneity $Z$ and the resultant distribution for $(Y_1,Y_2)$ .

Distribution of $(Y_1,Y_2)$	Distribution of $Z$	$b(\theta )$	${\xi }_0$	$c(z;\varphi )$
BNB	Gamma	$-\log (-\theta )$	$-1$	$\varphi \log \varphi -\log \mathrm{\Gamma }(\varphi )+(\varphi -1)\log z$
BPIG	Inverse Gaussian	$-\sqrt{-2\theta }$	$-{\displaystyle \frac{1}{2}}$	${\displaystyle \frac{1}{2}}\log ({\displaystyle \frac{\varphi }{2\pi z^3}})-{\displaystyle \frac{\varphi }{2z}}$

BNB: bivariate negative binomial; BPIG: bivariate Poisson-inverse Gaussian.

By using equation (8) from Proposition 2.3 and the quantities provided in Table 1, we obtain explicit forms for the joint pmf of the BNB (equation (12)) and the BPIG (equation (13)) distributions, which are, respectively, given by the following equations:

$$\begin{array}{rl}p(y_1,y_2)& ={\displaystyle \frac{{\mu }_1^{y_1}{\mu }_2^{y_2}{\varphi }^{\varphi }\mathrm{\Gamma }(y_1+y_2+\varphi )}{y_1!y_2!\mathrm{\Gamma }(\varphi )}}\{{\displaystyle \frac{1+\delta e^{-(y_1+y_2)}}{({\mu }_1+{\mu }_2+\varphi {)}^{y_1+y_2+\varphi }}}-{\displaystyle \frac{\delta e^{-y_1}}{[{\mu }_1+{\mu }_2(1+c)+\varphi {]}^{y_1+y_2+\varphi }}}\\ & -{\displaystyle \frac{\delta e^{-y_2}}{[{\mu }_1(1+c)+{\mu }_2+\varphi {]}^{y_1+y_2+\varphi }}}+{\displaystyle \frac{\delta }{[({\mu }_1+{\mu }_2)(1+c)+\varphi {]}^{y_1+y_2+\varphi }}}\}\end{array}$$

(12)

and

$$\begin{array}{rl}p(y_1,y_2)& ={\displaystyle \frac{{\mu }_1^{y_1}{\mu }_2^{y_2}{\varphi }^{y_1+y_2}{e}^{\varphi }}{y_1!y_2!}}\sqrt{{\displaystyle \frac{2}{\pi }}}\{{\displaystyle \frac{[1+\delta e^{-(y_1+y_2)}]K_{y_1+y_2-1/2}({\upsilon }_1)}{{\upsilon }_1^{y_1+y_2-1/2}}}-{\displaystyle \frac{\delta e^{-y_1}{K}_{y_1+y_2-1/2}({\upsilon }_2)}{{\upsilon }_2^{y_1+y_2-1/2}}}\\ & -{\displaystyle \frac{\delta e^{-y_2}{K}_{y_1+y_2-1/2}({\upsilon }_3)}{{\upsilon }_3^{y_1+y_2-1/2}}}+{\displaystyle \frac{\delta K_{y_1+y_2-1/2}({\upsilon }_4)}{{\upsilon }_4^{y_1+y_2-1/2}}}\}\end{array}$$

(13)

where ${\upsilon }_1=\sqrt{{\varphi }^2+2({\mu }_1+{\mu }_2)\varphi }$ , ${\upsilon }_2=\sqrt{{\varphi }^2+2[{\mu }_1+{\mu }_2(1+c)]\varphi }$ , ${\upsilon }_3=\sqrt{{\varphi }^2+2[{\mu }_1(1+c)+{\mu }_2]\varphi }$ , ${\upsilon }_4=\sqrt{{\varphi }^2+2[({\mu }_1+{\mu }_2)(1+c)]\varphi }$ , and

$$K_{\lambda }\left(\omega \right)\equiv {\displaystyle \frac{1}{2}}{\int }_0^{\mathrm{\infty }}{u}^{\lambda -1}\exp \{-{\displaystyle \frac{\omega }{2}}(u+{\displaystyle \frac{1}{u}})\}\mathrm{d}u,\omega >0,\lambda \in \mathbb{R}$$

with $K_{\lambda }(\cdot )$ being the modified Bessel function of the third kind.

In the next section, we introduce the class of BMP regression models and discuss how to perform inference, including Monte Carlo simulation results for checking the finite-sample behavior of the proposed estimators and related standard error evaluation.

3. Regression modeling and estimation

In this section, we discuss the maximum likelihood method for inference of the parameters of the BMP regression models. For bivariate count random variables $\{(Y_{1i},Y_{2i}){\}}_{i=1}^n$ being BMP distributed, with observed values $\{(y_{1i},y_{2i}){\}}_{i=1}^n$ , let us consider two regression structures for the marginal means $\{({\mu }_{1i},{\mu }_{2i}){\}}_{i=1}^n$ . Hence, consider the following regression structures:

$$\begin{array}{rl}g({\mu }_{1i})& ={\mathit{x}}_i^{\mathrm{\top }}\mathit{\alpha }\end{array}$$

(14)

$$\begin{array}{rl}h({\mu }_{2i})& ={\mathit{v}}_i^{\mathrm{\top }}\mathit{\beta }\end{array}$$

(15)

where ${\mathit{x}}_i$ and ${\mathit{v}}_i$ are the explanatory variables vectors with $p\times 1$ and $q\times 1$ dimensions, respectively, for $i=1,\dots ,n$ , with $n$ denoting the sample size. Regarding the unknown regression coefficients, $\mathit{\alpha }=({\alpha }_1,\dots ,{\alpha }_p{)}^{\mathrm{\top }}$ and $\mathit{\beta }=({\beta }_1,\dots ,{\beta }_q{)}^{\mathrm{\top }}$ are the vectors of the parameters related to the covariates. One has varied alternatives available for the functions $g(\cdot )$ and $h(\cdot )$ that establish the functional relationship between the modeled parameters and the linear predictors. Therefore, it is meaningful to state that one must pick up, preferably, a link function that guarantees the support of the model parameters. Regular choices are the logarithmic (log) and the square root link functions, to estimate positive parameters, and the logit and the probit link functions, to estimate parameters with domain in the range $(0,1)$ , to name some examples.

Considering $\mathit{\theta }=({\mathit{\alpha }}^{\mathrm{\top }},{\mathit{\beta }}^{\mathrm{\top }},\delta ,\varphi )$ as the vector of the model parameters, the log-likelihood function is given by $\sum _{i=1}^n\log p(y_{1i},y_{2i})$ , with $p(\cdot ,\cdot )$ stated in (8). More specifically, the log-likelihood function is proportional to the following expression:

$$\begin{array}{rl}\ell \left(\mathit{\theta }\right)& \propto \sum _{i=1}^n[y_{1i}\log {\mu }_{1i}+y_{2i}\log {\mu }_{2i}-\varphi b({\xi }_0)+\log \{{\mathrm{\Psi }}_1(y_{1i}+y_{2i};{\mu }_{1i},{\mu }_{2i},\varphi )\\ & +\delta [e^{-(y_{1i}+y_{2i})}{\mathrm{\Psi }}_1(y_{1i}+y_{2i};{\mu }_{1i},{\mu }_{2i},\varphi )-e^{-y_{1i}}{\mathrm{\Psi }}_2(y_{1i}+y_{2i};{\mu }_{1i},{\mu }_{2i},\varphi )\\ & -e^{-y_{2i}}{\mathrm{\Psi }}_3(y_{1i}+y_{2i};{\mu }_{1i},{\mu }_{2i},\varphi )+{\mathrm{\Psi }}_4(y_{1i}+y_{2i};{\mu }_{1i},{\mu }_{2i},\varphi )]\}]\end{array}$$

(16)

The log-likelihood function is, then, maximized over $\mathit{\theta }$ , denoting $\hat{\mathit{\theta }}={\text{argmax}}_{\mathit{\theta }}\ell (\mathit{\theta })$ as the MLE of $\mathit{\theta }$ , which can be obtained numerically through some optimization algorithm. Furthermore, the standard errors of the MLEs are directly obtained from the Hessian matrix. To achieve this, we use the optim function from the stats package in the R program, ²⁷ implementing the Nelder-Mead algorithm. The output of the generalized additive models for location, scale, and shape (GAMLSS) models, from the gamlss package, is used to provide initial parameter estimates, which are then used to initialize the algorithm. Specifically, in gamlss, we choose either the negative binomial model or the Poisson inverse-Gaussian model to first fit the data, and the resulting output is used to initialize the algorithm for fitting the BNB or BPIG models, respectively. We also set the initial guess for $\delta$ as zero.

3.1. Monte Carlo simulation

A Monte Carlo study is performed to evaluate the performance of MLEs for finite samples. For this simulation analysis, we have set, in expressions (14) and (15), $g(\cdot )=h(\cdot )=\log (\cdot )$ and applied it to the regression models since the logarithmic link function is a typical choice. However, it should be noted that other link functions can be used with a preference for those that guarantee the support of the model parameters. Hence, 5000 Monte Carlo replicas were run, with the following structure for the marginal means

$$\begin{array}{rl}\log {\mu }_{1i}& ={\alpha }_0+{\alpha }_1{x}_{1i}+{\alpha }_2{x}_{2i}\\ \log {\mu }_{2i}& ={\beta }_0+{\beta }_1{x}_{1i}+{\beta }_2{x}_{2i}\end{array}$$

for $i=1,\dots ,n$ , where $x_{1i}$ was generated from the Bernoulli distribution with a 0.5 success probability, while $x_{2i}$ was raised from the uniform distribution with $(0,0.8)$ range. Additionally, the parameter vector $\mathit{\theta }=({\alpha }_0,{\alpha }_1,{\alpha }_2,{\beta }_0,{\beta }_1,{\beta }_2,\varphi ,\delta {)}^{\mathrm{\top }}=(-2.0,0.3,1.4,-2.0,0.7,3.0,1.0,0.6{)}^{\mathrm{\top }}$ was defined, considering two regressors for the response variables. We emphasize that the values of the parameter vector were derived from the dataset described in the illustration Section 5. Initially, the number of consultations with a doctor and the number of prescribed medications were chosen as response variables, as other authors had previously analyzed them in similar contexts. The response variables were subsequently adjusted based on the individual’s gender and age, with the resulting coefficients serving as true parameter values in the simulation scenario. Hence, to justify what was previously established, gender was simulated using a Bernoulli distribution, with the success parameter representing the proportion of female individuals, which was approximately balanced with male individuals. Age was simulated using a uniform distribution, considering that, in the dataset, this variable is scaled by a factor of 100. We take into account samples with sizes $n=50,100,200,500,1000$ . In the Monte Carlo simulation, the values of all regressors were kept constant. Two particular cases of the BMP models were considered for the simulation study, the BNB and BPIG regression models.

We start the analysis of the simulation results from Table 2, which comprises the empirical mean and the root mean square error (RMSE) of the parameter estimates for the BNB regression model. We observe small bias and RMSE across most sample size configurations, except for sample sizes 50 and 100. Notably, the parameters $\delta$ and $\varphi$ exhibit the largest bias and RMSE in these smaller sample sizes. However, both bias and RMSE decrease as the sample size increases. Thus, we conclude that the MLEs perform well in moderate finite samples.

Table 2.

Empirical mean and root mean square error (in parentheses) of the estimates for the BNB regression model, with $n=50,100,200,500,1000$ .

BNB	$n$ = 50	$n$ = 100	$n$ = 200	$n$ = 500	$n$ = 1000
${\alpha }_0$	$-$ 2.318	$-$ 2.108	$-$ 2.046	$-$ 2.024	$-$ 2.014
	(1.938)	(0.841)	(0.582)	(0.349)	(0.240)
${\alpha }_1$	0.407	0.330	0.322	0.299	0.305
	(1.697)	(0.657)	(0.460)	(0.275)	(0.192)
${\alpha }_2$	1.599	1.467	1.404	1.417	1.411
	(2.104)	(1.496)	(0.980)	(0.618)	(0.421)
${\beta }_0$	$-$ 2.141	$-$ 2.071	$-$ 2.042	$-$ 2.017	$-$ 2.011
	(1.003)	(0.662)	(0.466)	(0.281)	(0.193)
${\beta }_1$	0.736	0.723	0.718	0.701	0.705
	(0.708)	(0.478)	(0.340)	(0.205)	(0.145)
${\beta }_2$	3.116	3.052	3.027	3.016	3.012
	(1.525)	(1.132)	(0.757)	(0.478)	(0.324)
$\varphi$	1.680	1.293	1.130	1.045	1.021
	(3.500)	(1.265)	(0.510)	(0.265)	(0.175)
$\delta$	2.031	0.941	0.762	0.651	0.606
	(13.559)	(2.830)	(1.621)	(0.999)	(0.719)

BNB: bivariate negative binomial.

The suitable performance of the estimation method based on the simulation study is reported in Figures 3 and 5 in the Supplemental Material, which covers the histograms of the standardized (mean-standard deviation) parameter estimates, seeming to be asymptotically normally distributed with the increase of the sample size, as they follow along with normal density curves.

Figure 1.

Descriptive data analysis on the number of doctor consultations versus the number of non-prescribed medications, from a sample of the Australian Health Survey (AHS) study, 1977–1978.

Figure 2.

Histograms of the randomized PIT under the BNB regression for the number of doctor consultations and the number of non-prescribed medications. PIT: probability integral transform; BNB: bivariate negative binomial.

Figure 3.

Histograms of the randomized PIT under the BPIG regression for the number of doctor consultations and the number of non-prescribed medications. PIT: probability integral transform; BPIG: bivariate Poisson-inverse Gaussian.

Analyzing the simulation study for the BPIG regression model, through Table 3, we anticipate results similar to that of the BNB regression model. Therefore, small bias and RMSE for almost all configurations of considered sample sizes and asymptotic normal distribution of the parameter estimates, as the sample size increases, except the scenarios with sample sizes 50 and 100. The histograms of the standardized parameter estimates, along with a normal density curve, have a pattern similar to those produced by the BNB regression model, as presented in Figures 4 and 6 from the Supplemental Material.

Figure 4.

Simulated envelopes for Pearson residuals under the bivariate negative binomial (BNB) regression model for the number of doctor consultations and the number of non-prescribed medications.

Figure 5.

Simulated envelopes for Pearson residuals under the bivariate Poisson-inverse Gaussian (BPIG) regression model for the number of doctor consultations and the number of non-prescribed medications.

Figure 6.

Histograms of the randomized PIT under the BP regression for the number of doctor consultations and the number of non-prescribed medications. PIT: probability integral transform; BP: bivariate Poisson.

Table 3.

Empirical mean and root mean square error (in parentheses) of the estimates for the BPIG regression model, with $n=50,100,200,500,1000$ .

BPIG	$n$ = 50	$n$ = 100	$n$ = 200	$n$ = 500	$n$ = 1000
${\alpha }_0$	$-$ 2.239	$-$ 2.097	$-$ 2.040	$-$ 2.023	$-$ 2.012
	(1.839)	(0.843)	(0.589)	(0.350)	(0.236)
${\alpha }_1$	0.331	0.314	0.308	0.303	0.301
	(1.523)	(0.661)	(0.458)	(0.277)	(0.193)
${\alpha }_2$	1.517	1.427	1.415	1.412	1.408
	(2.136)	(1.524)	(0.991)	(0.626)	(0.418)
${\beta }_0$	$-$ 2.130	$-$ 2.069	$-$ 2.028	$-$ 2.017	$-$ 2.009
	(0.976)	(0.660)	(0.461)	(0.285)	(0.191)
${\beta }_1$	0.718	0.715	0.705	0.706	0.700
	(0.692)	(0.475)	(0.336)	(0.205)	(0.142)
${\beta }_2$	3.092	3.046	3.016	3.012	3.009
	(1.495)	(1.141)	(0.749)	(0.476)	(0.318)
$\varphi$	1.853	1.378	1.178	1.061	1.032
	(3.191)	(1.222)	(0.712)	(0.341)	(0.225)
$\delta$	1.369	0.888	0.696	0.639	0.612
	(7.631)	(2.373)	(1.540)	(0.953)	(0.682)

BPIG: bivariate Poisson-inverse Gaussian.

We now present in Tables 4 and 5 the empirical standard error of the estimates stated in the Monte Carlo simulation and the average of the standard errors of these estimates obtained via the Hessian matrix. From the results of these tables, we conclude that the outcomes of the Hessian matrix are consistent for estimating the standard errors of the parameter estimates, for both the BNB and the BPIG regression models, except for sample size 50.

Table 4.

Empirical and Hessian derived standard errors of the parameter estimates for the BNB regression model, with $n=50,100,200,500,1000$ .

	${\alpha }_0$	${\alpha }_1$	${\alpha }_2$	${\beta }_0$	${\beta }_1$	${\beta }_2$	$\varphi$	$\delta$
$n=50$
Empirical	1.352	1.198	1.481	0.702	0.500	1.075	2.428	9.534
Hessian	3.436	3.256	1.342	0.625	0.464	0.979	1.494	3.038
$n=100$
Empirical	0.590	0.464	1.057	0.466	0.338	0.800	0.870	1.987
Hessian	0.556	0.448	1.031	0.445	0.327	0.779	0.649	1.682
$n=200$
Empirical	0.410	0.325	0.693	0.328	0.240	0.535	0.349	1.141
Hessian	0.397	0.314	0.684	0.318	0.233	0.523	0.328	1.143
$n=500$
Empirical	0.246	0.194	0.437	0.198	0.145	0.338	0.185	0.705
Hessian	0.241	0.193	0.432	0.197	0.143	0.332	0.179	0.725
$n=1000$
Empirical	0.170	0.136	0.298	0.136	0.103	0.229	0.123	0.509
Hessian	0.165	0.136	0.294	0.136	0.101	0.227	0.121	0.515

BNB: bivariate negative binomial.

Table 5.

Empirical and Hessian derived standard errors of the parameter estimates for the BPIG regression model, with $n=50,100,200,500,1000$ .

	${\alpha }_0$	${\alpha }_1$	${\alpha }_2$	${\beta }_0$	${\beta }_1$	${\beta }_2$	$\varphi$	$\delta$
$n=50$
Empirical	1.290	1.077	1.508	0.684	0.489	1.055	2.174	5.368
Hessian	15.06	14.88	1.335	0.613	0.453	0.955	1.802	2.612
$n=100$
Empirical	0.592	0.467	1.077	0.464	0.336	0.806	0.822	1.665
Hessian	0.554	0.447	1.029	0.440	0.322	0.765	0.805	1.566
$n=200$
Empirical	0.415	0.323	0.701	0.326	0.237	0.529	0.487	1.087
Hessian	0.395	0.312	0.681	0.314	0.230	0.515	0.436	1.076
$n=500$
Empirical	0.247	0.196	0.442	0.201	0.145	0.336	0.237	0.673
Hessian	0.241	0.192	0.430	0.195	0.141	0.328	0.229	0.682
$n=1000$
Empirical	0.166	0.137	0.295	0.135	0.100	0.225	0.158	0.482
Hessian	0.165	0.135	0.293	0.135	0.100	0.224	0.155	0.484

BPIG: bivariate Poisson-inverse Gaussian.

Remark 3.1

More simulated results are provided in the Supplemental Material (Section 2) about the adverse effects of incorrectly ignoring conditional dependence, that is assuming $\delta =0$ when in fact $\delta \ne 0$ .

4. Hypothesis testing, residual, and model diagnostics

The model specification encompasses the following stages: estimation, diagnostic checking, and inference. This section discusses some techniques for performing such tasks. In particular, we are interested in investigating possible sources of correlation. To this end, we conduct hypothesis testing involving the correlation parameter $\delta$ using the likelihood ratio test (LRT) which is discussed below. For model evaluation, residual analysis will be conducted based on simulated envelopes and randomized probability integral transformation as a diagnostic tool.

4.1. LRT and simulated results

Regarding some hypotheses on the BMP regression models, we are interested in determining when the parameter $\delta$ , associated with the correlation structure, is significantly non-zero. Hence, we have the null and the alternative hypotheses expressed as $H_0:\delta =0$ and $H_1:\delta \ne 0$ , respectively. Under the LRT, denote $\mathbf{\Omega }$ to be the parameter space of the model and $\ell ({\mathit{\theta }}_0)$ as the log-likelihood function assuming the null hypothesis to be true, for ${\mathit{\theta }}_0\in {\mathbf{\Omega }}_0$ , with ${\mathbf{\Omega }}_0\subset \mathbf{\Omega }$ , and $\ell (\hat{\mathit{\theta }})$ as the log-likelihood function assuming the null hypothesis to be not true, with $\hat{\mathit{\theta }}\in \mathbf{\Omega }$ and where $\hat{\mathit{\theta }}={\text{argmax}}_{\mathit{\theta }}\ell (\mathit{\theta })$ , the LRT statistic is expressed as follows:

$$\mathrm{LRT}(\mathit{\theta })=-2[\ell ({\mathit{\theta }}_0)-\ell (\hat{\mathit{\theta }})]$$

which is provided by Wilks ²⁸ to be asymptotically chi-squared distributed with degrees of freedom established by the difference of the dimensions of $\mathbf{\Omega }$ and ${\mathbf{\Omega }}_0$ . The finite-sample performance of the LRT is evaluated via a simulation study. We present the relative frequency of type-I error of this simulation study, considering nominal significance levels of 0.01, 0.05, and 0.10 under both BNB and BPIG regression models. We generate data from these models with the same settings considered in the Monte Carlo study provided in Subsection 3.1 with an exception for $\delta$ ; here we set $\delta =0$ . Tables 6 and 7 present the empirical type-I errors for the LRT under the BNB and BPIG regression models. From these results, we can observe that the LRT is yielding the desirable nominal level since the relative frequencies are close to the established significance levels.

Table 6.

Empirical type-I errors of the LRT under the BNB regression model, with $n=200,500,1000$ .

	Type-I error
Sample size	0.01	0.05	0.10
$n=$ 200	0.0120	0.0546	0.1002
$n=$ 500	0.0096	0.0474	0.0948
$n=$ 1000	0.0092	0.0494	0.0998

LRT: likelihood ratio test; BNB: bivariate negative binomial.

Table 7.

Empirical type-I errors of the LRT under the BPIG regression model, with $n=200,500,1000$ .

	Type-I error
Sample size	0.01	0.05	0.10
$n=$ 200	0.0102	0.0490	0.1006
$n=$ 500	0.0084	0.0494	0.0888
$n=$ 1000	0.0088	0.0458	0.0904

LRT: likelihood ratio test; BPIG: bivariate Poisson-inverse Gaussian.

4.2. Residual and diagnostic methods

Considering the model evaluation stage, one might conduct residual analysis and use goodness-of-fit measures. According to Cameron and Trivedi, ⁷ one carries the residual analysis for many purposes, such as to detect model misspecification, outliers, poor fit, and influential observations. Therefore, techniques to perform residual analysis will measure the departure between the fitted and the original values of the dependent variable, indicating that it is essential. In addition, a visual inspection may potentially indicate the nature of model misspecification and the magnitude of its effect. The literature has proposed miscellaneous residuals for count data, and one of the nominees presented by Cameron and Trivedi ⁷ is the Pearson residual, also known as the standardized residual. In order to check the goodness-of-fit of the joint model instead of a marginal approach, we consider the Pearson residual of the sum of the count variables as follows:

$$r_i={\displaystyle \frac{y_{1i}+y_{2i}-{\hat{\mu }}_{1i}-{\hat{\mu }}_{2i}}{\sqrt{{\hat{\sigma }}_{1i}^2+{\hat{\sigma }}_{2i}^2+2\widehat{\text{cov}}(y_{1i},y_{2i})}}},i=1,\dots ,n$$

(17)

where

$$\begin{array}{rl}& {\hat{\mu }}_{1i}=g^{-1}\left({\mathit{x}}_i^{\mathrm{\top }}\hat{\mathit{\alpha }}\right)\\ & {\hat{\mu }}_{2i}=h^{-1}\left({\mathit{v}}_i^{\mathrm{\top }}\hat{\mathit{\beta }}\right)\\ & {\hat{\sigma }}_{1i}^2=g^{-1}\left({\mathit{x}}_i^{\mathrm{\top }}\hat{\mathit{\alpha }}\right)\{1+g^{-1}\left({\mathit{x}}_i^{\mathrm{\top }}\hat{\mathit{\alpha }}\right){\displaystyle \frac{b^{″}({\xi }_0)}{\hat{\varphi }}}\}\\ & {\hat{\sigma }}_{2i}^2=h^{-1}\left({\mathit{v}}_i^{\mathrm{\top }}\hat{\mathit{\beta }}\right)\{1+h^{-1}\left({\mathit{v}}_i^{\mathrm{\top }}\hat{\mathit{\beta }}\right){\displaystyle \frac{b^{″}({\xi }_0)}{\hat{\varphi }}}\},\text{and}\\ & \widehat{\text{cov}}(y_{1i},y_{2i})=g^{-1}\left({\mathit{x}}_i^{\mathrm{\top }}\hat{\mathit{\alpha }}\right)h^{-1}\left({\mathit{v}}_i^{\mathrm{\top }}\hat{\mathit{\beta }}\right)\left[{\displaystyle \frac{(1+\hat{\delta }d{c}^2)(1+\hat{\varphi })-\hat{\varphi }}{\hat{\varphi }}}\right]& \end{array}$$

where $d=\exp \{-[g^{-1}({\mathit{x}}_i^{\mathrm{\top }}\hat{\mathit{\alpha }})+h^{-1}({\mathit{v}}_i^{\mathrm{\top }}\hat{\mathit{\beta }})]c\}$ , with $\hat{\mathit{\alpha }}$ , $\hat{\mathit{\beta }}$ and $\hat{\varphi }$ being the MLEs of $\mathit{\alpha }$ , $\mathit{\beta }$ , and $\varphi$ , respectively.

An ordinary way to employ residuals is to plot them against the normal quantiles. However, even though Pearson’s residuals have zero mean and unit variance for large samples, they are skewed in distribution. Therefore, we expect the normal approximation to be poor, even for moderate sample sizes. A path to overpower this barrier is to construct simulated envelopes for the residuals, as suggested by Atkinson, ²⁹ and Hinde and Demétrio. ³⁰ The steps to produce simulated envelopes for count regression models follow the description of Algorithm 2. Here, we will demonstrate the performance of these simulated envelopes in the empirical illustration in Section 5.

Regarding the calibration and sharpness of the model, Czado et al. ³¹ pointed out that a calibrated model refers to the consistency between the observed variable and its predictions. Hence, one can request the probability integral transform (PIT)^32,33 to check if it is possible to conclude that the observations of a continuous variable are derived from the suggested predictive distribution. For this, we can proceed graphically or through hypothesis tests since, for suitable cases, the PIT must be uniform in distribution.

However, when it comes to a discrete variable, it is well-known that the PIT does not have a uniform distribution. Consequently, the literature introduced the randomized PIT,^34–36 which follows the standard uniform distribution, to make it available for count data. Therefore, we suggest using the randomized probability integral transform histogram for calibration model control as a diagnostic tool:

$$\mathrm{PIT}(Y_k)=U{F}_{Y_k}(y)+(1-U)F_{Y_k}(y-1),k=1,2$$

(18)

where $Y_k$ is the count random variable with marginal distribution $F_{Y_k}(\cdot )$ , assuming $F_{Y_k}(-1)=0$ , and $U$ is the standard uniform distributed.

Residual analysis and diagnostics were conducted marginally to smooth the model evaluation process. The bivariate residual analysis is not straightforward since the idea of sample ordering is lost, as pointed out by Moral et al. ³⁷ Moreover, regarding model calibration, it is not trivial to determine if the observations can be said to be derived from the suggested predictive distribution, not to mention that the bivariate PIT is no longer uniformly distributed. For a thorough assessment of the model, which will be further developed in future research for the class of BMP regression models, we refer to the paper by Genest and Rivest, ³⁸ which proposes a bivariate probability integral transformation based on copulas and makes considerations for a multivariate extension. Also, we refer to the work by Moral et al., ³⁷ in which the authors propose a generalization of the simulated envelopes yielding simulated polygons for bivariate residual analysis. We illustrate the usefulness of the proposed forthright residual analysis and randomized probability integral transformation in the real data analysis in Section 5.

5. Analysis of doctor consultations and non-prescribed drugs intake

The proposed class of models was motivated by the problem of finding the association between the number of consultations with a doctor and consumption of non-prescription medication. The dataset regards information on health conditions and, as well, demographic and socioeconomic characteristics of the individuals, collected from the 1977–1978 AHS, conducted by the Australian Bureau of Statistics (ABS). Since 1989, the Australian health survey series has been part of the National Health Survey (NHS), with similar information, and one can grasp more about the NHS at https://www.abs.gov.au/statistics/health.

The employed dataset, a sample of the 1977–1978 AHS, was first analyzed by Cameron et al. ³⁹ through a microeconometric model for healthcare demand. This dataset is available in the book by Cameron and Trivedi ⁷ and, later, in the book by Faraway, ⁴⁰ which made it ready for use on the faraway package in the CRAN of $\mathtt{R}$ software. The AHS is a comprehensive source of healthcare data, and several researchers have used this dataset to illustrate their models for bivariate count analysis. For example, the AHS sample used here is the same dataset used by Gurmu and Elder, ¹¹ Karlis and Ntzoufras, ⁴¹ and Famoye, ¹³ where these authors jointly analyzed the number of consultations with a doctor and the number of consultations with non-doctor health professionals or the number of prescribed medications.

To investigate the association between the count of doctor consultations and the non-prescribed drug intake, we apply the proposed model on the sample of the 1977–1978 AHS which contains attributes of 5190 individuals aged 15 years or older. Unlike previous research, we considered, as count responses, variables that are negatively correlated—the number of consultations with a doctor and the number of non-prescribed medications. The explanatory variables in this dataset are enumerated as follows age, ¹ gender (0 = male, 1 = female), income, ² insurance coverage (0 = no health insurance, 1 = covered by public/private health insurance), number of illnesses, number of days of reduced activity, health score (high score indicates bad health), and chronic conditions (0 = no chronic conditions, 1 = chronic condition(s)).

Regarding covariates, there is more than one type of health insurance, such as private and government coverage, and, as well, there is more than one combination of the chronic conditions variable. Hence, for a straightforward interpretation, we decided to dichotomize these covariates. Moreover, for the number of illnesses covariate, it can be seen ³ that, for individuals with five or more illnesses, it was coded as 5, so it is more appropriate to treat this variable in categories: Illness1 (1 = 1 or 2 illnesses, 0 = otherwise); Illness2 (1 = 3 or 4 illnesses, 0 = otherwise); Illness3 (1 = 5 illnesses or more, 0 = otherwise); No illnesses (baseline).

Figure 1 shows the number of non-prescribed medications versus the number of consultations with a doctor, with the size of the points in this plot being proportional to the frequency of these paired count variables. In a first look at the plot, we notice small counts for both variables, with a higher proportion of zero doctor consultations and, as well, zero non-prescribed medications. Also, as anticipated, a lower number of non-prescribed medications would be expected for an individual with a higher count of consultations with a doctor, as this patient will be more likely to purchase prescribed medications. Conversely, a higher number of non-prescribed medications is associated with fewer doctor appointments. Again, this was expected since these individuals, without medical advice, will acquire non-prescribed medications. Thus, it seems that for this dataset, as we can see through Figure 1, an inverse relationship between the number of non-prescribed medications and the number of consultations with a doctor, indicating the existence of a negative correlation.

To fit statistical models, we consider the BNB and the BPIG regressions, which belong to our class of BMP models. The results of the BNB regression model are reported in Table 8 which presents the parameter estimates, the standard error estimates, z-values, and associated p-values. After initial data analysis, the model fit chose as significant covariates gender, age, income, number of illnesses, number of days of reduced activity, and health score. The results of this modeling lead us to conclude that all covariates are significant, considering a significance level of 5%. Also, through the LRT for the parameter $\delta$ , we can notice statistical evidence, regarding the standard 5% level of significance, to reject the null hypothesis that this parameter equals zero, as the associated p-value is less than the 0.1% level of significance.

With respect to the BPIG regression model, the preliminary analysis selected the same covariates as in the fit of the BNB regression model, with the addition of the explanatory variable chronic conditions. Table 9 reveals that all covariates are significant, considering the significance level of 5%. Similarly, the LRT for the parameter $\delta$ indicates the statistical evidence prone to a non-nullity of this parameter, as the associated p-value is less than the 0.1% level of significance.

Table 8.

Parameter estimates, standard errors, z-values, and p-values for the BNB regression model, applied to the number of doctor consultations and the number of non-prescribed medications.

	Est.	S.E.	z-value	p-value
Intercept	$-$ 3.074	0.113	$-$ 27.146	<0.001
Gender	0.148	0.061	2.444	0.015
Age	0.919	0.146	6.285	<0.001
Illness1	1.218	0.104	11.724	<0.001
Illness2	1.335	0.116	11.509	<0.001
Illness3	1.607	0.140	11.444	<0.001
R.actdays	0.117	0.006	18.813	<0.001
Hscore	0.562	0.116	4.848	<0.001
Intercept	$-$ 1.901	0.104	$-$ 18.359	<0.001
Gender	0.172	0.055	3.133	0.002
Age	$-$ 0.435	0.139	$-$ 3.121	0.002
Income	0.334	0.073	4.586	<0.001
Illness1	0.765	0.073	10.432	<0.001
Illness2	0.962	0.089	10.785	<0.001
Illness3	1.386	0.120	11.556	<0.001
Hscore	0.419	0.113	3.723	<0.001
$\delta$	$-$ 2.034	0.148	$-$ 13.734	<0.001
$\varphi$	1.820	0.177	$(-)$	$(-)$
LRT $(\delta )$	139.33			<0.001
AIC	14347.90

BNB: bivariate negative binomial; Est.: estimate; S.E.: standard error; LRT: likelihood ratio test; AIC: Akaike information criterion.

Table 9.

Parameter estimates, standard errors, z-values, and p-values for the BPIG regression model, applied to the number of doctor consultations and the number of non-prescribed medications.

	Est.	S.E.	z-value	p-value
Intercept	$-$ 3.206	0.120	$-$ 26.803	<0.001
Gender	0.221	0.063	3.498	<0.001
Age	0.861	0.151	5.709	<0.001
Illness1	1.298	0.109	11.860	<0.001
Illness2	1.510	0.121	12.454	<0.001
Illness3	1.629	0.148	11.042	<0.001
R.actdays	0.125	0.006	19.474	<0.001
Hscore	0.488	0.121	4.037	<0.001
Intercept	$-$ 1.978	0.106	$-$ 18.634	<0.001
Gender	0.232	0.057	4.094	<0.001
Age	$-$ 0.512	0.149	$-$ 3.429	0.001
Income	0.396	0.075	5.317	<0.001
Illness1	0.749	0.075	9.946	<0.001
Illness2	0.889	0.095	9.386	<0.001
Illness3	1.196	0.128	9.364	<0.001
Hscore	0.405	0.117	3.471	0.001
Chcond	0.147	0.061	2.422	0.015
$\delta$	$-$ 2.500	0.177	$-$ 14.125	<0.001
$\varphi$	1.298	0.117	$(-)$	$(-)$
LRT $(\delta )$	150.44			<0.001
AIC	14324.86

BPIG: bivariate Poisson-inverse Gaussian; Est.: estimate; S.E.: standard error; LRT: likelihood ratio test; AIC: Akaike information criterion.

To check potential flaws in the model and prediction, we first examined the BNB regression model. The randomized PIT histogram in Figure 2 appears to be approximately uniformly distributed, inducing a consistency for both responses, the number of consultations with a doctor, and the number of non-prescribed medications. Likewise, the BPIG regression model also suggests a consistency for the number of doctor consultations and the number of non-prescribed medications, as we observe calibration in Figure 3 through the histograms that appear to be nearly uniformly distributed.

Continuing the modeling cycle, we are now interested in verifying whether the assumed BMP distribution is adequate for the dataset considered here. Figures 4 and 5 deliver the simulated envelopes for the Pearson residual against the theoretical quantiles of the standard normal distribution, respectively, for the BNB and the BPIG regression models. Assuming the unobserved heterogeneity is either gamma-distributed or inverse-Gaussian distributed, both choices look suitable since almost all residuals remain within the simulated envelopes.

Note that around 98.6% of the residuals within the simulated envelopes for the BNB regression model, and almost 99.5% for the BPIG regression model, with even better performance for the latter.

Not considering overdispersion in the modeling process can compromise the data analysis. Thus, the BP regression model was applied to underline the need for a model that allows overdispersion; the summary fit is provided in Table 10. Regarding its randomized PIT, Figure 6 exposes a slight deviation from uniformity, implying a reduced calibration compared to the BMP regression models. Even more worrisome, Figure 7 points out that the BP distribution is not indicated for this dataset since only 56.1% of the residuals lie inside the simulated envelopes.

Figure 7.

Simulated envelopes for Pearson residuals under the bivariate Poisson (BP) regression model for the number of doctor consultations and the number of non-prescribed medications.

Table 10.

Parameter estimates, standard errors, z-values, and p-values for the BP regression model, applied to the number of doctor consultations and the number of non-prescribed medications.

	Est.	S.E.	z-value	p-value
Intercept	$-$ 2.867	0.102	$-$ 27.991	<0.001
Gender	0.147	0.054	2.733	0.006
Age	0.742	0.130	5.718	<0.001
Illness1	1.088	0.095	11.447	<0.001
Illness2	1.204	0.105	11.461	<0.001
Illness3	1.419	0.121	11.714	<0.001
R.actdays	0.121	0.005	24.811	<0.001
Hscore	0.397	0.097	4.099	<0.001
Intercept	$-$ 1.744	0.094	$-$ 18.638	<0.001
Gender	0.215	0.050	4.304	<0.001
Age	$-$ 0.653	0.131	$-$ 4.979	<0.001
Income	0.288	0.066	4.345	<0.001
Illness1	0.659	0.068	9.738	<0.001
Illness2	0.759	0.084	8.995	<0.001
Illness3	1.142	0.109	10.484	<0.001
Hscore	0.385	0.099	3.898	<0.001
Chcond	0.123	0.054	2.288	0.022
$\delta$	$-$ 0.472	0.162
LRT $(\delta )$	12.39			0.0004
AIC	14689.74

BP: bivariate Poisson; Est.: estimate; S.E.: standard error; LRT: likelihood ratio test; AIC: Akaike information criterion.

We conclude this section on the interpretation of the model, using the relativities provided in Table 11. These measures aim to compare a covariate’s category with its baseline in terms of predicted response.

Table 11.

Relativities of the explanatory variables, to estimate the number of doctor consultations and non-prescribed medications, under the BPIG regression model.

Explanatory	Doctor	Non-prescribed
variables	consultations	medications
Intercept	0.040	0.138
Female	1.247	1.261
Age	1.009 ${}^{\mathrm{a}}$	0.995 ${}^{\mathrm{b}}$
Income	$(-)$	1.486 ${}^{\mathrm{c}}$
Ilness1	3.664	2.115
Ilness2	4.525	2.432
Ilness3	5.098	3.306
R.actdays	1.133	1.499
Hscore	1.630	1.159

BPIG: bivariate Poisson-inverse Gaussian.

${}^{\mathrm{a}}$ Considering a one-year increase in life.

${}^{\mathrm{b}}$ Considering a one-year increase in life.

${}^{\mathrm{c}}$ Considering a $1000 increase in annual income.

To proceed with the interpretation of the model, we chose the BPIG regression as the final model for this empirical illustration since it presents the best performance due to the conjuncture of its results. Table 11 reveals that the expected number of doctor consultations and non-prescribed medications is about 25% and 26%, respectively, higher for females than males. For each one-year increase in age, one expects an almost 1% increase in the number of doctor consultations. Conversely, one expects a reduction of around 0.5% for the number of non-prescribed medications. For individuals with one or two illnesses, the expected number of doctor consultations and non-prescribed medications is nearly three times and twice, respectively, compared with those without illnesses. The interpretation of the relativities for other covariates follows similarly.

6. Concluding remarks and future research

We have introduced a new class of regression models to analyze bivariate count data. We have considered a latent component, with distribution belonging to the exponential family, to account for the unobserved heterogeneity inherent to such data. We proposed the maximum likelihood method for inference on the model parameters. We have evaluated the MLEs for two particular cases of our class, the BNB, and the BPIG models, concluding with a good finite-sample performance of the proposed estimators in both models. We also provided a general measure of model calibration and simulated envelopes for checking the model adequacy of our BMP regression models.

A complete statistical analysis of a random sample dataset of the 1977–1978 AHS from the Australian Bureau of Statistics has been conducted. We studied the number of doctor consultations by individuals and the number of non-prescribed medications acquired by these individuals, using these variables for joint modeling through the proposed BMP models. The bivariate analysis of this AHS count dataset is, to the best of our knowledge, a new contribution once we deliver a complete analysis of the data, showing that the BMP models provide an adequate joint fit to the number of doctor consultations and non-prescribed medications. We also illustrated that ignoring the unobserved heterogeneity in such data can lead to a poor model fit.

Towards future research, noteworthy issues that deserve further investigation are (i) the generalization of the model allowing for a varying dispersion parameter and a regression framework for the $\delta$ parameter, for an even more flexible correlation structure; (ii) the introduction of a zero-inflated version of the BMP models; and (iii) to design an R package for this class of BMP models.

Supplemental Material

sj-zip-1-smm-10.1177_09622802251345332 - Supplemental material for Paired count regressions for modeling the number of doctor consultations and non-prescribed drugs intake

Supplemental material, sj-zip-1-smm-10.1177_09622802251345332 for Paired count regressions for modeling the number of doctor consultations and non-prescribed drugs intake by Jussiane Nader Gonçalves,Wagner Barreto-Souza and Hernando Ombao in Medical Research

Supplemental Material

sj-pdf-2-smm-10.1177_09622802251345332 - Supplemental material for Paired count regressions for modeling the number of doctor consultations and non-prescribed drugs intake

Supplemental material, sj-pdf-2-smm-10.1177_09622802251345332 for Paired count regressions for modeling the number of doctor consultations and non-prescribed drugs intake by Jussiane Nader Gonçalves,Wagner Barreto-Souza and Hernando Ombao in Medical Research