Regression Analysis Of Censored Failure Time Data

In recent years,regression analysis of censored data has attracted great attention.Many semiparametric regression models and the corresponding estimating procedures have been proposed in the literature.There are several types of censored data,in this dissertation,we mainly focus on regression analysis of interval censored data and doubly censored data.By interval censored data,we mean that the failure time of interest can not be observed exactly but only known to occur in a time interval(L,R).Such data occur frequently in many fields including medical studies,sociology,demographic studies and tumorigenicity experiments(Zhang et al.,2005;Sun,2006;Wang et al.,2015).In general,there are two types of interval censored data,Case 1 interval censored data and Case 2 interval censored data.By Case 1 interval censored data,we know that each subject is observed only once and the failure time of interest occurs before or after the observation time.In other words,the failure time of interest is either left censored(L = 0)or right censored(R = ?)rather than observed exactly.In the literature,Case 1 interval censored data are often referred to as current status data(Rossini and Tsiatis,1996;MacMahan et al.,2013).By Case 2 interval censored data,we mean that the failure time of interest occurs in a finite time interval,0<L<R<?.Many authors have discussed regression analysis of current status data(Huang,1996;Rossini and Tsiatis,1996;Lin et al.,1998;Martinussen and Scheike,2002;Sun,2006;Chen et al.,2009;Hu et al.;2009;Wen and Chen,2011).All the methods mentioned above are based on the assumption that the failure time is independent of the observation time or censoring time.However,there may exist some dependence between the failure time and the censoring time in current status data,which is often referred to as dependent or informative censoring.For this problem,we propose to use frailty model to characterize the dependence of the two times.Consider a failure time study that consists of n independent subjects and only gives current status data.For subject i,let Ti denote the failure time of interest and suppose that there exists a d-dimensional vector of covariates denoted by Xi.Also for the subject,suppose that there exist two other associated times Ci and CiC related to the observation on Ti,where Ci may be related to Ti and CiC is independent of Ti.For the information on Ti,one only observes Ci = min(Ci,CiC),?i = I(Ci = Ci)and?i=I(Ti ? Ci).That is,the observed current status data have the formIn the following,we assume that the main goal is to estimate the covariate effects.To describe the covariate effects,suppose that there exists a latent variable bi with mean one and variance ? and given Xi and bi,the cumulative hazard function of Ti has the form?(t|Xi,bi)= ?1(t)exp(XiT?1)bi.(1)Here ?1(t)denotes an unknown baseline cumulative hazard function and ?1 is a d-dimensional vector of regression parameters.As Ti,in practice,the observation time Ci may depend on covariates too.For this,we assume that given Xi and bi,Ci follows the same model given by?c(t|Xi,bi)= ?2(t)exP(XiT?2)bi,(2)where ?2(t)and ?2 are defined as ?1(t)and ?1,respectively.In the following,as usual,we will assume that Ti and Ci are conditionally independent given the latent variable bi.Let S(t)= exp{-?1(t)exp(XiT ?1)bi},the survival function of Ti given Xi and bi,Sc(t)= exp{-?2(t)exp(XiT?2)bi} and fc(t)= d?2(t)exp(XiT ?2)bi Sc(t),the survival and density functions of Ci given Xi and bi,respectively.Then under the assumptions above,the likelihood function has the formwherep(·;?)denotes the density function of the bi's assumed to be known up to?.In particular,if the bi's follow the gamma distribution,one can easily show that Ln(?1,?2,?1,?2,?)has closed form,which will be given in Chapter 2.Now we consider estimation of regression parameters ?1 and ?2 as well as others.For this,we propose an EM algorithm by using some Poisson varibles(McMahan et al.,2013;Wang et al.,2015).For the sieve approximation,it will be supposed that the function ?1(t)can be approximated by monotone splines,which can be written asIn the above,{Il(t),l = 1,...,Kn} are integrated spline basis functions,each of which is nondecreasing and ranging from 0 to 1,and the ?l's are nonnegative coefficients that ensure monotonicity of ?1(t).For the function ?2(t),it is apparent that one can approximate it by using the same method above.On the other hand,it will be seen below that one can actually estimate it more easily and directly instead.Define ?=(?1T,?2T,?)T= with ?1=(?1T,?1,...,?Kn)T and ?2 =(?2T,?2)T.We propose to estimate ? by the value that maximizes the likeli-hood function Ln(?).To determine it,note that if the latent variables bi's were known,the likelihood function would have the formDefine a mapping between a new latent variable Zi and ?i such as ?i = I(Zi>0),where Zi follows the Poisson distribution with mean ?i(ci)exp(XiT?1)bi.Then it is easy to show that we can rewrite L1c(?)asif the Zi's were known,where p(Zi)denotes the probability function of Zi and 00 is assumed to be equal to 1.Next note that we can decompose Zi as the sum of Kn independent latent variables Zit's,where Zil follows the Poisson distribution with mean?l Il(Ci)exp(XiT ?1)bi,l = 1,...,Kn.Thus if the Zi's and Zil's were known,we would have the following complete data likelihood functionLet ?(m)denote the estimator of ? obtained in the m.th iteration.To obtain ?(m+1),in the E-step,we need to determine E{lc(?)|?(,m)},which can be written aswhere and Q4(?(m)does not involve ?.The conditional expectations above will be given in Chapter 2.In the M-step,we can calculate ?l's explicitly by setting which givesBy plugging(3)into Qi(?1,?(m)),we can derive the following estimating equation for?1To obtain the updated estimator of ?2,by treating A2 as a piecewise constant function between the uncensored observation times,one can get the score estimating functionwhere,Ni(t)= I(Ci?t,?i=1)andYi(t)= I(Ci?t).Given ?2(m+1),the updated estimator of A2(t)is given by the following breslow-type estimator.In summary,by combining all steps above,we can have the following algorithm.Step 1.Choose an initial estimator ?(0).Step 2.At the(m + 1)th iteration,first calculate the conditional expectations E{bi},E{log(bi)},E(Zi)and E(Zil).Step 3.Obtain ?1(m+1)by solving the score equation(4)and then determine ?l(m+1)with(3)for l = 1,...,Kn.Step 4.Update ?2(m+1)by solving the following score equation,and set ?2(m+1)to be the estimator given in(6).Step 5.Obtain ?(M+1)= argmax?Q3(?,?(m)).Step 6.Repeat Steps 2-5 until the convergence is achieved.In the following,we will establish the asymptotic properties of the estimators proposed above and denote the estimators by ?n,Define the sieve space where B is compact set in R2d+1,with ?c being the longest follow-up time and K1 and K2 some positive constants.Theorem 1 Suppose that regularity conditions given in Chapter 2 hold.Then as n ? ?,we have d(?n,?0)= O(n-(1-v)/2 + n-rv).Theorem 2 Suppose that regularity conditions given in Chapter 2 hold.Then as n ?? we havewhere b is any(2d+1)-dimensional vector with ||b||E?1,g is a function with bounded variation on[0,?c],and E is the semiparametric efficiency bound.For inference about regression parameters,it is apparent that we also need to estimate the asymptotic covariance matrix,say E =(?ij),of ?n(?n-?0).For this,we propose to employ the profile likelihood approach(Zeng et al.,2006;Wen and Chen.,2011).Specifically,let ei be a 2d-dimensional vector with 1 at the ith position and 0 elsewhere and hn is a positive constant with the same order as n-1/2.Then one can estimate ?ij bywhere PLn(/3)= sup?3 Ln(?),the profile likelihood function of ?,with?3 =(?2,?,?1,...,?Kn)T.Note that for fixed ?,one can use the simplified version of the EM algorithm described above to obtain the maximum likelihood estimator of ?3 and thus easily obtain the profile likelihood.Secondly,we discuss regression analysis of multivariate current status data.Con-sider a failure time study that involves n independent subjects and each subject con-tains K possibly correlated failure times of interest.Let Tik denote the failure time of kth event in ith subject and Xik be the corresponding d-dimensional vector of co-variates for Tik.Also suppose that each Tik can not be observed exactly but is known only to be smaller or greater than an observation time denoted by Cik.Here Tik is assumed to be independent of Cik given the covaraites,which is also often referred to as independent censoring or noninformative censoring.Thus,the information we have for Tik only consists of {Cik,?ik = I(Tik-Cik)}.Given Xik and the latent variable bi,we assume the conditional cumulative hazard function of Tik takes the form,where ?fk(t)denotes an unknown increasing baseline cumulative hazard function,? is a d-dimensional vector of regression parameters and Gk is a prespecified function that is also increasing.In addition,we assume that Ti1,...,TiK are assumed conditionally independent given bi and the bi's follow a parametric model with mean one and the density function p(bi|?),where 77 denotes the unknown parameter.Note that we can obtain many commonly used models by choosing different Gk.For example,by letting Gk(x)= x,we can obtain proportional hazards frailty model and it gives proportional odds frailty model when Gk(x)= log(1 + x).Then we can write down the likelihood function in the following formNote that the above transformation function can be derived by Laplace transfor-mation of frailty variable with support[0,?)as the following formwhere ?(t|rk)is the density function of the frailty.For instance,when ?(t|rk)is the density function of a gamma variable with mean 1 and variance rk,we can obtain Gk(x)= log(1+rk x)/rk,the logarithmic transformation function.Kosorok et al.(2004)gives many other instances of frailty-induced transformation functions.Therefore,one can convert the transformation frailty model into the proportional hazards model with two sets of random effects.Then the likelihood function above can be re-expressed asNext we will develop a nonparametric maximum likelihood estimation method to estimate the parameters involved in the model.For each k,let t1k,...,tnkk be the the distinct ordered censoring times,a subset of {cik,i=1,...,n},and assume ?k is a step function with non-negative jump sizes at t<...<tnk,denoted by ?lk for l = 1,...,nk.Let ?=(?,?,?)and ??(?lk:l = 1,...,nklk=1,...,K),then we can rewrite the likelihood function as the following formIn the following,we will describe the derivation of our proposed EM algorithm,which relies on two-stage data augmentation involving Poisson variables.In the first stage,it is nature to assume the latent variables bi's and ?ik'S were known,then the likelihood function would have the formIn the second stage,we define a mapping between ?ik and a new latent variable Zik according to ?ik = I(Zik>0),where Zik = ?tlk?cik Zilk and Zilk is a poisson random variable with the parameter Thus if the Zilk's were known,we would have the following complete data likelihood functionLet ?(m)denote the estimator of ? obtained in the mth iteration.To obtain ?(m+1),in the E-step,we first take conditional expectations with respect to all latent variables in the log-likelihood function lc(?)= log Lc(?).For notational simplicity,we will suppress the conditional arguments in all conditional expectations.The conditional expectations will be given in Chapter 3.In the M-step,we need to maximize the following objective function with respect to ?Setting(?)Q(?,?(m))/(?)?lk = 0,we can update ?lk with the following closed-form expressionBy plugging the estimator above into Q(?,?(m)),we can get the score equation for the regression parameter ?Finally,by setting(?)Q(?,?(m))/(?)?= 0,the estimator of ? can be obtained by solving the following score equation.In summary,by combining all preceding steps,we suggest the following estimating procedure.Step 1.Choose an initial estimator ?(0).Step 2.At the(m+1)th iteration,first calculate the conditional expectations E(?ik bi),E(Zilk)and E{h{bi)} at ?(m).Step 3.Update ?(m+1)by solving the equation(9)with one-step Newton-Raphson method.Step 4.?lk(m+1)can be obtained by the expression(8).Step 5.Calculate ?(m+1)by solving ?i=1n(?)E{log(p(bi|?))}/(?)??0,which can also be achieved by one-step Newton-Raphson method.Step 6.Repeat Steps 2-5 until the convergence is achieved.In the following,we will establish the asymptotic properties of the nonparametric maximum likelihood estimator ?n.For notationai simplicity,we let ? =(?T,?)T,?=(?1,...,?K)and ? =(?T,?T)T.Let ?0,?0,and ?k0 be the true values of ?,? and Ak,respectively.Theorem 3 Under conditions(A1)-(A5)given in Chapter 3,||?n-?0|| ?0 and?k=1K supt?|?1,?2|?kn(t)-?k0(t)| ? 0 in probability,where || · || is the Euclidean norm.Theorem 4 Under conditions(A1)-(A5)given in Chapter 3,we have that d(?n,?0)= Op(n-1/3),where d(?n,00)= ||?n-?0||2 +?k=1K?|?kn(c)-?k0(c)|2fk(c)dc,fk(c)is the density of Ck.Theorem 5 Under conditions(A1)-(A5)given in Chapter 3,?n(?n-?0)?d N(0,I0-1),where I0-1 is the inverse of the efficient Fisher information matrix I0.That is,?n is asymptotically efficient.For inference about the parameters of interest,which contain the regression coeffi-cients ? and frailty parameter 77,it is apparent that we need to estimate the asymptotic covariance matrix,say ?.Since it would be very difficult to derive a consistent esti-mator of ?,for this,we suggest to employ nonparametric bootstrap method(Efron,1981;Su and Wang,2016).Finally,we will discuss regression analysis of doubly censored data.Doubly cen-sored data occur in many areas including demographical studies,epidemiology studies,medical studies and tumorigenicity experiments(Cai and Cheng.,2004;Kim et al.,2013;Zhang and Jamshidian,2004).By doubly censored data,we mean that the fail-ure time of interest can be observed exactly only when the failure event of interest occurs within a certain interval or window,[L,R],and otherwise,it will be known only to be smaller than L or lager than R.In other words,the failure time of interest is either left-censored or right-censored if the event occurs outside the window.One common situation where doubly censoring occurs is that the outcome of the interest can only be measured within a certain or given range.In the following,we will discuss regression analysis of doubly censored data.Consider a failure time study that involves n independent subjects and for subject i,let Ti denote the failure time of interest and Xi(·)be the d-dimensional vector of possibly time-dependent covariates.To describe the possible covariate effects,in the following,we will assume that the cumulative hazard function of Ti has the formwhere G is a known increasing function,?(t)denotes an unknown increasing baseline cumulative hazard function,and ? is a d-dimensional vector of regression parameters.It is easy to see that the model above provides a general class of flexible models that include many commonly used models as special cases.For example,one can obtain the proportional hazards model by letting G(x)= x and it gives the proportional odds model when G(x)= log(1 + x).Also if assuming that the covariates are time-invariant,we will have the linear transformation modellog ?(T)=-XT? + ?,that has been discussed by many authors(Chen et al.,2002;Zhang et al.,2005;Lin et al.,2014;Zhou et al.,2014).In the above,? denotes the random error term with the known distribution function 1-exp{-G(exp(·))}.Suppose that for subject i,the failure time of interest Ti is observed to be left-censored at Li,right-censored at Ri,or known exactly if Ti is between Li and Ri with Li<Ri.That is,we have doubly censored data.Define ?i1=1 if Ti is left-censored and 0 otherwise,?i2=1 if Ti is observed exactly and 0 otherwise,and ?i3=1 if Ti is right-censored and 0 otherwise with the constraint ?i1 + ?i2 + ?i3 = 1 for subject i.We will assume that Li = 0 for a right-censored observation and R,is equal to ? if ?i1= 1.Then the observed data have the form {(Ti,?i1,?i2,?i3,Li,Ri,Xi);i = 1,...,n},where Ti = mar{Li,min(Ti,Ri)}.Let f and S denote the density and survival functions of the failure time,respectively,and assume that the failure time of interest Ti and(Li,Ri)are conditionally independent given covaxiates.Then the likelihood function can be written aswhereandwith G'(x)= dG(x)/dx.Note that one can transform or convert it to the type of proportional hazards frailty models by applying the Laplace transformation(Kosorok et al.2004;Zeng and Lin,2007).More specifically,suppose that there exist some latent variables ?i's with the density function ?(?|r),where r is some known constant.Then we will have the proportional hazards frailty model by settingCorrespondingly,the likelihood function can be rewritten asAmong others,one common choice for ?(?|r)is the density function of a gamma variable with mean 1 and variance r and in this case,we have G(x)= log(1 + rx)/r,the logarithmic transformation function.Now we discuss the estimation of unknown parameters by maximizing the like-lihood function given in(11).For this,we will assume that A(t)is a step function with non-negative jumps only at uncensored failure times and an additional set of left-censored observation times(Mykland and Ren,1996).We will denote these time points by ti<…<tKn with corresponding jump sizes ?(tk),k = 1,...,Kn.Then the likelihood function L becomesFor notational simplicity,in the following,we will denote ? k =?(tk)and Xik=Xi(tk)for k=1,...,Kn.To describe the EM algorithm for maximizing L,for subject i,let {Zik}k=1 Kn denote a set of subject-specific independent latent variables with Zik following the Poisson distribution with mean ?keXikT??i.Let p(Zik|?keXikT??i)be the density function of Zik and ? denote all unknown parameters.Then if assuming that the latent variables ?i's and Zik's were known,we would have the pseudo complete data likelihood functionwhere ?tk?ti Zik>0 if ?i1=1,?tk<ti,Zik=0 and Zik|tk=ti=1 if ?i2=1 and?tk?ti,Zik=0 if ?i3=1.In the E-step of the EM algorithm,we need to calculate the conditional expectation of the log-likelihood function lc(?)= log Lc(?)with respect to all latent variables,which yields Q(?,?(m)= ?i=1n ?k=1Kn XikT?E(Zik)+ log(?k)E(Zik)-?keXikT? E(?i).Note that here for notational simplicity,we have suppressed the conditional arguments in all conditional expectations.In other words,we need to evaluate the conditional expectations given in Chapter 4.In the M-step,by setting(?)Q(?,?(m))/(?)?k = 0,we can update ?k with the following closed-form expressionBy plugging the estimators above into equation Q(?,?(m))we can obtain the following estimating equation for the regression parameter ?In summary,the EM algorithm discussed above can be summarized as follows.Step 1.Choose an initial estimator ?(0).Step 2.At the(m + 1)th iteration,first calculate the conditional expectations E(Zik)and E(?i)by setting ? = ?(m).Step 3.Update ?(m+1)by solving the score equation(13)with one step Newton-Raphson method.Step 4.Update ?k(m+1)by(12)after replacing ?by ?(m+1).Step 5.Repeat Steps 2-4 until the convergence is achieved.In the following,we will discuss the asymptotic properties of the estimators ?n.Let ?0 and AO(·)denote the true values of ? and A(·),respectively.Theorem 5 Assume that conditions(A1)-(A5)given in Chapter 4 hold.Then we have that as n??,||?n-?0||?0 and supt?[0,?1]|?n(t)-?0(t)| ? 0 almost surely,where ||·|| denotes the Euclidean norm.Theorem 6 Assume that conditions(A1)-(A5)given in Chapter 4 hold.Then we have that as n??,the random element ?n(?n-?0)T,?n(·)-?0(·))T converges weakly to a zero-mean Gaussian process in the metric space Rd×l00[[0,?1],where l?[0,?i]denotes the normed space consisting of all the bounded functions with the norm defined to be the supremum norm on[0,?1].Furthermore,?n is asymptotically efficient.Since it would be very difficult to derive a consistent estimator of the asymptotic covariance matrix of ?,for this,we suggest to employ nonparametric bootstrap method(Efron,1981;Su and Wang,2016).