The Log-exponentiated-weibull Regression Models with Cure Rate: Local Influence and Residual Analysis

In this paper the log-exponentiated-Weibull regression model is modified to allow the possibility that long term survivors are present in the data. The modification leads to a log-exponentiated-Weibull regression model with cure rate, encompassing as special cases the log-exponencial regression and log-Weibull regression models with cure rate typically used to model such data. The models attempt to estimate simultaneously the effects of covariates on the acceleration/deceleration of the timing of a given event and the surviving fraction; that is, the proportion of the population for which the event never occurs. Assuming censored data, we consider a classic analysis and Bayesian analysis for the parameters of the proposed model. The normal curvatures of local influence are derived under various perturbation schemes and two deviance-type residuals are proposed to assess departures from the log-exponentiated-Weibull error assumption as well as to detect outlying observations. Finally, a data set from the medical area is analyzed.


Introduction
Models for survival analysis typically assume that all units under study are susceptible to the event and will eventually experience this event if the follow-up is sufficiently long.However, there are situations where a fraction of individuals are not expected to experience the event of interest; that is, those individuals are cured or insusceptible.For example, researchers may be interested in analyzing the recurrence of a disease.Many individuals may never experience a recurrence; therefore, a cured fraction of the population exists.Cure rate models have been applied to estimate the possibility of a cured fraction.These models extend the understanding of time-to-event data by allowing the formulation of more accurate and informative conclusions.These conclusions are otherwise unobtainable from an analysis which fails to account for a cured or insusceptible fraction of the population.If a cured component is not present, the analysis reduces to standard approaches of survival analysis.Cure rate models have been used for modeling time-to-event data for various types of cancer, including breast cancer, non-Hodgkin lymphoma, leukemia, prostate cancer and melanoma.Thus, a cure rate model is suitable for modeling data from cancer clinical trials.Berkson and Gage (1952) introduce the mixture cure rate models, Maller and Zhou (1996) give an extensive discussion of classic methods of inference for mixture cure rate models.An alternative formulation of the parametric cure rate models is discussed in Yakovlev and Tsodikov (1996).A Bayesian formulation of this model is given in Chen et al. (1999).Tsodikov et al. (2003) give an excellent review of such methods.Outside the applications in this area, Hoggart and Griffin (2001) focused on the time to a customer leaving a bank and Yamaguchi (1992) applied the mixture cure rate models to the analysis of permanent employment in Japan.
The literature presents many applications of the survival models with cure rate considering the Weibull family of distributions (see, Ibrahim et al., 2001;Maller and Zhou , 1996).This family is suitable in situations where the failure rate function is constant or monotone.This paper examines the statistical inference aspects and the modeling of the presence of a cure rate of a given data set by using the log-exponentiated-Weibull regression model.The inferential part was carried out using the asymptotic distribution of the maximum likelihood estimators, which in situations when the sample is small, the normality is more difficult to be justified.As an alternative for classic analysis, we explore the use of the Bayesian methods via Markov Chain Monte Carlo.
Influence diagnostics is an important part in the analysis of a data set, as it provides us with an indication of bad model fitting or of influential observations.Cook (1986) proposed a diagnostic approach, named local influence, to assess the effect of small perturbations in the model and/or data on the parameter estimates.Several authors have applied the local influence methodology in regression models more general than the normal regression model (see, for example, Paula 1993, Galea et al., 2000, Díaz-García et al., 2004, and Le et al., 2006).Moreover, some authors have investigated the assessment of local influence in survival analysis models: for instance, Pettitt and Bin Daud (1989) investigated local influence in proportional hazard regression models, Escobar and Meeker (1992) adapted local influence methods to regression analysis with censoring, Ortega et al. (2003) considered the problem of assessing local influence in log-exponentiated-Weibull regression models with censored observations and Ortega et al. (2008) investigated local influence in the Weibull mixture cure models.
In Section 2 we briefly describe the cure rate model.In Section 3 we suggest a log-exponentiated-Weibull regression model with cure fraction, in addition with the maximum likelihood estimators and the Bayes estimator.In Section 4, we discuss the global and local influence method, the likelihood displacement is used to evaluate the influence of observations on the maximum likelihood estimators.We present residual from a fitted model using the martingale residual proposed by Therneau et al. (1990) and we proposed a modified deviance residual for the generalized log-gamma regression model with cure fraction for the log-exponentiated-Weibull regression model with cure rate in Section 5.In Section 6 presents the results of an analysis with a real data set and analysis residual some and conclusions appear in Section 7.

The Cure Rate Model
As in Yakovlev and Tsodikov (1996) and Chen et al. (1999), we formulate the model within a biological context.The promotion time for the jth tumor cell is denoted by R j , j = 1, . . ., N, where N is random variable unobservable that denotes the number of competing causes that can produce a detectable cancer.If N = 0, we define P (R 0 = ∞) = 1 to represent a cure.Hence, the observable time to relapse of cancer or failure time is defined as If N is a Poisson random variable with mean φ independent of the sequence R j , j = 1, 2, . . ., also assumed independent random variable with the same cumulative distribution function (c.d.f.) F (.) and S(.) = 1−F (.), we have that the population survival function is given by where which is a proper survival function and that is a proper probability density function.Suppose we have n subjects and let N i latent variables that denote the number of competing causes that can produce a detectable cancer for the ith subject.Further, we assume that N i 's are independent with Poisson distributions with means φ i , i = 1, . . ., n. Further, suppose R i1 , . . ., R iN i are the promotion times for the N i competing causes for the ith subject, which are unobserved, and all have proper cumulative distribution function, F (•|ψ) and survival function S(•|ψ) = 1 − F (•|ψ), where ψ is the parameter vector.
As in Chen et al. (1999), we assume that the mean of N i is such that, φ(x i ) = exp(x T i γ), where x T i = (x i1 , . . ., x ip 1 ) denotes the p 1 × 1 vector of covariate values for the ith subject, and γ T = (γ 1 , . . ., γ p 1 ) denotes the corresponding vector of regression coefficients.Let t i denote the survival time for subject i, t i = min(T i , C i ), with T i = min(R i0 , R i1 , . . ., R iN i ) and C i is the censoring time whereas δ i is the censoring indicator, assuming 1 if t i = T i and 0 otherwise.The observed data is D = (n, t, δ, X), where t The complete data is given by D c = (n, t, δ, X, N), where N is an unobserved vector of latent variables.Chen et al.(1999) show that the likelihood function for the complete data is Assuming out the unobserved latent variable N in (2.2) the marginal loglikelihood function for the observable data is given by Recently De Castro et al. (2007) show that the likelihood function in (2.3) can be represented by the following expression:

The Log-exponentiated-Weibull Regression Models with Cure Rate
The random variable T nonnegative has a exponentiated Weibull (WE) distribution is its probability density is of form and survival function is given by, where t > 0 and α > 0, θ > 0, are shape parameters and λ > 0 is scale parameter.
It can be shown that the hazard function is given by The great flexibility of this model to fit lifetime data, is given by the different forms that the hazard function (3.3) can take, that is, (i) if α ≥ 1 and αθ ≥ 1, then the hazard function is monotonically increasing; (ii) if α ≤ 1 and αθ ≤ 1, then the hazard function is monotonically decreasing; (iii) if α > 1 and αθ < 1, then the hazard function is bathtub shaped and (iv) if α < 1 and αθ > 1, then we have a unimodal hazard function.In Figure 1, we have some case specials for hazard function (3.3).Applications of the WE distribution in reliability and survival studies were investigated by Mudolkar et al. (1995) and Cancho et al. (1999).Cancho and Bolfarine (2001) proposed the exponentiated-Weibull mixture model to modeling the presence of cure fraction in lifetime data.Some properties of this distribution have been studied more detailed in Mudholkar and Hutson (1996) and Nassar and Eissa (2003).
We assume that the promotion time in the cure rate model in (2.1) follows a exponentiated-Weibull distribution with density function given by (3.1).We consider the following log-exponentiated-Weibull regression models where Y i is log survival time for subject i, β = (β 1 , . . ., β p 2 ) is a vector of unknown parameters to be estimated, z i = (z i1 , z i2 , . . ., z ip 2 ) is the explanatory vector and σ > 0 is unknown parameter and ε i are random variables assumed to be identically distributed with common probability density function where θ > 0 is unknown parameter.It is supposition implicates that the variable Y i has been log-exponentiated-Weibull distribution (see, Cancho et al.,1999).
On the other hard, by considering the model with cure rate in (2.1) and by assuming that such log(R ij ) follow log-exponentiated-Weibull distribution with . This model will be referred to as the log-exponentiated-Weibull regression model with cure rate (LEWR-CR).This model is an extension of an accelerated failure time model using the exponentiated-Weibull distribution and it allows to determine the effect of covariates both on the failure time and on the cure rate itself.The corresponding density function is given by where f pop (y) = − d dy S pop (y).Note that f pop (y) is not a proper probability density function since S pop (y) is not a proper survival function.

Likelihood inference for the LEWR-CR
Considering the result in (2.4) the log-likelihood function for α = ( θ, σ, γ T , β T ) T , corresponding to the observed data is given by where F and C denote, respectively, the set of uncensored and censored individuals, and r is number uncensored observations.This model allows to determine the effect of covariates both on the failure time and on the cure rate itself.The maximum likelihood estimates for the parameter vector α = ( θ, σ, γ T , β T ) T can be obtained by maximizing the log-likelihood function.In this paper, the software Ox (MAXBFGS subroutine) (see Doornik, 2001) was used to compute maximum likelihood estimates (MLE).Covariance estimates for the maximum likelihood estimators α can also be obtained using the Hessian matrix.Confidence intervals and hypothesis testing can be conducted by using the large sample distribution of the MLE which is a normal distribution with the covariance matrix as the inverse of the Fisher information since regularity conditions are satisfied.More specifically, the asymptotic covariance matrix is given by Since it is not possible to compute the Fisher information matrix I(α) due to the censored observations (censoring is random and noninformative), it is possible to use the matrix of second derivatives of the log likelihood, − L(α), evaluated at the MLE α = α, which is a consistent estimator.The required second derivatives are computed numerically.
For testing the adequacy of log-Weibull regression model with cure rate, that is, H 0 : θ = 1, we can consider the likelihood ratio statistics, which in this case is given by where α is the maximum likelihood estimator that follows by maximizing the loglikelihood in (3.8) and α the restricted maximum likelihood estimator computer under H 0 , that is, with θ = 1.For large samples, Λ n is distributed approximately like the chi-square distribution with one degree of freedom.For testing the adequacy of the log-exponential regression model with cure rate, that is, H 0 : (σ, θ) = (1, 1) , the likelihood deviance Λ n is as given in equation (3.9) but α being the restricted maximum likelihood ratio estimator computed under H 0 , that is, with σ = 1 and θ = 1.In this case Λ n is distributed in large samples approximately like the chi-square distribution with two degrees of freedom.

Bayesian inference for the LEWR-CR
The use of the Bayesian method besides being an alternative analysis, allows the incorporation of previous knowledge on the parameters through informative priori densities.When there is not this information one can consider noninformative priors.In the Bayesian approach, the relevant information to the model parameters is obtained through posterior marginal distributions.As such, two difficulties arise.The first refers to the attainment of the required marginal posterior distributions and the second to the calculation of the moments of interest.Both cases require solving integrals that many times do not present analytical solution.In this paper we have used the simulation method of Markov Chain Monte Carlo (MCMC), such as the Gibbs sampler and Metropolis-Hasting algorithm to encompass such difficulties.
We consider the joint prior density for α = ( θ, σ, γ T , β T ) T of the form where denoting the Inverse Gamma distribution with shape parameter a > 0 and scale parameter b > 0 and G(e, f ) denoting the Gamma distribution with mean e/f and variance a/f 2 .We assume that the hyperparameters are specified.
To implement the MCMC methodology, the full conditionals of the parameters are given by Since the above conditional posteriors do not present standard forms, the use of the Metropolis-Hasting sampler is required.

Influence Diagnostics
Local influence calculation can be carried out for model ( 6).If likelihood displacementLD(ω) = 2{l( α) − l( αω )} is used, where αω denotes the MLE under the perturbed model, the normal curvature for α at direction d, d = 1, is given by ] −1 ∆d|, where ∆ is a (p 1 + p 2 + 4) matrix that depends on the perturbation scheme and whose elements are given by ∆ vi = ∂ 2 l(α|ω)/∂θ v ∂ω i , i = 1, 2, . . ., n and v = 1, 2, . . ., p 1 + p 2 + 4 evaluated at α and ω 0 , where ω 0 is the no perturbation vector (see Cook, 1986).For the logexponentiated-Weibull regression models with cure rate, the elements of L(α) are given in appendix A. We can also calculate normal curvatures C d (θ),C d (σ), C d (γ) and C d (β) to perform various index plots, for instance, the index plot of d max , the eigenvector corresponding to C dmax , the largest eigenvalue of the matrix C d i (β) named total local influence (see, for example, Lesaffre and Verbeke, 1998), where d i denotes an n= vector of zeros with one at the i − th position.Thus, the curvature at direction d i assumes the form denotes the ith row of ∆.It is usual to point out those cases such that

Case-weights perturbation
Consider the vector of weights ω = (ω 1 , ω 2 , . . ., ω n ) T .In this case the log-likelihood function takes the form where 0 Then the elements of vector ∆ θ take the form On the other hand, the elements of the vector ∆ σ can be shown to be given by The matrix ∆ γ (p 1 × n) is expressed as The matrix ∆ β (p 2 × n) is expressed as

Response perturbation
We will consider here that each y i is perturbed as y iw = y i + ω i S y , where S y is a scale factor that may be the estimated standard deviation of Y and ω i ∈ R.
Here the perturbed log-likelihood function becomes expressed as where h * i = σ −1 (y i + ω i S y − z T i β).In addition, the elements of the vector ∆ θ take form On the other hand, the elements of the vector ∆ σ can be shown to be given by The entries of the matrix ∆ γ (p 1 × n) can be expressed as Furthermore, the elements the matrix ∆ β , (p 2 × n) can be expressed as

Explanatory variable perturbation (Cure rate)
Consider now an additive perturbation on a particular continuous explanatory variable (rate cure), namely X t , by making x itω = x it + ω i S x , where S x is a scaled factor, ω i ∈ R.This perturbation scheme leads to the following expressions for the log-likelihood function and for the elements of the matrix ∆: In this case the log-likelihood function takes the form In addition, the elements of the vector ∆ θ are expressed as The elements of vector ∆ σ are expressed as The elements the matrix ∆ γ may be expressed when k = t The elements of the vector ∆γ, when k = t are given by The elements the matrix ∆ β can be expressed

Explanatory variable perturbation (Failure Time T )
Consider now an additive perturbation on a particular continuous explanatory variable, namely Z t , by making z itω = z it + ω i S t , where S t is a scaled factor, ω i ∈ R.This perturbation scheme leads to the following expressions for the log-likelihood function and for the elements of the matrix ∆: where h In addition, the elements of the vector ∆ θ are expressed as The elements of vector ∆ σ are expressed as The elements the matrix ∆ β may be expressed when j = t, The elements of the vector ∆ β , when j = t are given by The elements the matrix ∆ γ may be expressed where ĝi = 1 − exp{− exp( ĥi )}, ĥi = , k = 0, 1, . . ., p 1 , j = 0, 1, . . ., p 2 , and i = 1, 2, . . ., n.

Residual Analysis
In order to study departures from the error assumptions as the well as presence of outlying observations, we will consider two kinds of residuals: deviance component residual (see, for instance, McCullagh and Nelder, 1989) and martingaletype residual (see for instance, Barlow and Prentice, 1988;Therneau et al., 1990).More details can be found in Ortega et al. (2003).

Martingale-type and deviance modified components residuals
This martingale-type modified residual was introduced in counting processes and can be written in log-exponentiated-Weibull regression models with rate cure as (5.1) More details about counting processes can be found, for instance, in Fleming and Harrington (1991) and Ortega (2001).These authors show that the distribution of the deviance component residual based on the martingale residual has very close asymptotic distribution to the normal distribution.
Therefore, we have that the deviance component residual for log-exponentiated-Weibull regression models with rate cure becomes (5.2)

Residual Analysis
In order to study departures from the error assumptions as the well as presence of outlying observations, we will consider two kinds of residuals: deviance component residual (see, for instance, McCullagh and Nelder, 1989) and martingaletype residual (see for instance, Barlow and Prentice, 1988;Therneau et al., 1990).More details can be found in Ortega et al. (2003).

Martingale-type and deviance modified components residuals
This martingale-type modified residual was introduced in counting processes and can be written in log-exponentiated-Weibull regression models with rate cure as (5.1) More details about counting processes can be found, for instance, in Fleming and Harrington (1991) and Ortega (2001).These authors show that the distribution of the deviance component residual based on the martingale residual has very close asymptotic distribution to the normal distribution.

Maximum likelihood results
To obtain the maximum likelihood estimates (MLEs) for the parameters in the log-exponentiate-Weibull regression model we use the subroutine MAXBFGS in Ox, whose results are given in the following Table 1.
We note the covariate treatment is significative (at 1%) in the log of time T , the predictor nodule is significative (at 5%) for both in the log of time and cure fraction, also the predictor sex is significative (at 10%).
The value of the likelihood ratio to test the null hypothesis H 01 : θ = 1 in (3.9) provides Λ n = 18.772, which clearly is significante at the 5%, with critical value χ 2 1,0.05 = 3.841.Clearly, the LER-CR model is also not adequate, since for testing H 0 : θ = 1, σ = 1, the observed value of likelihood ratio statistics is Λ n = 49.104(2 degrees of freedom) with p-value ≈ 0. Thus, the likelihood ratio test indicates that the LWER-CR model presents a much better fit that the LWR-CR to the data set under considerations.

Bayesian analysis
We consider now a Bayesian analysis for the data set, using minimal prior information.The following independent priors were considered to perform the Gibbs sampler: β i ∼ N (0, 100), i = 1, . . ., 6 γ i ∼ N (0, 100), i = 0, 1, . . ., 6 σ ∼ IG(1, 0.01), θ ∼ G(1, 0.01), so that we have a vague prior distribution.Considering these prior densities we generated two parallel independent runs of the Gibbs sampler chain with size 50000 for each parameter, disregarding the first 10000 iterations to eliminate the effect of the initial values and to avoid correlation problems, we considered a spacing of size 10, obtaining a sample of size 4000 from each chain.To monitor the convergence of the Gibbs samples we used the between and within sequence information, following the approach developed in Gelman and Rubin (1992) to obtain the potential scale reduction, R. In all cases, these values were close to one, indicating chain convergence.In Table 2 we report posterior summaries for the parameters of the log-exponentiated-Weibull regression model with cure rate.
To compare the LEWR-CR model and LWR-CR model fits by inspecting the estimated of the expected value of Akaike's Information Criterion (EAIC), the expected value of Bayesian Information Criterion (EBIC) and Deviance information criterion (DIC) (see, Spiegelhalter et al., 2002), all put together in Table 3.According to EAIC, EBIC and DIC, the log-exponentiated-Weibull regression (LEWR-CR) model improves the corresponding log-Weibull regression (LWR-CR) model.Similar conclusions can be made analyzing the credibility interval for the parameter θ under LEWR-CR model since it does not contain θ = 1.

Local influence analysis
In this section, we will make an analysis of local influence for the cancer data.

Case-weights perturbation
By applying the local influence theory developed in Section 4, where caseweight perturbation is used, value C dmax (α) = 1.5853,C dmax (γ) = 1.2163 and C dmax (β) = 1.4495 was obtained as maximum curvature.In Figure 3, the graph for the eigenvector corresponding to |d max (α)|, |d max (γ)| and |d max (β)| for all points is presented.Clearly, the most influential is observation 341 on α (See Figure 3(a)).But marginally we noticed that the observation 47 can be influential on γ and the observation 341 on β.We also mentioned that the observation 47 introduce one of the largest lifetimes and the observation 341 presents the smallest lifetime.

Influence using response variable perturbation
Next, the influence of perturbations on the observed survival times will be analyzed.The value for the maximum curvature calculated was C d max (α) = 19.303,C d max (γ) = 1.3153 and C d max (β) = 6.6636.Figure 4, containing the graph for |d max (γ)|, |d max (α)| and |d max (β)| for all points.Results in Figure 4(d) suggest that the observation 176 is the most influential on α.But marginally we noticed that the observation 386 can be influential on γ and the observations 134 and 369 on β.

Influence using explanatory variable perturbation
The perturbation of the covariables age(x 2 ) and tumor(x 6 ) is investigated here.After perturbation of covariable age, value

Residual analysis
In order to detect possible outlying observations as well as departures from the assumptions of the log-exponentiated-Weibull regression models with rate cure, we present in Figure 7, the graphs of r Di against the order observations.By analyzing the residuals deviance graph, a random behavior is observed for the datas.As we can observe through the local influence analysis and analysis and the residual analysis it doesn't exist jointly influential points.Thus, the final model becomes the one given by y i = β 1 x i1 + β 3 x i3 + β 4 x i4 + σε i , φ(x i ) = exp{γ 0 + γ 1 x i1 + γ 3 x i3 + γ 4 x i4 } (6.1)The parameters MLEs are reported in the Table 4.We note the covariate treatment and nodule is significative (at 5%) in the log of time T , the predictor nodule and sexo is significative (at 5%) in the cure fraction.The covariate nodule is significative for both in the log of time and cure fraction, also the predictor sex and treatment is significative (at 10%) in the log of time and cure fraction respectively.
We may interpret the estimated coefficients of the final model as following.The predictor nodule decelerates the lifetime of individual's and alters the cured proportion significantly, and that significant difference exists among the levels of treatments (observation and interferon) in relation log of time.

Concluding Remarks
In this study, the log-exponentiated-Weibull regression (LEWR) model was modified in order to include long-term individuals.In the proposal under consideration, log-linear parametric modeling was taken as a basis for survival time.The model attempts to estimate simultaneously the covariates effects on the accelaration/decelaration of the timing of a given event and surviving fraction, that is, the proportion of the population for which the event never occurs.
Continuing with modeling investigation, we applied local influence theory (Cook (1986) and Thomas and Cook (1990)) and conducted a study based on martingale and deviance residuals in a survival model with a cure fraction.The necessary matrices for application of the technique were obtained by taking into account various types of perturbations to the data elements and to the model.By applying such results to a data set, indication was found of which observations or set of observations would sensitively influence the analysis results.This fact is illustrated in application (Section 6).By means of a real data set, it was observed that, for some perturbation schemes, the presence of certain observations could considerably change the levels of significance of certain variables.The results of the applications indicate that the use of the local influence technique in models with a cure fraction may be rather useful in the detection of possibly influential points by admitting two types of estimation methodology: maximum likelihood and Bayesian.In order to measure the quality of fitting, martingale and deviance residuals were used, which showed that the model fitting was correct.Finally we can observe that the log-exponentiated-Weibull regression models with rate considered is a robust model.

Figure 2 :
Figure 2: Plot of the Survivor Function for the melanoma data Firstly, we consider the following regression model:

Figure 3 :Figure 4 :
Figure 3: Case-weights perturbation.(a) Index plot of d max for α.(b) Index plot of d max for γ.(c) Index plot of d max for β.(d) (e) (f)

Figure 5 :Figure 6 :
Figure 5: Age explanatory variable perturbation.(g) Index plot of d max for α.(h) Index plot of d max for γ. (i) Index plot of d max for β.

Figure 7 :
Figure 7: Index plot of the deviance

Table 1 :
Maximum likelihood estimates for the log-exponentiated-Weibull regression model with cure rate

Table 2 :
Posterior summaries for the log-exponentiated-Weibull regression model with cure rate.

Table 4 :
Maximum likelihood estimates for the log-exponentiated-Weibull regression model with cure rate