Bivariate Weibull Distributions Derived From Copula Functions In The Presence Of Cure Fraction And Censored Data

In this paper we introduce bivariate Weibull distributions derived from copula functions in presence of cure fraction, censored data and covariates. Two copula functions are explored: the FGM (Farlie Gumbel Morgenstern) copula and the Gumbel copula. Inferences for the proposed models are obtained under the Bayesian approach, using standard MCMC (Markov Chain Monte Carlo) methods. An illustration of the proposed methodology is given considering a medical data set.


Introduction
A long-term survivor mixture model, also known as standard cure rate model, assumes that the studied population is a mixture of susceptible individuals, who experience the event of interest and non susceptible individuals that will never experience it.These individuals are not at risk with respect to the event of interest and are considered immune, non susceptible or cured Maller and Zhou (1996).Different approaches, parametric and non-parametric, have been considered to model the proportion of immunes and interested readers can refer, for example, to Boag (1949), Berkson (1952), Haybittle (1965), Farewell (1982Farewell ( , 1986)), Meeker (1987), Dunsmuir et al., (1989), Taylor, (1995), Gamel et al. (1990Gamel et al. ( , 1999)), Ghitany and Maller (1992), Yamaguchi, (1992), Copas and Heydary (1997), Ng and McLachlan (1998), De Angelis et al. (1999), Peng and Dear (2000), Sy and Taylor (2000), Cancho and Bolfarine, (2001), Yu et al. (2004), Kannan et al., (2010).Cure models for paired and clustered survival data have been considered by a number of authors, Wienke et al. (2003Wienke et al. ( ,2006) ) introduced a model for a cure fraction in bivariate time-to-event data motivated by the article of Chatterjee and Shih (2001).Chatterjee and Shih (2003) used a correlated gamma frailty model by a copula function to analyse bivariate survival data.
Considering univariate lifetimes and following Maller and Zhou (1996), the standard cure rate model assumes that a certain fraction p in the population is cured or never fail with respect to the specific cause of death or failure, while the remaining (1 − p) fraction of the individuals is not cured, leading to the survival function for the entire population written as, S(t) = p + (1 − p)S0(t), (1) where p ∈ (0,1) is the mixing parameter and S0(t) denotes a proper survival function for the non cured group in the population.
From the mixture survival function, (1), the probability density function is obtained from f (t) = −dt d S(t) and given by, where f0 (t) is the probability density function for a susceptible individual.From , the hazard function is given by, .
(3) An alternative to a standard long-term survivor mixture model is the longterm survivor non-mixture model suggested by Yakovlev et al. (1996); Tsodikov et al. (2003); Lambert et al. (2010) which defines an asymptote for the cumulative hazard and hence for the cure fraction.The non-mixture model or the promotion time cure fraction has been used by Lambert et al. (2007Lambert et al. ( , 2010) ) to estimate the probability of cure fraction in cancer lifetime data.
In this paper we consider a generalization of the mixture cure fraction model (1) for the bivariate case considering different copula functions to capture the possible existing dependence between two lifetimes T1 and T2 for the susceptible individuals, censored data and the presence of covariates.Inferences for the proposed models are obtained using Bayesian methods where samples of the joint posterior distribution of interest is simulated using MCMC (Markov Chain Monte Carlo) methods, as the popular Gibbs sampling algorithm Gelfand and Smith (1990); Casella and George (1992) and the Metropolis-Hastings algorithm Chib and Greenberg (1995).
The paper is organized as follows: in Section 2 we introduce bivariate lifetime data assuming mixture cure fraction models; in Section 3 we present the likelihood function considering bivariate lifetime data in presence of cure fraction and censored data considering two special copula functions: Farlie-GumbelMorgenstern (FGM) copula and Gumbel copula; in Section 4 we introduce the Weibull marginal distributions for the lifetimes T1 and T2; in Section 5 we introduce an illustrative example.Finally, in Section 6 we present some concluding remarks.

Bivariate Lifetime Data Assuming Mixture Cure Fraction Models
In some areas of application, especially in medical and engineering studies, we could have two lifetimes T1 and T2 associated with each individual or unit.Usually, these data are assumed to be independent, but in many cases, the lifetime of one component could affect the lifetime of the other component.This is the case as an example in the medical area of paired organs like kidneys, lungs, eyes, ears, dental implants among many others.In the literature we observe many papers related to bivariate lifetime parametric models Gumbel (1960); Freund (1961); Marshall and Olkin (1967); Downton (1970); Hawkes (1972); Block and Basu (1974); Hougaard (1986); Sarkar (1987); Arnold and Strauss (1988).An alternative is the use of copula functions Nelsen (2006); Trivedi andZimmer (2005, 2006) to model bivariate lifetime data assuming different marginal lifetime distributions.
We also could have the presence of cure fraction or immunes for both lifetimes T1 and T2.In this way, let us assume mixture models for T1 and T2, given, respectively, by, where S10 (t1) and S01 (t2) are the survival functions for the susceptible individuals in the lifetimes T1 and T2, respectively; p1 and p2 are, respectively, the cure probabilities in T1 and T2; (1 − p1) and (1 − p2) are, respectively, the probabilities of susceptible in lifetimes T1 and T2.
It is important to point out that any multivariate distribution function F can be written in the form of a copula function Sklar (2010), that is, if F (t1,t2,...,tm) is a joint multivariate distribution function with univariate marginal distribution functions F1 (t1),F2 (t2),...,Fm (tm), thus there exist a copula function C (u1,u2,...,um), such that, If every Fl is continuous, then C is unique.For the special case of bivariate distributions, we have m = 2.
The approach to formulate a multivariate distribution using a copula is based on the idea that a simple transformation can be made of each marginal variable in such a way that each transformed marginal variable has a uniform distribution.Once this is done, the dependence structure can be expressed as a multivariate distribution on the obtained uniforms, and a copula is precisely a multivariate distribution on marginally uniform random variables.In this way, there are many families of copulas which differ in the detail of the dependence they represent.In the bivariate case, let T1 and T2 be two random variables with continuous distribution functions F1 and F2.
The probability integral transform can be applied separately to the two random variables to define U = F1 (t1) and V = F2 (t2), where U and V have uniform (0,1) distributions, but are usually dependent if T1 and T2 are dependent (T1 and T2 independent, implies that U and V are independent).Specifying dependence between T1 and T2 is the same as specifying dependence between U and V .With U and V uniform random variables, the problem reduces to specifying a bivariate distribution between two uniforms, that is, a copula.

Bivariate Lifetime Data In Presence Of Cure Fraction And Censored Data
Let us assume n pairs of a random sample of lifetimes t1i and t2i, i = 1,...,n, that can be right censored where censoring is independent of the lifetimes divided into four classes: C1 : both lifetimes t1i and t2i are observed; C2 : t1i is a lifetime and t2i is a censored time (that is, we only know that T2i ≥ t2i); C3 : t1i is a censoring time and t2i is a lifetime (that is, we only know that T1i ≥ t1i); C4 : both t1i and t2i are censoring times.
It is important to note that we could have other censoring schemes but this will be not considered in this paper.
The contribution of the i th individual for the likelihood function Lawless (1982) is given by, , (21) where f (t1i,t2i) is the joint probability density function for T1 and T2; S (t1i,t2i) is the joint survival function; and are the partial derivatives of S (t1i,t2i), with respect to t1i and t2i, respectively.
Let us define the indicator variables δ1i and δ2i, by, In this way, we rewrite the likelihood function (21) as, (23)
Different marginal parametric lifetime distributions could be used for the lifetimes T1 and T2, as exponential, Weibull, log-normal, gamma or generalized gamma distributions.

Weibull Marginal Distributions For The Lifetimes T1 And T2
Considering the bivariate lifetime distributions in presence of cure fractions derived from mixture models with a dependence structure given by FGM or Gumbel copula functions, we assume as a special model, a popular lifetime model extensively used in medical or engineering applications: the Weibull distribution.In this way, we consider marginal Weibull distributions for the lifetimes T1 and T2 of susceptible individuals with densities and survival functions given, respectively, by, ; ; ; .
The Weibull distribution is very popular to analyse lifetime data, since its hazard function h(t) = f (t)/S (t) is increasing if the shape parameter λ > 1; is decreasing if λ < 1 or is constant (an exponential distribution) if λ = 1 Lawless (1982).

An Example
In this application, we consider a data set Huster et al. (1989) with 197 patients were a 50% random sample of the 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 built-in lag time of approximately 6 months (visits were every 3 months).Survival times in this data set 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.This data set can be obtained in http://www.mayo.edu/research/documents/diabeteshtml/DOC-10027460 In the analysis considered here, we assume as lifetimes, the times (in months) to blindness for the eye randomized to laser treatment (T1), with 143 censored observations and 54 notcensored observations, and the time (in months) to blindness for the eye randomized that not received the treatment (T2), with 96 censored observations and 101 not-censored observations.Two covariates, the age at diagnosis of diabetes (in years) and the type of diabetes (juvenile, adult) were considered.Observe that among the censored observations, some patients will never have the occurrence of the event of interest, that is, they are immunes.
As a preliminary analysis of the data, we present in Figure 1, the KaplanMeier nonparametric estimates for the survival functions considering the two lifetimes associated to each patient (treated and untreated eyes).From these plots, we have some indication of cure fraction (a large number of censoring times at the end of study).From the Kaplan-Meier plots of Figure 1, we also observe better results for the treated eyes, that is, larger times to blindness in comparison with the untreated eyes.From the Kaplan-Meier estimates we have estimates for the mean time to blindness: 57.83 months for the treated eyes and 43.53 months for the untreated eyes.
In medical research, it is common to assume independence between T1 and T2, which could be not appropriated.As a first model, we assume a bivariate Weibull lifetime distribution derived from the FGM copula in presence of cure fraction presented in Section 3.1 with independent priors λl ∼ Gamma(0.001,0.001),µl ∼ Gamma(0.001,0.001),l = 1,2, and θ1 ∼ U (−1,1), were Gamma(a,b) denotes a gamma distribution with mean a/b and variance a/b 2 and U (c,d) denotes an uniform distribution defined in the interval (c,d).We also assume a Dirichlet prior given by ( 38) for φ11, φ10, φ01 and φ00 with hyperparameter values w1 = w2 = w3 = w4 = 1 (noninformative priors for the incidence parameters).
Considering the rjags library Plummer (2011) in R software R Development Core Team (2011), a single chain has been used considering the simulation for each parameter for 200,000 times with a burn sample size of 50,000 to eliminate the possible effect of the initial values.We choose the simulated values from 100 to 100, to get approximated uncorrelated values which result in a final chain of size 2,000.We automatically generate initial values by sampling from the prior for the parameters, since the priors are not totally diffuse (vague).Convergence of the Gibbs Sampling algorithm was monitored using standard graphical methods, as the trace plots and the estimated posterior densities of the simulated samples Brooks and Gelman (1998).In Figure 2 we have an example considering the parameter λ1 from the FGM copula (first fitted model), all the other parameters have similar behavior.Note we did not show the full conditional distributions since we consider the rjags library, making it unnecessary to use, however in appendix A we show the required codes to obtain the Bayesian estimates.In Table 1, we have the posterior summaries of interest (Monte Carlo estimates of the posterior means for each parameter).We also have in Table 1,   In Table 2, we have the posterior summaries of interest assuming independent marginal Weibull distributions for the lifetimes T1 and T2 in presence of cure fraction, that is, with θ1 = 0, assuming the same priors for the parameters of interest and the same Gibbs sampling steps using the rjags library in R software used for the bivariate Weibull distribution derived from FGM copula in presence of cure fraction (results in Table 1).
From the obtained Monte Carlo estimates for DIC assuming the two models given in Tables 1 and 2 (presence or not of dependence between T1 and T2) we observe that the bivariate Weibull distribution derived from a FGM copula function gives better fit for the data, since the obtained DIC value is smaller for the dependent model than for the independent model.Considering "FGM copula model" in the presence of cure fraction and the covariates age and type of diabetes (see (40)), that is, where i = 1,...,n (n = 197) x1i denotes the age at diagnosis of diabetes (in years) and x2i is a binary variable denoting the type of diabetes (1 = juvenile; 0 = adult).We consider the same priors for the parameters, λ1, λ2, θ1, φ11, φ10, φ01 and φ00, and the same Gibbs sampling steps using the rjags library in R software used for the bivariate Weibull distribution derived from FGM copula in presence of cure fraction.For the others parameters we consider independent priors αl ∼ Gamma(0.001,0.001),βl1 ∼ N (0,1000) and βl2 ∼ N (0,1000), l = 1,2, were N (e,f) denotes a normal distribution with mean e and variance f.In Table 3, we have the posterior summaries of interest assuming "FGM copula model" in the presence of cure fraction and the covariates age and type of diabetes.
From the results of Table 3, we observe that both covariates (age at diagnosis and type of diabetes) do not present significant effects on the response since the 95% credible intervals contain the zero value.In fact the Monte Carlo estimates of DIC are very similar considering the dependent model in presence or not of covariates.Considering the independent model, we observe a DIC value a little bit larger, an indication that the dependence model is better fitted by the data.In this case, we also could to model the data considering other existing copula functions since the FGM copula is appropriated only to model weak dependence Nelsen (2006).Observe that the 95% credible interval of the dependence parameter θ1 includes the zero value, implying that θ1 is insignificant and it is enough to justify that T1 and T2 are not dependent at least through FGM copula function.Now, let us assume a bivariate Weibull lifetime distribution derived from the Gumbel copula in presence of cure fraction presented in Section 3.2 with independent priors λl ∼ Gamma(0.001,0.001),µl ∼ Gamma(0.001,0.001),l = 1,2, and θ2 ∼ U (0,1).We also assume a Dirichlet prior given by (38) for φ11, φ10, φ01 and φ00 with hyperparameter values w1 = w2 = w3 = w4 = 1 (noninformative priors for the incidence parameters).
Considering the rjags library Plummer (2011) in R software R Development Core Team (2011), a single chain has been used considering the simulation for each parameter for 210,000 times with a burn sample size of 50,000 to eliminate the possible effect of the initial values.We choose the simulated values from 100 to 100, to get approximated uncorrelated values which result in a final chain of size 2,000.Convergence of the Gibbs Sampling algorithm was monitored using standard graphical methods, as the trace plots of the simulated samples Brooks and Gelman (1998).
In Table 4, we have the posterior summaries of interest (Monte Carlo estimates of the posterior means for each parameter).We also have in  Finally, assuming "Gumbel copula model" in presence of cure fraction and covariates, we assumed the same regression model (43) and consider the same priors for the parameters, λ1, λ2, θ1, φ11, φ10, φ01 and φ00, and the same Gibbs sampling steps using the rjags library in R software used for the bivariate Weibull distribution derived from Gumbel copula in presence of cure fraction.For the others parameters we consider independent priors αl ∼ Gamma(0.001,0.001),βl1 ∼ N (0,1000) and βl2 ∼ N (0,1000), l = 1,2.In Table 5, we have the posterior summaries of interest assuming "Gumbel copula model" in the presence of cure fraction and the covariates age and type of diabetes.From the results of Table 5, we also conclude that both the covariates do not have significant effects on the response, but the dependent model derived from Gumbel copula gives smaller values of DIC, possibly an indication of better fit of the model for the data.Since, in this case, it is not possible to compare the Kaplan-Meier curves with the estimated survival functions, we could use other discrimination methods to choose the best model for the data set.A possibility to discriminate the proposed models in terms of better fit for the data not considering the presence of covariates is to compare the Monte Carlo estimates of the means for the lifetimes T1 and T2 using the bivariate Weibull distributions derived from copula functions.Standard non-parametric estimates for the means based on Kaplan-Meier estimates are given, respectively, by E (T1) = 57.83and E (T2) = 43.53(use of the R software).Using the bivariate Weibull distributions derived from copula functions in presence of cure fraction, we get estimates for the means of T1 and T2 from (37), given by:  FGM copula in presence of cure fraction not considering the presence of covariates: E (T1) = 36.08and E (T2) = 25.62. Independent lifetimes with Weibull distributions in presence of cure fraction not considering the presence of covariates: E (T1) = 22.42 and E (T2) =23.87. Gumbel copula in presence of cure fraction not considering the presence of covariates: E (T1) = 36.87and E (T2) = 19.28.From these results, we also observe that the dependent models based on the copula functions are better fitted by the data, as compared to the independent Weibull model, since the obtained Bayesian estimates using copulas for the means of the marginal lifetimes are more close to the Kaplan-Meier estimates.

Concluding Remarks
The use of copula functions could be a good alternative to analyse bivariate lifetime data in presence of censored data, cure fraction and covariates.Observe that in many applications of lifetime modeling we could have the presence of a cure fraction for individuals that are "long term survivors" or "cured individuals".
In this paper, we have used two copula functions, but we also could consider any other existing copula function to build new bivariate lifetime models.We also have used marginal Weibull distributions, since this distribution has a good flexibility of fit for lifetime data, as constant, increasing or decreasing hazard functions.Other usual existing univariate lifetime distributions could be used (Gamma, log-normal, exponential, exponentiated Weibull distributions, among many others) to model the marginal distributions.This approach could be of great interest in many areas of interest, especially in medical or engineering studies.
The use of standard existing MCMC methods, especially with the rjags library in R software, could give a great simplification to obtain inferences for the proposed models.Using a real data set example, we observed that using both proposed models (FGM or Gumbel copulas) give similar results when observing the Monte Carlo estimates for the posterior means and credible intervals for all parameters, except for the association parameters θ1 and θ2, which have different interpretations Nelsen (2006).A small improvement for the fit of the models for the data is observed considering "Gumbel copula model" where the DIC values are smaller than considering "FGM copula model" (see Tables 1-5).Observe that prior information on the association parameter of copula functions usually it is not available from preliminary studies.
It is interesting to observe that we have great difficulties to get standard classical inferences for the parameters of the proposed models.In our case, we tried to find maximum likelihood estimators for the parameters of the models in our application using SAS software, but we did not get good results.Other problem with the use of classical approach: confidence intervals and hypothesis tests are based on asymptotic results.In our case, these results could be not accurate even using large sample sizes.
It is important to point out that any existing parametric lifetime distribution could be considered for the lifetimes of the susceptible individuals by choosing different parametric forms for f10 (t1), f01 (t2), S10 (t1) and S01 (t2) used in Sections 3.1 and 3.2 in the analysis of lifetimes in presence of cure fraction and censored data.

Figure 2 :
Figure 2: Trace plot and estimated posterior density of the simulated samples considering the parameter λ1 from the FGM copula.

where
the Monte Carlo estimate of DIC (Deviance Information Criterion), introduced by Spiegelhalter et al. (2002), used as a discrimination criterion for different models.Smaller values of DIC indicates better models.The deviance can be expressed as, D(θ) = −2logL(θ | y) + c, L(θ | y) is the likelihood function for the unknown parameters in θ given the observed data y and c is a constant not required for comparing models.Spiegelhalter et al. (2002) defined the DIC criterion by, DIC = D(θ ˆ) + 2nD, (42) where D(θˆ) is the deviance evaluated at the posterior mean θˆ and nD is the effective number of parameters in the model, namely nD = D¯ − D(θˆ), where D¯ = E[D(θ)] is the posterior deviance measuring the quality of the goodness-offit of the current model to the data.

Table 1 :
Posterior summaries (bivariate Weibull distribution derived from FGM copula in presence of cure fraction)

Table 2 :
Posterior summaries (independent Weibull distributions in presence of cure fraction)