Bivariate lifetime geometric distribution in presence of cure fractions

In this paper, we introduce a Bayesian analysis for bivariate geometric distributions applied to lifetime data in the presence of covariates, censored data and cure fraction using Markov Chain Monte Carlo (MCMC) methods. We show that the use of a discrete bivariate geometric distribution could bring us some computational advantages when compared to standard existing bivariate exponential lifetime distributions introduced in the literature assuming continuous lifetime data as for example, the exponential Block and Basu bivariate distribution. Posterior summaries of interest are obtained using the popular OpenBUGS software. A numerical illustration is introduced considering a medical data set related to the analysis of a diabetic retinopathy data set.


Introduction
In medical, engineering or other lifetime data applications, we could have more than one lifetime associated to each unit.A special situation is the presence of two lifetimes X1 and X2 associated to each unit.This is the case, for example, considering X1 and X2 the timing of failure of paired organs like kidney, lungs, eyes, ears, dental implants among many others.In this situation, we could consider some existing bivariate lifetime distributions that has been introduced in the literature for continuous bivariate lifetime data (see for example, Arnold and Strauss, 1988;Marshall and Olkin, 1967;Sarkar, 1987;Block and Basu, 1974).
Usually these bivariate lifetime distributions generalize some popular existing lifetime distributions as exponential, Weibull, gamma or a log-normal distributions (see for example, Lawless, 1982).
Other parametrical lifetime distributions also could be generalized for the bivariate case.One of these models is given by the generalized exponential distribution (see for example, Gupta andKundu, 2007, Mohsin et al., 2013).A bivariate generalized exponential distribution also was introduced by Kundu and Gupta (2008).Among all these continuous bivariate lifetime distributions, one model has been very well explored in the literature: the Block and Basu (1974) bivariate exponential distribution, whose marginals are mixtures of exponentials and having an absolutely continuous joint distribution with positive correlation coefficient.This distribution was derived originally by omitting the singular part of the Marshall and Olkin (1967) distribution.The Marshall and Olkin distribution is not absolutely continuous but satisfies the lack of memory property.
It is important to point out that censoring is a key aspect of survival analysis.In this way, many recent papers related to the analysis of bivariate lifetime data are related to the presence of censored data (see for example As an alternative for the use of bivariate parametrical distributions for bivariate continuous lifetime data, we could assume that the lifetimes X1 and X2 are discrete random variables assuming any positive integer that is, we round the lifetimes given with decimal part to an integer.It is important to point out, that despite discrete measuring of medical or engineering lifetime data is very common in applications, few papers related to discrete lifetime data are introduced in the literature, (see for example, Grimshaw and et al., 2005;Davarzani and Parsian, 2013).
In this way, Arnold (1975) introduced a general multivariate geometric distribution which showed that it leads in a natural way to the Marshall-Olkin multivariate exponential distribution.Nair and Nair (1988) studied the characterizations of bivariate exponential distributions and geometric distributions.
In this paper, we explore the use of a bivariate geometric distribution to analyze bivariate lifetime data.The bivariate geometric distribution proposed by Arnold (1975) has joint mass probability given by, (1) where the marginal probability mass functions are respectively, given by standard geometric distributions starting at one, that is, and which means, variances, covariance and correlation given respectively, by where r = 1 -θ1 -θ2 ; 0 < θ1 < 1 and 0 < θ2 < 1.
From (1), we observe that this bivariate discrete model assumes that P (X1 = X2) = 0, but in applications we could have this probability very small, but different of zero, a restrictive property of the model.In some cases, however, like studying the time of blindness of eyes or infection of kidneys, this probability practically is zero.In such cases the proposed model will be more fit to the data rather that the models with positive probability of X1 = X2.
In the presence of right censored data and covariates, a common situation in applications, we could have some computational difficulties to get standard classical inferences for the parameters of bivariate lifetime distributions, as for example, using maximum likelihood approach.As an alternative for the use of standard classical inference methods, the use of a Bayesian approach is becoming very popular in the analysis of bivariate lifetime data (see for example, Achcar and Leandro, 1998; or Santos and Achcar, 2011).
This paper is organized as follows: in section 2, we introduce the presence of censored data; in section 3, we introduce the presence of cure fraction; in section 4, we introduce an analysis of a diabetic retinopathy data set; finally, in section 5, we present some concluding remarks.

Presence of cure fraction
In lifetime data analysis the presence of censored observations is very common in applications.The censored data could be related to individuals lost to follow-up or that will never experience the event of interest.This situation could occur in different areas, such as in cancer studies where the researchers are interested in the proportion of patients cured and where many individuals will die of other causes.In practical work we should carefully interpret the results using a cure fraction model when this situation is not true.As an example of this case, disease recurrence in primary breast cancer; even after many years, the condition can recur.A detailed discussion on the presence of cure fraction is introduced by Lambert (2007).As pointed out by this author, information of cure at the individual level will rarely be available, and so in these models we are concerned with population (or statistical) cure.If reliable information on cause of death is available, then one can perform a cause-specific analysis where deaths not due to the disease of interest can be treated as censored observations.In many situations, the cause of death may either not be recorded or obtained from death certificates, which are often inaccurately recorded (Begg and Schrag, 2002).A possibility in this case: to obtain the expected survival and/or the expected mortality rate from national mortality statistics, and such is usually calculated after matching for age, sex, year of diagnosis, and possibly other covariates (Coleman et al,1999).
In this situation, after a careful preliminary data analysis, we could have a fraction of individuals not expecting to experience the event of interest, that is, these individuals are not at risk ("long term survivors" or "cured individuals").Different approaches have been considered to model cure fraction, especially for univariate lifetime data (see for example

Analysis of a diabetic retinopathy data set
In this application, we consider a data set introduced by Huster et al (1989).In this study, we have 197 patients with "high-risk" diabetic retinopathy as defined by the Diabetic Retinopathy Study (DRS).Each patient had one eye randomized to laser treatment and the other eye received no treatment.For each eye, the event of interest was the time from initiation of treatment to the time when visual acuity dropped below 5/200 two visits in a row (call it "blindness").Thus there is a builtin lag time of approximately 6 months (visits were every 3 months).Survival times in this dataset are therefore the actual time to blindness in months, minus the minimum possible time to event (6.5 months).Censoring was caused by death, dropout, or end of the study (Huster et al, 1989).Associated to each individual we have two covariates: Z1 denoting the age at diagnosis of diabetes and Z2 denoting the type of diabetes where 1= juvenile (age at diagnosis < 20) and 0 = for adults.

A preliminary data analysis
In Figure 1  From the plots of Figure 1, we observe some indication of cure fraction, or an indication that a great number of individuals that will never become blind.

A Bayesian analysis
To analyse the diabetic retinopathy data set under a Bayesian approach, let us first assume the bivariate geometric distribution with density (1) since both responses (X1 = time to blindness for the treated eye and X2 = time to blindness for the untreated eye) are associated to the same individual.Observe that in this case, we are transforming both continuous lifetimes to discrete lifetimes to assume a bivariate geometric distribution.As a first analysis, we do not consider the presence of the covariates Z1 denoting the age at diagnosis of diabetes and Z2 denoting the type of diabetes where 1= juvenile (age at diagnosis < 20) and 0 = for adults and not considering the presence of cure fraction.Let us denote this model as "model 1".
For a Bayesian analysis of "model 1", let us assume a Dirichlet prior distribution for θ1, θ2 and r with hyper parameter values α1 = α2 = α3 = 1, given by; (6) Using the OpenBUGS software (Spiegelhalter et al, 2003) in the simulation of samples of the joint posterior distribution for θ1, θ2 and r we discarded the first 10,000 simulated Gibbs samples (burn-in-samples) to eliminate the effect of the initial values for the parameters of the model.Choosing every 50th simulated Gibbs sample, we obtained a final sample of size 1,000 to get the posterior summaries of interest.Convergence of the Gibbs sampling algorithm was monitored using existing methods as time series plots for the simulated samples.In Table 1, we have the posterior means, the posterior standard-deviations and 95% credible intervals for the parameters for θ1, θ2 and r.We also have in Table 1, the posterior summaries for the means, standard-deviations and for the correlation coefficient (see ( 1)) for the lifetimes X1 and X2.
In the presence of the covariates Z1 denoting the age at diagnosis of diabetes and Z2 denoting the type of diabetes where 1=juvenile (age at diagnosis < 20) and 0 = for adults, we now assume the bivariate geometric distribution with density (1) and the following logistic regression model for the parameters θ1 and θ2, given by, (7) where i = 1, 2, ..., 197; logit(θ) = log[θ/(1 − θ)], Z ̅ 1 is the average of the covariate age.Let us denote this model as"model 2".Observe that ri = 1 − θ1i -θ2i.Assuming normal prior distributions N (0, 10) for the regression parameters βj0, j = 1, 2, and normal prior distributions N (0,1) for the regression parameters βjl , j = 1, 2; l = 1, 2, and using the OpenBUGS software also considering a "burn-in-sample" of size 10,000 and final sample of size 1,000 taking every 50th sample, we have in Table 2, the posterior summaries of interest.
It is important to point out that using the approximately normal non-informative N(0,10) for the regression parameters βj0, j = 1, 2, and N (0, 1) for the regression parameters βjl , j = 1, 2; l = 1, 2, using the OpenBugs software, we obtained very good convergence of the MCMC simulation algorithm as observed in trace plots of the simulated Gibbs samples.Other very non-informative priors also have been considered as normal N (0, 10) priors for all regression parameters but the obtained posterior summaries were very close to the obtained results given in Table 2, that is, the obtained posterior summaries are not sensible to the choice of other priors, leading to the same inference results.From the results of Table 2, we observe that zero is included in the 95% credible interval for β11, β12, β21, and β22, that is, the lifetime X1 and X2 are not affected by the covariates Z1 denoting the age at diagnosis of diabetes and Z2 denoting the type of diabetes where 1=juvenile (age at diagnosis < 20) and 0 = for adults.
As a third modeling approach, not considering the presence of the covariates Z1 denoting the age at diagnosis of diabetes and Z2 denoting the type of diabetes where 1= juvenile (age at diagnosis < 20) and 0 = for adults, we assume the bivariate geometric distribution with density (1) considering the presence of cure fraction (see, section 3).Let us denote this model as"model 3".Table 3: Poterior summaries "model 3".
For a Bayesian analysis of "model 3", let us assume a Dirichlet prior distribution (6) for θ1 , θ2 and r with hiperparameter values α1 = α2 = α3 = 1 and another Dirichlet prior distribution for ϕ11, ϕ10 , ϕ01 and ϕ00 with hyper-parameter values α1 = α2 = α3 = α4 = 1.Note that with these values of hyper-parameters, we are assuming approximately non-informative priors for θ1, θ2 and r and for ϕ11, ϕ10 , ϕ01 and ϕ00 .Using the OpenBUGS software also considering a "burn-in-sample" of size 10,000 and final Gibbs sample of size 1,000 taking every 50th sample, we have in Table 3, the posterior summaries of interest.
From the results of Table 3, we observe that the Monte Carlo estimate for the posterior mean for ϕ10 (probability of observations pairs with V1 = 1 for susceptible and V2 = 0 for immune) based on the 1,000 simulated Gibbs samples is approximately equal to 0.67 and for ϕ11 (probability of V1 = 1 and V2 = 1) is approximately equal to 0.28 with large 95% credible intervals.It is important to point out that more accurate Bayesian inferences could be obtained using more informative priors for the parameters ϕ0 , ϕ01, ϕ10 and ϕ11 based on expert opinions or using empirical Bayesian methods.
Finally, considering the presence of the covariates Z1 denoting the age at diagnosis of diabetes and Z2 denoting the type of diabetes where 1 = juvenile (age at diagnosis < 20) and 0 = for adults, we assume the bivariate geometric distribution with density (1) and regression model (3) considering the presence of cure fraction (see, section 3).Let us denote this model as "model 4".
For a Bayesian analysis of "model 4", let us assume a Dirichlet prior distribution for ϕ11, ϕ10, ϕ01 and ϕ00 with hyper parameter values α1 = α2 = α3 = α4 = 1 and the regression model (3) with the following priors for the regression parameters: β10 ∼ N (−5, 1), β11 ∼ N (0.005, 1), β12∼ N (0.5, 1), β20 ∼ N (−4, 1), β21 ∼ N (−0.02, 1) and β22 ∼ N (−0.4,1).Observe that we are using informative priors for the regression parameters based on the information obtained from "model 2" (Use of empirical Bayesian methods, see for example, Carlin and Louis, 2000.)Using the OpenBUGS software also considering a "burn-in-sample" of size 10,000 and a final Gibbs sample of size 1,000 taking every 50th sample, we haveing Table 4, the posterior summaries of interest.From the results of Table 4, we also observe that zero is included in the 95% credible interval for β11 , β12 , β21, and β22 , that is, the lifetime X1 and X2 are not affected by the covariates Z1 denoting the age at diagnosis of diabetes and Z2 denoting the type of diabetes where 1= juvenile (age at diagnosis < 20) and 0 = for adults.Similar Bayesian inferences were obtained for the parameters ϕ00, ϕ01, ϕ10 and ϕ11 (see Table 3).
As a comparison for the four proposed models, we could use a Bayesian discrimination criterion, as for example, the DIC (Deviance Information Criterion) introduced by Spiegelhalter et al. ( 2002), (see the appendix) and given automatically by the OpenBUGS software.Assuming "model 1" (bivariate geometric distribution without the presence of covariates and cure fraction), we have a Monte Carlo estimate for DIC given by the value 1684.0;assuming "model 2" (bivariate geometric distribution in presence of covariates but not considering presence of cure fraction) the DIC value is given by 1679.0;assuming "model 3" (bivariate geometric distribution not considering the presence of covariates but considering the presence of cure fraction) the DIC value is given by 1553.0 and assuming "model 4" (bivariate geometric distribution considering the presence of covariates and the presence of cure fraction) the DIC value is given by 1556.0.That is, "model 3" and "model 4" (presence of cure fraction) give better fit for the diabetic retinopathy data set (smaller value for DIC).Since both covariates do not show significant effects on the lifetimes X1 and X2, and the DIC values for "model 3" and "model 4" are very similar, we could conclude that "model 3" is the best model to be fitted by the data set.

Concluding remarks
The search for new lifetime models is of great interest for researchers in statistics and in different areas of applications as medical or engineering studies.In some situations, we have bivariate lifetime data, that is, two lifetimes are measured for each unit where it is important, to have a model that captures the possible dependence between both responses.Assuming continuous lifetimes, the literature presents many bivariate parametrical distributions, as for example, the popular Block-Basu (1974) and the Marshall Olkin (1967) distributions among many others.A Bayesian analysis of the Block and Basu distribution in presence of cure fraction is introduced by Achcar et al., 2013.As an alternative, we propose the use of a bivariate geometric distribution which could be a new and suitable alternative for the analysis of bivariate lifetime data, especially under a Bayesian approach and using MCMC methods.The use of discrete probability distributions still is not well explored in the literature for the analysis of lifetime data.This family of distributions could bring us some improvements in the computational work to obtain the inferences of interest.In many applications, especially in medical studies, we could have the presence of censored data, the presence of one or more covariates and the presence of cure fraction, a situation that could bring us great computational difficulties to get the usual classical inferences for the parameters of the proposed model like the convergence of the iterative numerical algorithm used to get the maximum likelihood estimators for the parameters of the model.It is also important to point out that these classical inference procedures are based on the asymptotical normality of maximum likelihood estimators (see for example, Lawless,1982).This could be a problem when we have small data sizes, a situation very common in medical studies.In this way, we use Bayesian methods to analyze bivariate lifetime data assuming a bivariate geometric distribution.Also it is important to point out that the use of available Bayesian software like the OpenBUGS software, only requires the specification of the distribution for the data and prior distributions for the parameters of the model, giving us a great simplification to obtain posterior summaries of interest, as it was observed in the application considering a medical data set introduced in section 4.These results could be of great interest for applications in medical and engineering studies.
, Yu et al, 2004; Gamel et al, 1999; Yamaguchi, 1992; Cancho and Bolfarine, 2001; Taylor, 1995; Kannan et al, 1999).Wienke et al (2006) introduced a model for a cure fraction in bivariate time-to-event data.Let us assume that the population is divided in two groups of individuals: a group of cured individuals with probability 1 − ϕ and a group of susceptible individuals with a proper survival function S(x) = P (X > x), where X denotes the discrete lifetime of the individual with probability ϕ.In this way, we have a survival function equals to one for all X considering the cured individuals.Considering univariate lifetimes, a model that incorporates a cure fraction (see for example, Wienke et al., 2006) is given by, (4) where p ∈ (0, 1) is the mixing parameter and S0(x) denotes a proper survival function for the non-cured group in the population.A generalization of (3) considering bivariate lifetimes X1 and X2 is (see Wienke et al., 2006) given by, (5) where S(x1, x2) = P (X1 > x1, X2 > x2 ) is the joint survival function; S(x1) = P (X1 > x1) is the marginal survival function for X1; S(x2) = P (X2 > x2) is the marginal survival function for X2; ϕ11 = P (V1 = 1, V2 = 1); ϕ10 = P (V1 = 1, V2 = 0); ϕ01 = P (V1 = 0, V2 = 1); ϕ00 = P (V1 = 0, V2 = 0); ϕ11 + ϕ10 + ϕ01 + ϕ00 = 1; V1 and V2 are binary variables such that V1 = 1 if the individual is susceptible for lifetime X1 and that V1 = 0 if the individual is immune; in the same way, V2 = 1 if the individual is susceptible for lifetime X2 and that V2 = 0 if the individual is immune.From the above definitions, the likelihood function (2) in presence of cure function and right censoring is given by, analysis, we assume a Dirichlet prior for the parameter ϕ11, ϕ10 , ϕ01 and ϕ00.The Dirichlet prior is used since ϕ11 + ϕ10 + ϕ01 + ϕ00 = 1.
, we have the plots of the non-parametrical Kaplan and Meier (1958) estimators of the survival functions for both times (X1 = time to blindness for the treated eye and X2 = time to blindness for the untreated eye).The non-parametrical estimates for the means (MTTF) of both survival functions obtained from the Kaplan-Meier estimates are given, respectively, by 53.7284 (for X1, with 54 uncensored observations and 143 right censored observations) and 43.5258 (for X2 , with 101 uncensored observations and 96 right censored observations).