Analyzing dependence in incidence of diabetes and heart problem using generalized bivariate geometric models with covariates

ABSTRACT

For analyzing incidence data on diabetes and health problems, the bivariate geometric probability distribution is a natural choice but remained unexplored largely due to lack of models linking covariates with the probabilities of bivariate incidence of correlated outcomes. In this paper, bivariate geometric models are proposed for two correlated incidence outcomes. The extended generalized linear models are developed to take into account covariate dependence of the bivariate probabilities of correlated incidence outcomes for diabetes and heart diseases for the elderly population. The estimation and test procedures are illustrated using the Health and Retirement Study data. Two models are shown in this paper, one based on conditional-marginal approach and the other one based on the joint probability distribution with an association parameter. The joint model with association parameter appears to be a very good choice for analyzing the covariate dependence of the joint incidence of diabetes and heart diseases. Bootstrapping is performed to measure the accuracy of estimates and the results indicate very small bias.

KEYWORDS:

• Bivariate geometric model
• generalized linear models
• dependence in outcomes
• incidence of diabetes and heart diseases
• conditional model
• marginal model
• joint model

1. Introduction

In health data, often there are occurrences of two correlated outcomes for the first time or incidence of two symptoms, diseases or conditions. The repeated measures data on the outcomes at different times observed longitudinally provide us with the scope for modeling correlated outcomes as functions of potential risk factors. There is association between incidence of diabetes and heart diseases. However, we are yet to know the relationship between the incidence of diabetes mellitus and heart diseases fully. It is known that there is an increased risk of diabetes mellitus among the elderly people. Manson et al. [18] examined the relationship between diabetes mellitus and incidence of heart diseases (coronary heart disease, stroke, total cardiovascular diseases) using a cohort study data for women 30–55 years of age. They found the association between diabetes and increased risk of nonfatal myocardial infarction and fatal coronary heart disease. Similarly, according to Haffner et al. [13] the non-insulin dependent type 2 diabetes is associated with a marked increase in the risk of coronary heart disease. There is a very important issue of concern (debate) regarding whether the delayed onset of fatal and nonfatal coronary heart disease in women as compared to that of the men is offset (eradicated) attributable to the reduced risk as a result of estrogen therapy in postmenopausal period among the diabetic women [5–7]. Another possible link can be examined between the obesity, insulin resistance, and heart disease. Caballero [9] observed that there is no fully established cause–effect relationship but we can come closer to more conclusive evidence if combined together. In other words, still there is ambiguity in understanding the underlying mechanism properly. One of the reasons for the lack of understanding in such relationships is the dependence in outcomes as well as the dependence between risk factors and outcomes in the presence of dependence in outcomes. In most of the models, these diseases are not considered as joint outcomes with dependence and hence the relationships under such dependence between risk factors and outcomes cannot be assessed properly.

In the Health and Retirement Study [23], there are questions regarding occurrence of diabetes and heart problem during successive follow-ups or waves. These outcome variables are expected to be correlated. At this backdrop, the joint model for incidence of diabetes and heart problem considering the association between the outcome variables as well as the relationship between risk factors and outcome variables are of great interest to epidemiologists and biostatisticians as well as to the medical professionals in general. The analysis based on the dependence between diabetes and heart disease at the time of incidence (first occurrence) for each along with their dependence on potential risk factors may reveal the nature of realistic underlying relationships which may not be found from analysis based on univariate models for incidence of diabetes and heart problem separately. The conditional models alone cannot provide any realistic explanation to these relationships in true sense due to complexity associated between two outcome variables.

Dependence in outcome variables may arise in various fields such as epidemiology, time series, environment, public health, economics, and anthropology. In longitudinal studies, we deal with outcomes from more than one time point on each subject under study. If the observations are recorded on the same subjects at different times or observations recorded at the same time by different persons then there may be dependence in the outcome variables. We can use marginal, conditional and joint models to address the dependence in outcomes with covariate dependence. In this context, the occurrence of two events for the same individuals for the first time is of interest in many situations. As the dependence in two outcomes results in a bivariate geometric distribution, an attempt is made to propose a generalized linear model based on bivariate geometric distribution. A bivariate Bernoulli generalized linear model has been proposed by Islam et al. [15]. It is known that the geometric distribution is a discrete counterpart of the exponential distribution. In many instances, this can be applied to survival or reliability analysis for discrete time analysis for follow up or panel data where time to event data are not available but occurrence or non-occurrence of events are reported at different times. Examples of the application of reliability and survival analysis concepts in discrete time are discussed in [12,20,27]. Nair and Nair [21] provided a form of the bivariate geometric distribution which is employed in this paper. Other forms can be seen in [2,14,22,25]. Analogous to bivariate exponential model by Marshal and Olkin [19], Basu and Dhar [8] proposed a bivariate geometric model. Characterization results are shown by Sun and Basu [26] and Sreehari [24] and Sreehari and Vasudeva [25]. Recently characterizations of bivariate geometric distributions were also shown by Jayakumar and Mundassery [16] and Li and Dhar [17]. An application of bivariate lifetime geometric distribution in the presence of cure fraction was studied by Davarzani et al. [10]. Arnold [1] and Arnold and Strauss [3,4] proposed bivariate distributions with exponential conditionals rather than using the joint distributions directly in order to make the modeling easier. In this paper, we have proposed two generalized linear models, for bivariate geometric and the link functions are identified. The first model is based on the marginal and conditional distributions to obtain joint distribution [22,25]. The second model is based on a joint distribution form of bivariate geometric proposed by Nair and Nair [21]. It is noteworthy that there have been no attempts so far to develop generalized linear models based on the bivariate geometric forms. In this paper, the link functions are obtained by extending the generalized linear model for the bivariate geometric model using the bivariate exponential form. Similar extended approach with theoretical justification was illustrated for the generalized bivariate Bernoulli model by Islam et al. [15]. The proposed work may initiate extensive potential application of the bivariate geometric models.

2. Univariate geometric distribution

Let us consider a geometric distribution of the form (1) $p\left(y\right)=p(1-p{)}^{y-1}=p{q}^{y-1},y=1,2,\dots ,0<q<1,p=1-q,$(1) and the moment-generating function is given by $M_Y\left(t\right)=\frac{p{e}^t}{1-q{e}^t},|q{e}^t|<1.$

The expected value and variance are $E\left(Y\right)={\mu }_Y=\frac{1}{p}\mathrm{and}\mathrm{Var}\left(Y\right)=\frac{1-p}{p^2}=\frac{q}{p^2}.$

3. Bivariate geometric distribution: marginal and conditional models

The marginal or conditional distributions alone may not be sufficient to obtain the joint distribution. However, the knowledge about the marginal and conditional distributions can provide helpful strategy in determining the joint distribution. For the bivariate geometric distribution, we can employ the joint distribution based on marginal and conditional distributions introduced by Phatak and Sreehari [22] and Sreehari and Vasudeva [25] as follows: (2) $P(Y_1=y_1,Y_2=y_2)=\left(\begin{array}{c}{y}_1+y_2-2\\ y_2-1\end{array}\right)q_1^{y_1-1}{q}_2^{y_2-1}(1-q_1-q_2),$(2) (2) $y_1,y_2=1,2,\dots ,0<q_k<q_1+q_2<1,k=1,2.$(2)

Let us denote this as Model 1. The marginal and conditional distributions are (3) $\begin{array}{r}P(Y_1=y_1)=\left(\frac{1-q_1-q_2}{1-q_2}\right){\left(\frac{q_1}{1-q_2}\right)}^{y_1-1},y_1=1,2,\dots \end{array}$(3) (4) $\begin{array}{r}P(Y_2=y_2\mid Y_1=y_1)=\left(\begin{array}{c}{y}_1+y_2\\ y_2-1\end{array}\right)q_2^{y_2-1}(1-q_2{)}^{y_1},y_1,y_2=1,2,\dots \end{array}$(4) The likelihood function for the marginal model in Equation (3(2) $y_1,y_2=1,2,\dots ,0<q_k<q_1+q_2<1,k=1,2.$(2) ) is as follows: (5) $\begin{array}{rl}L& =\prod _{y_1=1}P(Y_1=y_1)=\prod _{y_1=1}^{}\left(\frac{1-q_1-q_2}{1-q_2}\right){\left(\frac{q_1}{1-q_2}\right)}^{y_1-1}\\ & =\prod _{y_1=1}\exp \left[\right(y_1-1)\ln q_1-y_1\ln (1-q_2)+\ln (1-q_1-q_2\left)\right].\end{array}$(5) The log likelihood for the marginal model is (6) $\ln L=\sum _{y_1}^{}\left[\right(y_1-1)\ln q_1-y_1\ln (1-q_2)+\ln (1-q_1-q_2\left)\right].$(6)

The exponential form for Model 1 is (7) $\begin{array}{rl}P(Y_1=y_1,Y_2=y_2)& =\exp \left[\right(y_1-1)\ln q_1+(y_2-1)\ln q_2+\ln (1-q_1-q_2)\\ & +\ln (y_1+y_2-2)!-\ln (y_1-1)!-\ln (y_2-1)!],\end{array}$(7) which can be expressed as (8) $\begin{array}{rl}P(Y_1=y_1,Y_2=y_2)& =\exp [y_1\ln q_1+y_2\ln q_2+\ln (1-q_1-q_2)-\ln q_1-\ln q_2\\ & +\ln (y_1+y_2-2)!-\ln (y_1-1)!-\ln (y_2-1)!],\end{array}$(8) where the natural parameters are $\ln q_1={\theta }_1$ and $\ln q_2={\theta }_2$. The link functions can be expressed as follows: ${\theta }_1=\ln q_1=g\left({\mu }_1\right)=X_1^{\prime }{\beta }_1,{\theta }_2=\ln q_2=g\left({\mu }_2\right)=X_2^{\prime }{\beta }_2,q_1=e^{{X^{\prime }}_1{\beta }_1},q_2=e^{{X^{\prime }}_2{\beta }_2},$ where $X_k^{\prime }=(1,x_{k1},x_{k2},\dots ,x_{kp});{\beta }_k^{\prime }=({\beta }_{k0},{\beta }_{k1},\dots ,{\beta }_{kp}),k=1,2$.

Model 1 as shown in Equation (7(7) $\begin{array}{rl}P(Y_1=y_1,Y_2=y_2)& =\exp \left[\right(y_1-1)\ln q_1+(y_2-1)\ln q_2+\ln (1-q_1-q_2)\\ & +\ln (y_1+y_2-2)!-\ln (y_1-1)!-\ln (y_2-1)!],\end{array}$(7) ) provides link functions for the outcomes but it does not contain any parameter pertaining to the association between the dependent outcomes.

The log likelihood function using Equation (8(8) $\begin{array}{rl}P(Y_1=y_1,Y_2=y_2)& =\exp [y_1\ln q_1+y_2\ln q_2+\ln (1-q_1-q_2)-\ln q_1-\ln q_2\\ & +\ln (y_1+y_2-2)!-\ln (y_1-1)!-\ln (y_2-1)!],\end{array}$(8) ) is (9) $\begin{array}{rl}\ell =\ln \left(L\right)& =\sum _{i=1}^n[y_{1i}{x}_{1i}{\beta }_1+y_{2i}{x}_{2i}{\beta }_2+\ln (1-{\mathrm{e}}^{x_{1i}{\beta }_1}-{\mathrm{e}}^{x_{2i}{\beta }_2})-x_{1i}{\beta }_1-x_{2i}{\beta }_2\\ & +\ln \left(\left(y_{1i}+y_{2i}-2\right)!\right)-\ln \left(\left(y_{1i}-1\right)!\right)-\ln \left(\left(y_{2i}-1\right)!\right)]\\ & =\sum _{i=1}^n[y_{1i}{\theta }_{1i}+y_{2i}{\theta }_{2i}+\ln (1-{\mathrm{e}}^{{\theta }_{1i}}-{\mathrm{e}}^{{\theta }_{2i}})-{\theta }_{1i}-{\theta }_{2i}\\ & +\ln \left(\left(y_{1i}+y_{2i}-2\right)!\right)-\ln \left(\left(y_{1i}-1\right)!\right)-\ln \left(\left(y_{2i}-1\right)!\right)].\end{array}$(9)

For parameter estimation, we used estimating equations and information matrix. The estimating equations are (10a) $\begin{array}{r}\frac{\mathrm{\partial }\ell }{\mathrm{\partial }{\beta }_{1j}}=\sum _{i=1}^n\left(y_{1i}-\frac{e^{{\theta }_{1i}}}{1-e^{{\theta }_{1i}}-e^{{\theta }_{2i}}}\right)X_{1i}=0,j=0,1,\dots ,p;\end{array}$(10a) (10b) $\begin{array}{r}\frac{\mathrm{\partial }\ell }{\mathrm{\partial }{\beta }_{2j}}=\sum _{i=1}^n\left(y_{2i}-\frac{e^{{\theta }_{2i}}}{1-e^{{\theta }_{1i}}-e^{{\theta }_{2i}}}\right)X_{2i}=0,j=0,1,\dots ,p.\end{array}$(10b) The second derivatives are (11a) $\begin{array}{r}\frac{{\mathrm{\partial }}^2\ell }{\mathrm{\partial }{\beta }_{1j}\mathrm{\partial }{\beta }_{1j^{\prime }}}=-\sum _{i=1}^n\left(\frac{e^{{\theta }_{1i}}}{1-e^{{\theta }_{1i}}-e^{{\theta }_{2i}}}+\frac{(e^{{\theta }_{1i}}{)}^2}{(1-e^{{\theta }_{1i}}-e^{{\theta }_{2i}}{)}^2}\right)X_{1i}^{\prime }{X}_{1i},j,j^{\prime }=0,1,\dots ,p,\end{array}$(11a) (11b) $\begin{array}{r}\frac{{\mathrm{\partial }}^2\ell }{\mathrm{\partial }{\beta }_{2j}\mathrm{\partial }{\beta }_{2j^{\prime }}}=-\sum _{i=1}^n\left(\frac{e^{{\theta }_{2i}}}{1-e^{{\theta }_{1i}}-e^{{\theta }_{2i}}}+\frac{(e^{{\theta }_{2i}}{)}^2}{(1-e^{{\theta }_{1i}}-e^{{\theta }_{2i}}{)}^2}\right)X_{2i}^{\prime }{X}_{2i},j,j^{\prime }=0,1,\dots ,p,\end{array}$(11b) (11c) $\begin{array}{r}\frac{{\mathrm{\partial }}^2\ell }{\mathrm{\partial }{\beta }_{1j}\mathrm{\partial }{\beta }_{2j}}=-\sum _{i=1}^n\left(\frac{e^{{\theta }_{1i}}{e}^{{\theta }_{2i}}}{(1-e^{{\theta }_{1i}}-e^{{\theta }_{2i}}{)}^2}\right)X_{1i}^{\prime }{X}_{2i},j=0,1,\dots ,p,\end{array}$(11c) (11d) $\begin{array}{r}\frac{{\mathrm{\partial }}^2\ell }{\mathrm{\partial }{\beta }_{2j}\mathrm{\partial }{\beta }_{1j}}=-\sum _{i=1}^n\left(\frac{e^{{\theta }_{1i}}{e}^{{\theta }_{2i}}}{(1-e^{{\theta }_{1i}}-e^{{\theta }_{2i}}{)}^2}\right)X_{2i}^{\prime }{X}_{1i},j=0,1,\dots ,p.\end{array}$(11d)

The estimates of the standard errors are obtained from the observed information based on the second derivatives as shown above. An application of the models is illustrated in Section 5.

We can use the likelihood ratio tests for $H_0:\beta =0$ against the alternative $H_1:\beta \ne 0$ for the significance of the full model where ${\beta }^{\prime }=({\beta }_{11},\dots ,{\beta }_{1p},{\beta }_{21},\dots ,{\beta }_{2p})$. Let us denote ${\beta }_f^{\prime }=({\beta }_{10},{\beta }_{11},\dots ,{\beta }_{1p},{\beta }_{20.},{\beta }_{21},\dots ,{\beta }_{2p})$ for the parameters in the full model and ${\beta }_0^{\prime }=({\beta }_{10},{\beta }_{20})$ for the parameters in the null model then the test statistic is (12) $-2[\ln (L({\hat{\beta }}_0,\beta =0)-\ln \left(L\right({\hat{\beta }}_f\left)\right]\sim {\chi }_{2p}^2.$(12)

4. Bivariate geometric distribution: joint model

In Section 3, a bivariate model is proposed on the basis of marginal and conditional distributions to specify the joint distribution. Although the proposed model in Section 3 provides the estimates of the parameters of the linear models, the association between outcome variables cannot be estimated for this model. Hence, a joint distribution for correlated geometric outcomes [21] is employed in this section which includes an association parameter. The bivariate geometric distribution of $Y_1$ and $Y_2$ can be expressed as (13) $\begin{array}{rl}P(Y_1=y_1,Y_2=y_2)& =q_1^{y_1-1}{q}_2^{y_2-1}{\theta }^{(y_1-1)(y_2-1)-1}\left[\right(1-q_1{\theta }^{y_2}\left)\right(1-q_2{\theta }^{y_1})+\theta -1\\ & =\exp [y_1\ln q_1+y_2\ln q_2+y_1{y}_2\ln \theta \\ & +\ln \left[\right\{(1-q_1{\theta }^{y_2+1})(1-q_2{\theta }^{y_1+1})+\theta -1\left\}/\theta \right]],\end{array}$(13) where $y_1=1,2,\dots ,y_2=1,2,\dots ,0<q_1,q_2<1,1+q_1{q}_2\theta \ge q_1+q_2$.

The marginal expectations for bivariate geometric distribution of $Y=\left(Y_1,Y_2\right):$ $E\left(Y_1\right)={\mu }_1=\frac{1}{1-q_1},q_1=1-\frac{1}{{\mu }_1},E\left(Y_2\right)={\mu }_2=\frac{1}{1-q_2},q_2=1-\frac{1}{{\mu }_2}.$ Let $\ln \left(q_1\right)={\theta }_1\mathrm{and}\ln \left(q_2\right)={\theta }_2,$ then $q_1=1-\frac{1}{{\mu }_1}=e^{{\theta }_1}\mathrm{and}{q}_2=1-\frac{1}{{\mu }_2}=e^{{\theta }_2}.$ From the condition $1+q_1{q}_2\theta \ge q_1+q_2,$ we have $\theta \ge \frac{q_1+q_2-1}{q_1{q}_2}=\frac{e^{{\theta }_1}+e^{{\theta }_2}-1}{e^{{\theta }_1+{\theta }_2}}.$

We can express the bivariate geometric in Equation (13(13) $\begin{array}{rl}P(Y_1=y_1,Y_2=y_2)& =q_1^{y_1-1}{q}_2^{y_2-1}{\theta }^{(y_1-1)(y_2-1)-1}\left[\right(1-q_1{\theta }^{y_2}\left)\right(1-q_2{\theta }^{y_1})+\theta -1\\ & =\exp [y_1\ln q_1+y_2\ln q_2+y_1{y}_2\ln \theta \\ & +\ln \left[\right\{(1-q_1{\theta }^{y_2+1})(1-q_2{\theta }^{y_1+1})+\theta -1\left\}/\theta \right]],\end{array}$(13) ) in terms of exponential family as (14) $\begin{array}{rl}P(Y_1=y_1,Y_2=y_2)& =\exp [y_1{\theta }_1+y_2{\theta }_2+y_1{y}_2\ln \theta -{\theta }_1-{\theta }_2-(y_1+y_2-1)\ln \theta \\ & +\ln \left[\right\{(1-e^{{\theta }_1}{\theta }^{y_2})(1-e^{{\theta }_2}{\theta }^{y_1})+\theta -1\left\}/\theta \right]],\end{array}$(14) where ${\theta }_1=\ln q_1=g\left({\mu }_1\right)=X_1^{\prime }{\beta }_1$, ${\theta }_2=\ln q_2=g\left({\mu }_2\right)=X_2^{\prime }{\beta }_2,$ $X_k^{\prime }=(1,x_{k1},x_{k2},\dots ,x_{kp})$, ${\beta }_k^{\prime }=({\beta }_{k0},{\beta }_{k1},\dots ,{\beta }_{kp}),k=1,2$, and $\theta =\frac{e^{X_1^{\prime }{\beta }_1}+e^{X_2^{\prime }{\beta }_2}-1}{e^{X_1^{\prime }{\beta }_1+X_2^{\prime }{\beta }_2}}.$

Here, ${\theta }_1$ and ${\theta }_2$ are natural link functions and θ is the association parameter of outcome variables $Y_1$ and $Y_2$ which is a function of natural link functions, ${\theta }_1$ and ${\theta }_2$. This is Model 2.

The log likelihood function using Equation (14(14) $\begin{array}{rl}P(Y_1=y_1,Y_2=y_2)& =\exp [y_1{\theta }_1+y_2{\theta }_2+y_1{y}_2\ln \theta -{\theta }_1-{\theta }_2-(y_1+y_2-1)\ln \theta \\ & +\ln \left[\right\{(1-e^{{\theta }_1}{\theta }^{y_2})(1-e^{{\theta }_2}{\theta }^{y_1})+\theta -1\left\}/\theta \right]],\end{array}$(14) ) is (15) $\begin{array}{rl}\ell & =\ln \left(L\right)=\sum _{i=1}^n{\ell }_i=\sum _{i=1}^n[y_{1i}\ln q_{1i}+y_{2i}\ln q_{2i}+y_{1i}{y}_{2i}\ln {\theta }_i\\ & +\ln \left[\right\{(1-q_{1i}{\theta }_i^{y_{2i}+1})(1-q_{2i}{\theta }_i^{y_{1i}+1})+{\theta }_i-1\left\}/{\theta }_i\right]]\\ & =\sum _{i=1}^n{y}_{1i}{\theta }_{1i}+y_{2i}{\theta }_{2i}+y_{1i}{y}_{2i}{\theta }_{3i}+\ln \left[\right\{(1-e^{{\theta }_{1i}}{\theta }_i^{y_{2i}+1})(1-e^{{\theta }_{2i}}{\theta }^{y_{1i}+1})+{\theta }_i-1\left\}/{\theta }_i\right],\end{array}$(15) where ${\theta }_{3i}=\ln {\theta }_i$.

Then the estimating equations for Model 2 are (16a) $\begin{array}{rl}\frac{\mathrm{\partial }\ell }{\mathrm{\partial }{\beta }_{1j}}& =\sum _{i=1}^n\left(y_{1i}-\frac{e^{{\theta }_{1i}}{\theta }^{y_{2i}+1}\left(1-e^{{\theta }_{2i}}{\theta }^{y_{1i}+1}\right)}{\left(1-e^{{\theta }_{1i}}{\theta }^{y_{2i}+1}\right)\left(1-e^{{\theta }_{2i}}{\theta }^{y_{1i}+1}\right)+\theta -1}\right)X_{1i}=0,j=0,1,\dots ,p;\end{array}$(16a) (16b) $\begin{array}{rl}\frac{\mathrm{\partial }\ell }{\mathrm{\partial }{\beta }_{2j}}& =\sum _{i=1}^n\left(y_{2i}-\frac{e^{{\theta }_{2i}}{\theta }^{y_{1i}+1}\left(1-e^{{\theta }_{1i}}{\theta }^{y_{2i}+1}\right)}{\left(1-e^{{\theta }_{1i}}{\theta }^{y_{2i}+1}\right)\left(1-e^{{\theta }_{2i}}{\theta }^{y_{1i}+1}\right)+\theta -1}\right)X_{2i}=0,j=0,1,\dots ,p;\end{array}$(16b) (16c) $\begin{array}{rl}\frac{\mathrm{\partial }\ell }{\mathrm{\partial }\theta }& =\sum _{i=1}^n\left[\frac{y_{1i}{y}_{2i}}{\theta }+\frac{\theta }{(1-e^{{\theta }_{1i}}{\theta }^{y_{2i}+1})(1-e^{{\theta }_{2i}}{\theta }^{y_{1i}+1})+\theta -1}\right.\end{array}$(16c) (16d) $\begin{array}{rl}& \times \left\{\frac{1}{\theta }\left(-\frac{e^{{\theta }_{1i}}{\theta }^{y_{2i}+1}(y_{2i}+1)(1-e^{{\theta }_{2i}}{\theta }^{y_{1i}+1})}{\theta }\right.\right.\\ & \left.-\frac{(1-e^{{\theta }_{1i}}{\theta }^{y_{2i}+1})e^{{\theta }_{2i}}{\theta }^{y_{1i}+1}(y_{1i}+1)}{\theta }+1\right)\\ & \left.\left.-\frac{(1-e^{{\theta }_{1i}}{\theta }^{y_{2i}+1})(1-e^{{\theta }_{2i}}{\theta }^{y_{1i}+1})+\theta -1}{{\theta }^2}\right\}\right]=0.\end{array}$(16d)

Similarly, we can obtain $\frac{{\mathrm{\partial }}^2\ell }{\mathrm{\partial }{\beta }_{1j}\mathrm{\partial }{\beta }_{1j^{\prime }}},\frac{{\mathrm{\partial }}^2\ell }{\mathrm{\partial }{\beta }_{2j}\mathrm{\partial }{\beta }_{2j^{\prime }}},\frac{{\mathrm{\partial }}^2\ell }{\mathrm{\partial }{\beta }_{1j}\mathrm{\partial }{\beta }_{2j^{\prime }}}\mathrm{and}\frac{{\mathrm{\partial }}^2\ell }{\mathrm{\partial }{\beta }_{2j}\mathrm{\partial }{\beta }_{1j^{\prime }}}.$

Similar as the likelihood ratio test as shown in Equation (12(12) $-2[\ln (L({\hat{\beta }}_0,\beta =0)-\ln \left(L\right({\hat{\beta }}_f\left)\right]\sim {\chi }_{2p}^2.$(12) ), we can test for $H_0:(\beta ,\theta )=0$ against the alternative $H_0:(\beta ,\theta )\ne 0$ for significance of the full model where $(\beta ,\theta {)}^{\prime }=({\beta }_{11},\dots ,{\beta }_{1p},{\beta }_{21},\dots ,{\beta }_{2p},\theta )$. Let us denote $({\beta }_f,\theta {)}^{\prime }=({\beta }_{10},{\beta }_{11},\dots ,{\beta }_{1p},{\beta }_{20.},{\beta }_{21},\dots ,{\beta }_{2p},\theta )$ for the parameters in the full model and ${\beta }_0^{\prime }=({\beta }_{10},{\beta }_{20})$ for the parameters in the null model then the test statistic is (17) $-2[\ln (L({\hat{\beta }}_0,\beta =0,\theta =0)-\ln \left(L\right({\hat{\beta }}_f,\theta \left)\right]\sim {\chi }_{2p+1}^2.$(17)

5. Applications

In this paper, the panel data from the Health and Retirement Study (HRS), sponsored by the National Institute of Aging (grant number NIA U01AG09740), conducted by the University of Michigan [23] is used. A total of 11 waves of the RAND version of the data from 1992 to 2012, two consecutive waves are conducted two years apart, is considered for this application. In wave one, a total of 12,762 subjects were included in the HRS cohort. Out of all these subjects, 598 cases were identified (3 cases were considered as missing) with the incidence of diabetes first followed by incidence of heart problems at a later wave. The outcome variables considered are: reported diabetes based on the diagnosis by doctor in the current or any prior wave ($Y_1$) and reported heart problem in the current or any prior wave ($Y_2$) followed by the incidence of diabetes in any wave. To make sure of the incidence cases, we have started from the respondents who did not report conditions of diabetes and heart problems in wave 1 to exclude all the prevalent cases of diabetes and heart problems in the beginning of the study. Then the cases of incidence are identified when the reporting of diabetes is found in any of the subsequent waves for the first time. Similarly, we have identified the incidence cases of heart problems when it is reported first time followed by diabetes at a later wave. It may be noted here that the incidence cases of heart problems are considered here for only those who reported health problems as diagnosed by doctors after incidence of diabetes. The heart problems include conditions such as heart attack, coronary heart disease, angina, congestive heart failure, or other heart problems. The explanatory variables considered are: age (in years), BMI (Body Mass Index), whether drink ($\mathrm{yes}=1$, $\mathrm{no}=0$) and sex ($\mathrm{male}=1$, $\mathrm{female}=0$). BMI is calculated as weight divided by the square of height where weight is converted to kilograms and height to meters. The variable drinks indicate whether the respondent drinks alcoholic beverages. After removal of cases with missing covariate values, the number of subjects is 598 who had reported diabetes in any wave followed by reported heart problems. Table 1 displays the bivariate frequency distribution of incidence of diabetes and heart problem for different waves.

Table 1. Bivariate frequency distribution of incidence of Diabetes and Heart Diseases from the HRS Data, Waves 1–10.

	Heart diseases ($Y_2$)
Diabetes ($Y_1$)
Waves	1	2	3	4	5	6	7	8	9	10	Total
1	69	42	27	36	33	35	27	26	21	10	326
2	0	2	11	6	9	10	2	5	4	3	52
3	0	0	4	4	9	12	6	4	3	5	47
4	0	0	0	4	5	6	4	5	3	5	32
5	0	0	0	0	7	3	8	6	12	4	40
6	0	0	0	0	0	6	12	6	7	1	32
7	0	0	0	0	0	0	7	9	5	7	28
8	0	0	0	0	0	0	0	10	4	8	22
9	0	0	0	0	0	0	0	0	8	8	16
10	0	0	0	0	0	0	0	0	0	3	3
Total	69	44	42	50	63	72	66	71	67	54	598

In Tables 2–4, we have shown estimates from marginal, conditional-marginal, and joint models. We have considered two models, one with the interaction term between age and sex and another without the interaction term. We have considered additional models with interaction term for age and sex to examine the role of age for females in explaining the time of incidence of diabetes and subsequently its influence in explaining the time of incidence of heart problems. We have used the time of incidence to represent the follow-up or wave number as a discrete-time variable. In the univariate geometric models, presented in Table 2, we observe that age is associated positively with both times (follow-up number represents time here) of incidence of diabetes and heart problem (p=0.000 ). The joint model based on the marginal and conditional probabilities, illustrated in Section 3, is used to estimate the parameters of the generalized linear model for the bivariate geometric model which are presented in Table 3. In this case, the marginal model is fitted for diabetes incidence and conditional model for reported heart problem for given reported diabetes. The result for the reported diabetes and heart disease appears to be similar to that of the marginal models as displayed in Table 2 for age. As we have mentioned earlier that the model shown in Section 3 does not provide any estimate for the measure of association, we can use the alternative model, shown in Section 4, for joint distribution of reported diabetes and heart problems which includes an association parameter, θ. The estimated association parameter shows a positive association between the times of incidence of diabetes and heart problems (p<0.001 ). We now observe a marked variation in the estimates for both diabetes and heart problems. The variables age, BMI and sex appear to have positive association while drinks and interaction between age and sex display statistically significant negative association for incidence of diabetes. In other words, the risk of diabetes for the males increases at higher ages as compared to that of the females in this new joint model unlike all the previous models presented in Tables 2 and 3. However, the males, compared to the females, show a significant negative association between age and time of incidence of diabetes indicating an increase in the risk for females at higher ages although still the males have a relatively higher overall risk.

Table 2. Estimates of the parameters of the univariate geometric models with or without interaction terms.

	Full model without interaction			Full model with interaction
Variables	Estimate	S.E.	p-Value	Estimate	S.E.	p-Value
Y1:Constant	−6.121	0.2401	0.000	−6.212	0.3113	0.000
Age	0.115	0.0036	0.000	0.117	0.0049	0.000
BMI	0.003	0.0032	0.317	0.003	0.0032	0.313
Drink	−0.069	0.0424	0.106	−0.069	0.0424	0.102
Sex	0.039	0.0419	0.352	0.245	0.4410	0.578
Age*Sex				−0.003	0.0073	0.639
Y2:Constant	−4.151	0.1729	0.000	−4.340	0.2230	0.000
Age	0.084	0.0023	0.000	0.087	0.0031	0.000
BMI	0.004	0.0023	0.061	0.005	0.0023	0.052
Drink	0.037	0.0313	0.240	0.037	0.0313	0.238
Sex	0.028	0.0288	0.324	0.458	0.3179	0.150
Age*Sex				−0.006	0.0047	0.175

Table 3. Estimates of the parameters of the bivariate geometric based on conditional-marginal models with or without interaction terms.

	Full model without interaction			Full model with interaction
Variables	Estimate	S.E.	p-value	Estimate	S.E.	p-value
Y1:Constant	−2.346	0.1417	0.000	−2.336	0.1874	0.000
Age	0.017	0.0020	0.000	0.017	0.0027	0.000
BMI	0.001	0.0019	0.540	0.001	0.0019	0.549
Drink	−0.029	0.0274	0.281	−0.028	0.0280	0.321
Sex	−0.007	0.0404	0.853	−0.034	0.2601	0.897
Age*Sex				0.000	0.0040	0.906
Y2:Constant	−0.782	0.1026	0.000	−0.810	0.1299	0.000
Age	0.004	0.0013	0.005	0.004	0.0017	0.019
BMI	0.001	0.0013	0.616	0.001	0.0013	0.617
Drink	0.0004	0.0153	0.978	0.001	0.0155	0.968
Sex	0.013	0.0218	0.547	0.082	0.1821	0.652
Age*Sex				−0.001	0.0026	0.704

Table 4. Estimates of the parameters of the bivariate geometric models with association parameter with or without interaction terms.

	Full model without interaction			Full model with interaction
Variables	Estimate	S.E.	p-Value	Estimate	S.E.	p-Value
Y1:Constant	−4.636	0.1782	0.000	−8.151	0.1866	0.000
Age	0.075	0.0029	0.000	0.138	0.0032	0.000
BMI	0.004	0.0007	0.000	−0.0004	0.0006	0.493
Drink	−0.082	0.0119	0.000	−0.109	0.0055	0.000
Sex	0.017	0.0081	0.033	1.373	0.1265	0.000
Age*Sex				−0.020	0.0019	0.000
Y2:Constant	−2.138	0.0814	0.000	−3.582	0.0852	0.000
Age	0.028	0.0010	0.000	0.046	0.0011	0.000
BMI	0.004	0.0004	0.000	0.013	0.0004	0.000
Drink	0.063	0.0041	0.000	0.185	0.0063	0.000
Sex	0.005	0.0041	0.257	−1.056	0.1356	0.000
Age*Sex				0.015	0.0019	0.000
$\hat{\theta }$	0.878	0.0045	0.000	0.8110	0.0041	0.000

For heart diseases, the results are similar except for the risk factor sex and the interaction between age and sex. The females appear to have higher incidence of heart diseases at a later time compared to males while there is positive association of the age and sex interaction term indicating increased incidence of heart diseases among the males compared to females as they become older. In this case, the overall risk remains higher for females, compared to males. These results reveal that the incidence of diabetes is increased at older ages but is increased for heart disease among the males compared to the females although these effects are offset to some extent due to reversal in the interaction term in both the cases.

These findings provide new insights to explain some of the ongoing debates in the relationship between incidence of diabetes and heart problem in relation to age, BMI and sex as we have indicated in the introduction. The positive association of age with the incidence of diabetes is well established and this is confirmed by the estimates obtained from the joint model as shown in Table 4. It is observed that an increase in incidence occurs with aging. Similar roles of age and BMI are also found for incidence of heart problem, although BMI does not show statistically significant relationship in the model with interaction between age and sex.

Table 5 summarizes the log likelihood and AIC for the conditional-marginal and joint models with the association parameter with or without interaction terms. In both the cases, the joint models with association parameter appear to be better in terms of goodness of fit. A comparison between the joint models with association parameter, the model with interaction term shows the best fit to the dependent outcomes for incidence of diabetes and heart diseases.

Table 5. Log Likelihood and AIC for bivariate geometric models with or without interaction terms.

	Full Model without Interaction		Full Model with Interaction
Model Type	Log Likelihood	AIC	Log Likelihood	AIC
Joint model (marginal /conditional)	−2731.54	5483.07	−2731.435	5486.87
Joint model with $\hat{\theta }$	−2576.60	5173.20	−2233.19	4490.38

6. Bootstrapping

To measure the accuracy of sample estimates, bootstrapping is used. We performed 10,000 bootstraps and computed bias, standard error, and mean-squared error for estimates. Estimates from Tables 3 and 4 are considered as population parameters while bias, standard error and mean-squared error are computed. Bias is very small generally (less than 1%) for the estimators of parameters of bivariate geometric model based on conditional marginal distributions. Standard error and mean-squared errors are also found to be very small. For the model based on the marginal/conditional approach, convergence is achieved for all 10,000 bootstrap samples. Results are presented in Table 6. Density plots for all bootstrap estimates for Model 1 are shown in Figures 1 and 2.

Figure 1. Density plot of bootstrap estimates of model without interaction for Model 1.

Figure 2. Density plot of bootstrap estimates for the model with interaction for Model 1.

Table 6. Bootstrap results from the bivariate geometric model based on conditional–marginal models with or without interaction terms.

	Full model without interaction			Full model with interaction
Variables	$\widehat{\mathrm{Bias}}$	$\widehat{\mathrm{S}.\mathrm{E}.}$	$\widehat{\mathrm{MSE}}$	$\widehat{\mathrm{Bias}}$	$\widehat{\mathrm{S}.\mathrm{E}.}$	$\widehat{\mathrm{MSE}}$
Y1:Constant	−0.02757	0.07953	0.0070848	−0.03268	0.09612	0.0103078
Age	0.00045	0.00100	0.0000012	0.00063	0.00138	0.0000023
BMI	−0.00008	0.00144	0.0000021	−0.00024	0.00161	0.0000027
Drink	0.00254	0.01746	0.0003113	0.00292	0.01782	0.0003261
Sex	0.00031	0.03548	0.0012586	−0.00028	0.13590	0.0184694
Age*Sex				−0.00003	0.00203	0.0000041
Y2:Constant	−0.00578	0.04616	0.0021638	−0.01093	0.07935	0.0064153
Age	0.00011	0.00065	0.0000004	0.00021	0.00109	0.0000012
BMI	−0.00007	0.00066	0.0000004	−0.00011	0.00080	0.0000007
Drink	−0.00022	0.00769	0.0000592	−0.00031	0.00883	0.0000780
Sex	0.00005	0.01603	0.0002570	−0.00439	0.11564	0.0133920
Age*Sex				0.00007	0.00169	0.0000028

Bootstrap results for the joint model with association parameter are presented in Table 7. For the model without interaction, 483 bootstrap samples did not converge. On the other hand, for the model with an interaction term, convergence is not achieved in case of 1839 samples. Convergence problem arises because of very few or no cell frequencies in case of bootstrap samples. Density plots for all bootstrap estimates for Model 2 are shown in Figures 3^{Figure 4}–5. Bias for the full model without interaction appears to be very small for Model 2, however, it is higher for some estimators in the full model with interaction term.

Figure 3. Density plot of bootstrap estimates for the model without interaction for Model 2.

Figure 4. Density plot of bootstrap estimates for the model with interaction for Model 2.

Figure 5. Density plot of bootstrap estimates for θ for Model 2.

Table 7. Bootstrap results from the bivariate geometric model with association parameter with or without interaction terms.

	Full model without interaction			Full model with interaction
Variables	$\widehat{\mathrm{Bias}}$	$\widehat{\mathrm{S}.\mathrm{E}.}$	$\widehat{\mathrm{MSE}}$	$\widehat{\mathrm{Bias}}$	$\widehat{\mathrm{S}.\mathrm{E}.}$	$\widehat{\mathrm{MSE}}$
Y1:Constant	−0.49503	0.74927	0.80646	−1.32840	1.71866	4.71844
Age	0.00884	00.01274	00.00024	0.02579	0.02981	0.00155
BMI	−0.00002	00.00324	00.00001	−0.00134	0.00566	0.00003
Drink	−0.00603	00.03953	00.00160	−0.03063	0.07183	0.00610
Sex	−0.01649	00.04427	00.00223	0.06377	0.93617	0.88048
Age*Sex				−0.00195	0.01474	0.00022
Y2:Constant	−0.35065	0.52177	0.39520	−1.60784	1.31771	4.32152
Age	0.00457	0.00663	0.00006	0.02079	0.01723	0.00073
BMI	0.00127	0.00345	0.00001	0.00521	0.00560	0.00006
Drink	0.02987	0.05219	0.00362	0.16579	0.12584	0.04332
Sex	0.01419	0.03391	0.00135	−0.09339	1.42762	2.04683
Age*Sex				0.00183	0.02062	0.00043
$\hat{\theta }$	−0.01140	0.01965	0.00052	−0.03403	0.03862	0.00265

7. Conclusions

In reliability or survival analysis, modeling of incidence data is a challenging issue of concern. In this paper, two bivariate models are proposed where correlated outcomes are of interest. As the geometric distribution is a natural choice for analyzing the incidence data, which has not been explored for correlated outcomes in the past, the proposed models may provide insights to the researchers and potential users with models which can be applied to the practical life situations very extensively. In addition, the estimation of parameters of the proposed models are based on generalized linear models and thus the procedure becomes easier to the potential users. Two bivariate generalized linear models are proposed in this paper for bivariate geometric distribution where two correlated geometric outcomes are considered and the corresponding link functions are shown for mean as function of covariates. An association parameter is also introduced in the bivariate generalized linear model based on a joint model which is a function of link functions for both the outcomes. An alternative model is also demonstrated using conditional and marginal distributions. However, this model does not provide any estimate of the measure of association between the outcome variables which is a major limitation, although there is no difficulty with estimation of the bivariate model for parameters for dependent outcomes. The proposed model has practical applications in survival and reliability analysis where correlated geometric outcomes are of interest.

We have applied the proposed models for analyzing the dependence in outcomes diabetes and heart disease. An effort is made to examine the effect of the explanatory variables age, BMI, drinking habit, and sex. Findings of this study provide new insights to explain some of the ongoing debates in the relationship between incidence of diabetes and heart problem in relation to age, BMI, and sex. The positive association of age with the incidence of diabetes is well established and this is confirmed by the estimates obtained from the joint model. An increased incidence of diabetes is observed with aging. Similar results are also found for age and BMI for incidence of heart problem. From our findings based on the data from only elderly population, we observe that although there is positive association between incidence of diabetes and sex, indicating higher incidence among elderly males, the age and sex interaction shows negative association. This is indicative of the fact that females start to have higher risk of diabetes compared to males with aging. It is noted by Gale and Gillespie [11] that men seem more susceptible than women to the consequences of indolence and obesity, possibly due to differences in insulin sensitivity and regional fat deposition. Among the elderly, the higher risk among males starts to offset over time for incidence of diabetes among those who end up with incidence of heart disease later. In case of heart disease for those who are already suffering from diabetes, the results are opposite. Females have higher risk of heart disease among the elderly subjects of this study which might be due to menopausal effect when the risk for heart disease increases for females. However, there is a reversal in the association between sex and time of incidence of the two diseases. From these results based on the dependence model, we can conclude that there might be delayed onset of fatal and nonfatal heart diseases in women as compared to the men [5–7] and may not necessarily be eradicated and women may still experience a higher risk of heart problem in general in the presence of diabetes mellitus. Based on the log likelihood and AIC for the conditional–marginal and joint models with the association parameter, with or without interaction, the joint models with association parameter appears to be better.

Acknowledgments

We acknowledge gratefully that the study is supported by the HEQEP subproject 3293 of the University Grants Commission of Bangladesh and the World Bank. The authors acknowledge gratefully to the HRS (Health and Retirement Study) which is sponsored by the National Institute of Aging (grant number NIA U01AG09740) and conducted by the University of Michigan for the data.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the University Grants Commission of Bangladesh and the World Bank [grant number HEQEP subproject 3293].