HM Medical Clinic

 

Jirss.irstat.ir

JIRSS (2013)
Vol. 12, No. 1, pp 71-112
Modeling and Inferential Thoughts for
Consecutive Gap Times Observed with Death
e de Strasbourg, Strasbourg, France.
Abstract. In the perspective of biomedical applications, consider a re-
current event situation with a relatively low degree of recurrence. In this
setting, the focus is placed on successive inter-event gap times which are
observed in the presence of both a terminal event like death and inde-
pendent censoring. The terminal event is potentially related to recurrent
events while the censoring process is an independent nuisance that bears
on the total observation time i.e. on the sum of the successive gap times.
We review different modeling and inferential strategies. We also present
a nonparametric estimation method of joint distribution functions and
outline the need for future developments.
Keywords. Consecutive gap times, independent censoring, joint mod-
eling, recurrent events, terminal event.
MSC: 62G05, 62N01, 62N02, 62P10.
Many clinical and epidemiological cohort studies involve health out-comes that a participant may experience a few times during the follow-up period. Interest is particularly centered on life-threatening recurrentevents that each patient may experience at most a very few times. Clas- Received: June 23, 2011; Accepted: July 17, 2012 sical examples of such recurrent events from longitudinal studies includesolid tumor recurrences in cancer patients, transient ischaemic attacksin atherosclerotic patients and recurrent leukaemia in patients undergo-ing allogeneic haemopoietic stem cell transplantation. The occurrenceof recurrent events as serious as tumor recurrence or ischaemic attack isassociated with a high risk of death so that the subject may die duringthe study. Obviously, the death of a patient stops any subsequent occur-rence of recurrent events, hence the name "terminal event". In addition,during a trial, right-censoring phenomenon is common and should alsobe accounted for. Studies such as CARE (Pfeffer et al., 1996), ASPIRE(Bowker et al., 1996), CAPRIE (Gent et al., 1996), LIPID (The LIPIDStudy Group, 1995, 1998 and Marschner et al., 2001) fall within thisframework.
For a given patient, the occurrence of a recurrent event often im- pacts the risk of novel recurrent event and even of death. Therefore theassumption of independence among the gap times of an individual isoften violated for recurrent event data and the death time is also likelyto be dependent on the recurrent event history. This dependence shouldbe taken care of in the inferential procedures and accounted for in thejoint modeling of recurrent and terminal events. On the opposite, whenno covariates are available, the censoring process is often assumed to beindependent of both the recurrent and terminal processes. Possible rea-sons for such an independent censoring are loss to follow-up (non-relatedto side-effects of the treatment under study) or end of study. When co-variates are available, it is often assumed that the censoring process isindependent of both the recurrent and terminal processes conditionallyon covariates.
Throughout, the recurrent and the terminal processes are all assumed to have distribution such that recurrent event and death cannot happenat the same time.
The kind of data collected during such a trial is illustrated on Fig- ure 1 for six subjects labeled S1 to S6. The follow-up time for a givenpatient is represented by a straight line along which the different eventsare indicated. Recurrence time data can be regarded as multivariatedata that have specific characteristics: • the different recurrence times of a given subject are stochastically • the different patients do not experience the same number of events and the number of observed events for each subject is not knownin advance, Modeling and Inferential Thoughts for Consecutive Gap Times • the number of patients still alive and still under study decreases as the events occur, • the terminal event stops the further occurrence of recurrent events, • the last recurrence time of a subject is either a censoring time or • censoring occurs at most once for a given patient and prevents any further event from being observed.
Figure 1: Example of data (RE = recurrent event).
Traditional statistical methods for the analysis of cohort study have been focused either on the survival time or on the first occurrence ofa composite outcome due to lack of appropriate methodology. For in-stance, the primary pre-specified endpoint in the LIPID (1995) studywas coronary heart disease related death and for secondary analyzes acomposite of coronary heart disease related death or non-fatal myocar-dial infarction.
Focusing only on the most serious issue i.e. on the fatal one and comparing treatments only with respect to total lifetime would lead toefficacy problems. Coping with these efficacy problems would requirelonger trial duration and larger sample size. It would also result in aconsiderable loss of information that would obscure the issue of seriousnon-fatal event recurrence which is also a major concern even if the re-currence degree is relatively low. This was a serious matter in Jokhadar et al. (2004) who outlined that, as myocardial infarction hospital fatal-ities decline, survivors are candidates for recurrent events. This paperalso states that the question of morbidity after non-fatal myocardialinfarction and how it may have changed over time with the arrival ofcontemporary treatments is still of interest. Assessment of prognostic forfurther recurrence is a critical step in evaluating the need for treatmentand lifestyle modifications to manage the risk of future events.
On the other hand, the most common analysis consists of restrict- ing attention to first event occurrence and of focusing on a predefinedcomposite endpoint that combine fatal and non-fatal events to demon-strate treatment efficacy. The first consequence of this is bias problemstoward shorter lifetimes. The second consequence is more subtle andhas been outlined in a number of paper including Pocock (1997), Mah´ and Chevret (1999), Ferreira-Gonz ales et al. (2007) and Kleist (2007) among many others. All outlined the fact that composite endpointsshould be clinically meaningful and that the expected effects on eachendpoint component should be similar, based on biological plausibility.
All components of the composite endpoint need to be analyzed sepa-rately. Difficulties in interpretation then arise when the results on singlecomponents of the composite endpoint go in opposite directions andwhen hard clinical outcomes are combined with soft endpoints, particu-larly if the latter occur more frequently but are of inferior relevance.
As a conclusion to this discussion, patients need to be followed up on assigned treatment until death or end of planned follow-up in the absenceof events and must not be regarded as 'a trial completer' after occurrenceof the first component event. More specific regulatory guidelines, betterreporting standards and appropriate statistical methodology are neededto this aim. In this set-up, some progress are still to be made to supportclinical decision making.
Various modeling approaches have been considered with recurrent event data to address different types of questions. The appropriatenessof a selected model depends on the nature of the recurrent event dataas well as on the interest of the study. In the analysis of recurrent eventdata, the focus can be placed either on the times to recurrent event oron the gap times between successive events or on the recurrent eventprocess N ∗(.) where, for t ≥ 0, the process N ∗(t) records the number ofrecurrent events occurring in the time interval [0, t].
Note that a connection can be made with multi-state models by considering the multi-state model depicted in Figure 2 with boxes rep-resenting the states and arrows the possible transitions. In this model, Modeling and Inferential Thoughts for Consecutive Gap Times the first (k0 + 1) states represent the cumulative number of events expe-rienced, the last state is absorbing and represents death. Only forwardtransitions are possible. The "gap time" timescale can be linked to time-homogeneous semi-Markov models in which the transition probabilitybetween two states only depends on the gap time whereas the "time-to-event" scale can be linked to time-inhomogeneous Markov models inwhich the transition probability between two states only depends on thetime since inclusion in the study. Dealing in full details with multistatemodels is beyond the scope of this paper. Nevertheless, we refer tothe recent papers by Andersen and Pohar Perme (2008) and by Meira-Machado et al. (2009).
Figure 2: Multistate model view (RE = recurrent event).
From now on, interest is specifically focused on successive gap times since this approach is rather suited to studying recurrent events dynam-ics. Modeling and analyzing the waiting times between successive eventsis attractive in specific settings. First of all, analyzes based on waitingtimes are often useful when the recurrence degree is relatively low. More-over, when evaluating the efficacy of a treatment on a life-threateningillness where only a few recurrences are expected, it is important to as-sess whether or not the treatment delays the time from treatment startto the first episode, the time from the first episode to the second episodeand so on. Indeed, several phenomenons may affect the understandingof the illness mechanism. On the one hand, a treatment which delays thefirst episode will inevitably lengthen the total time to the second episodeeven if it becomes ineffective after the first episode. On the other hand,in some cases, a compensating phenomenon between the different stagesof a disease or of a treatment may exist. For example, a treatment maydelay the occurrence of the first recurrent event but have the reverseeffect on the occurrence of the subsequent recurrent events. It is impor-tant to detect such a phenomenon. The distribution of successive gaptimes between recurrent events is then a valuable information.
At last, joint modeling of successive gap times including a possible fatal event and independent censoring is needed for practical purposesas stated in the paper of Cui et al. (2010). This problem is currentlynot fully addressed, up to our knowledge.
In Section 2, different modeling strategies frequently used in the literature for gap times inference are exposed. In Section 3, we inves-tigate relevant cause-specific distribution functions. Some perspectivesare given in Section 4.
A Review of Modeling Strategies for Gap
Times
In the literature, various modeling approaches have been adopted for theanalysis of recurrent event data. Amongst these are renewal processes,frailty models or multi-state models. Regression models that incorpo-rate either the past event history as covariates or explanatory covariatesor both are also of much use. At last, a purely non-parametric approachwith as few assumptions as possible regarding either the gap times dis-tribution or association structure is also of much interest.
The literature about gap times distribution function or hazard func- tion inference may be broadly classified into three categories accordingto whether authors focused on univariate functions, joint functions orconditional-on-past-event-history functions. These three analyzes canalso be carried out conditional on explanatory covariates if some areavailable. As we will see in the sequel, the two dominant methods to in-corporate explanatory covariates are Cox's proportional hazards modeland the accelerated failure time model. We now review some importantcontributions.
In the discussion to come, let us denote by Y [k] for k = 1, 2, . the successive gap times between successive recurrent events and by D the
death time if considered. We adopt the convention that Y [0] = 0 and that
Z is a vector of explanatory covariates when available. Time-dependent
vectors of covariates are written as Z(.). Throughout, the use of y is
dedicated to the gap timescale while the use of t is dedicated to the
calendar timescale.
Models for Univariate Functions
Models for univariate functions are useful when interest lies in under-standing the evolution through time of separate gap times even though Modeling and Inferential Thoughts for Consecutive Gap Times possible association is accounted for by some authors. In the absence of(within-subject) gap time independence, two major statistical issues inthe marginal analysis of gap times are identifiability and induced depen-dent censoring. On the one hand, since the study duration is typicallyless than the support of the first failure time, the marginal distributionof, say, the second gap time is not identifiable unless (within-subject)successive gap times are independent as discussed for instance by Wangand Wells (1998), Wang (1999), Lin et al. (1999). This explains theneed for more or less strong modeling assumptions even though marginalmethods are often said to be robust to the subject-specific correlationstructure between gap times. On the other hand, even if censoring actsindependently on the first gap time and on the total observation time,the second and subsequent gap times are subject to induced depen-dent censoring. For example, a greater first event time implies highercensoring probabilities for the second and subsequent gap times sinceindependent censoring bears on the times-to-event ie on the sums of thesuccessive gap times. Failures to account for this association may leadto substantial bias in dealing with gap times after the first.
We now review different approaches to univariate modeling.
Renewal processes have the property that the gaps between succes- sive events are independent and identically distributed. Even thoughmuch less suited to biomedical applications, one can nonetheless men-tion the work of Pe˜ na et al. (2001) under this renewal process assump- tion. The authors established Nelson-Aalen-like and Kaplan-Meier-likeestimators of the gap times marginal cumulative hazard function and dis-tribution function using the method of moments and then argued thatthey are NPMLE respectively for the gap times marginal cumulativehazard function and the marginal distribution function. The authorsshowed that any deviation from the independence assumption provokesbias problems that logically increase as the level of association betweengap times increases.
The assumption that gap times are independent and identically dis- tributed is very strong when no covariates are present. Therefore it isimportant to consider diagnostic checks. An important way of modelchecking is by fitting models that include renewal processes as specialcase. The independence assumption can also be checked informally whenno covariates are present by, for example, looking at scatter plots of suc-cessive gap times within individuals. There should be an absence oftrend if the renewal assumption is valid. If the gap times are indepen-dent, then informal checks on the assumption of a common distribution can also be made by comparing separate empirical distributions for thedifferent gaps. Anyway, the independence assumption is a very strongcondition and renewal processes are mainly useful in reliability when asubject is repaired in some sense after each event. The renewal assump-tion is clearly untenable in most biomedical applications.
Specific extensions of renewal processes can be found in Cook and Law-less (2008) and in Gill and Keiding (2010).
A classical generalization of renewal models that allows association between gap times consists of considering a frailty model in which a la-tent variable is used to take into consideration a subject-specific randomeffect. Specifically, it is supposed that, for any given subject, there existsan unobserved random variable U , called the frailty, with distributionFU such that given U = u the successive gap times of the individual arei.i.d. with some distribution F (. u). This unmeasured effect is assumedto follow a distribution with mean equal to one and unknown finite vari-ance. A distribution such as the Gamma (most popular frailty model)or Inverse Gaussian or positive stable or log-normal can be assumed forthe frailty. The goal of the frailty model analysis is generally to esti- mate the distribution function F (.) = F (. u)dFU (u) unconditional on the frailty or the hazard function also unconditional on frailty. Tech-nically speaking, frailty models can be fitted either with a frequentistapproach by maximizing the marginal likelihood or with a Bayesian ap-proach by computing parameter posteriors densities. Hougaard (2000)and Duchateau and Janssen (2008) provided a comprehensive coverageof this area.
Wang and Chang (1999) focused on a marginal approach for the gap times between successive events using this less restrictive frailtyapproach. They derived a weighted moment estimator for the marginalgap times distribution under a nonparametric frailty model assumingthat, given the frailty, the successive duration times are independentand identically distributed. Their correlation structure is quite generaland contains both the i.i.d. and multiplicative (hence Gamma) frailtymodel as special cases.
na et al. (2001) also proposed an estimator when the gap times follow a Gamma frailty model and compared the performance of theirestimator to Wang and Chang's (1999) estimator by simulations. Theauthors found out that, when applied to i.i.d. gap times, their estimatoris expected to be more efficient than that of Wang and Chang (1999).
From a practical viewopint, it is interesting to note that both Pe˜ et al. (2001) and Wang and Chang (1999) estimators are implemented Modeling and Inferential Thoughts for Consecutive Gap Times in the R package survrec.
One limitation of many types of random effect models is that only one or two parameters are used to model association among a numberof successive gap times. This can be inadequate when association struc-tures are complex or changing over time. Moreover these approaches donot readily deal with negative associations. At last, it should be notedthat misspecification of the frailty distribution can cause severe bias inestimation procedures with recurrent event data, see e.g. Kessing et al.
(1998). This remark entails that the possibility of testing model ade-quacy is an important issue that should be dealt with. Kvist et al. (2007)developed a procedure for checking the adequacy of Gamma frailty torecurrent events. To apply their model checking procedure, a consis-tent non-parametric estimator for the marginal gap time distributionsis needed. The performance of the model checking procedure dependsheavily on this estimator. Moreover, the authors concluded that theirprocedure in its current state only works when the within-subject asso-ciation between gap times is weak. They suggested possible future im-provement of their methods consisting of checking of the Gamma frailtymodel for recurrent events from a comparison of conditional distribu-tions instead of marginal distributions. Up to now, there is still a roomfor improvements on that issue.
When covariates are available, marginal methods are also helpful to understand how either population-level characteristics (e.g. treatmentgroup) or subject-specific (e.g. sex) or gap-time specific characteristics(e.g. some biological marker) influence the marginal gap time distribu-tion. Explanatory covariates are often incorporated through Cox-like oraccelerated-failure-time-like assumptions for their ease of interpretation.
Assuming that the successive gap times of each individual are i.i.d.
(unconditional on covariates), Huang and Chen (2003) considered a pro-portional hazards assumption of the form Λ(y Z) = Λ0(y) exp(β′.Z)
to assess the effect of a vector Z of time-independent covariates on the
common baseline cumulative hazard function Λ0(.) of the successive gap
times. The authors developed an inferential procedure that improves the
functional formulation of Cox regression by Huang and Wang (2000)
with respect to efficiency. To this aim, the authors noticed that the
uncensored gap times are exchangeable provided the model assumptions
are valid, then constructed specific clustered data. For each cluster, the
first gap time is chosen if the subject has only one censored gap time,
otherwise all the uncensored gap times are selected. Then they got theirnew estimates of β and Λ0(.) through a modified estimating equationobtained from the clustered data. This procedure is shown to performwell for practical sample sizes. The authors also noted that their modeland inferential procedure still apply for covariates that depend on timefrom earlier episode and have uniform effects across all gap times butthat difficulties arise if time-varying covariates are episode-specific.
Obviously, the validity of statistical inference depends on the ad- equacy of the model.
Recent progress have been made in Cox-type model checking for gap times in Huang et al. (2010) who proposed bothgraphical techniques and formal tests for checking the Cox model withrecurrent gap time data to assess different aspects of goodness-of-fit forthis model.
In the same setting, ie also assuming that the successive gap times of each individual are i.i.d. (unconditional on covariates), Sun et al. (2006)considered an alternative model under the form of an additive hazardsmodel defined as λ(y Z) = λ0(y) + β′.Z
where λ0(.) is the common baseline instantaneous hazard function of thesuccessive gap times. The authors used the same inferential procedureas in Huang and Chen (2003) to ensure satisfying efficiency properties.
Strawderman (2005) also proposed a marginal regression model for consecutive gap times of the accelerated failure type but alleviated the
i.i.d. property of gap times. Specifically, he assumed that, conditional
on a vector of covariates Z, the variables Y [k] exp(β′.Z) for k = 1, 2, .
are i.i.d. or equivalently that the common hazard function of Y [k] con-
ditional on Z is of the form
where λ0(.) is an unspecified baseline function. Similarly to the accel-erated failure time model, explanatory population-level (ie not episode-specific) covariates serve to accelerate or decelerate a baseline gap timehazard function. The problem of obtaining an efficient estimation of βis investigated.
However, the same restrictions and pitfalls as previously also apply to regression models in which the gap times are independent conditionalon covariates. Here again, the conditional independence is questionable.
We also mention that questions exist concerning the interpretation ofbaseline functions when there is association between the successive gaptimes.
Modeling and Inferential Thoughts for Consecutive Gap Times A possible extension consists of incorporating episode-specific covari- Chen, Wang and Huang (2004) considered a situation in which episode- specific vectors of covariates, say Z[k] for k = 1, 2, . are available and
assumed that the gap times are i.i.d. conditional on covariates. Account-
ing for right-truncation phenomenon in the observation of successive gap
times, they worked with a subject-specific reverse-time hazards function
defined for subject i by
(y Z ) =
P y − h ≤ Y ≤ y, Z
h→0+ h[ d log P Y ≤ y Z
Their approach relies on modeling proportional reverse-time hazardsfunctions so that, for individual i, we have (y Z ) = κ
i.e. each individual is assumed to have its own baseline reverse-timehazard function κi,0(.). Thus the model copes with high heterogeneityacross the patient population. The prize to pay for this generality is po-tential identifiability problems. The authors suggested as a special casethat the baseline reverse-time hazard function could be modeled using afrailty Ui for subject i as κi,0(y) = Uiκ0(y). Note that the effect of theepisode-specific covariates on the baseline reverse-time hazard functionis constrained to be identical across recurrences. The interpretation ofsuch a constraint may be tricky but has the advantage of allowing moreefficient estimation of β.
Similarly, Du (2009) assumed that each gap time Y [k] depends only on an episode-specific covariate Z[k] such that the Y [k] conditionally on
the Z[k] are i.i.d. Stating that the history of a subject before each recur-
rence conveys information for that recurrence, Du (2009) suggested to
include the number of past recurrences in the episode-specific covariates.
The author investigated a nonparametric estimator for the marginal gap
time hazard function of Y [k] conditional on Z[k] = z, denoted by λ(. ),
using a functional ANOVA decomposition of log λ of the form
log λ(y z) = η0 + ηgap(y) + ηcov(z) + ηgap,cov(y, z)
where η0 represents the grand mean, ηgap(.) represents the main gap
time effect, ηcov(z) represents the main covariates effect and ηgap,cov(y, z)
represents the interaction effect between gap time and covariates, pro-vided some identifiability conditions are ensured. As a consequence, thismodel can be applied to assess the validity of the proportional hazardsassumption by examining the interaction between gap times and covari-ates. The inferential procedure is based on non-parametric penalizedlikelihood with a cross-validation step to select smoothing parameter.
This model has the advantage of allowing greater flexibility for the func-tional form of the different effects and of generalizing multiplicative formof the hazard at the expense of loosing the usual interpretation in termsof risk ratio.
An alternative approach to these models consists of accounting for as- sociation between within subject gap times via random effects to lightenthe i.i.d. assumption on successive gap times, conditionally on covariatesin the regression setting.
Therneau and Grambsch (2000) considered that the successive gap times are i.i.d. conditionally on both an (unobservable) frailty U and
an (observed) episode-specific vector of covariates Z[k]. The authors
discussed the fitting of a model for the conditional hazard function of
Y [k] of the form
λ[k](y Z[k], U ) = U λ (y) exp(β′
(.) is an episode-specific baseline hazard function and where βk is the episode-specific effect on the episode-specific baseline hazard
function of the episode-specific vector of covariates Z[k]. Note that here
the Y [k] are conditionally independent but not identically distributed.
This enlarged flexibility may lead to inconsistency problems as k grows
if too few patients experience k events. Except for the case where U has
a positive stable distribution, these models do not give unconditional
(on U ) distributions for Y [k] given Z[k] of proportional hazards form.
Chang (2004) considered a marginal accelerated failure time frailty model of the form log Y [k] = U + β′.Z + ε[k], k = 1, 2, .
where the variables ε[k] for k = 1, 2, . are i.i.d. This model assumes thatthe covariates effect and the subject-specific frailty U are additive on thegap time logarithm and that the covariates effect remains the same overdistinct episodes. The distributions of the frailty and the random errorin the model are left unspecified which decreases adequacy issues. Theauthor developed two estimation methods, the second of which beingrobust to deviation from the hypothesis that the ε[k] for k = 1, 2, . are Modeling and Inferential Thoughts for Consecutive Gap Times identically distributed. The authors mentioned the possibility to extendtheir model to allow the incorporation of an episode-specific covariateeffect log Y [k] = U + β[k]′ .Z + ε[k], k = 1, 2, .
even though consistent estimation of β[k] may not be possible if thenumber of subjects experiencing k events is not large enough.
All these methods, however, assume that recurrent events are not terminated by death during the study.
Rondeau et al. (2007) accounted for death in their analysis. Specif- ically, they jointly modeled the association between survival time and
within subject gap times through a Gamma frailty U . Conditional on
U and on an external time-dependent vector of covariates Z(.), they as-
sumed that the Y [k] for k = 1, 2, . are i.i.d. with conditional hazard
function given by
λ(y U, Z(.)) = U λ0(y) exp(β′.Z(t))
and are independent of the death time with conditional hazard functiongiven by λD(t U, Z(.)) = U αλ0,D(t) exp(γ′.Z(t)).
The frailty effect on recurrent events and death is different unless α = 1.
When α > 1, the recurrent rate ad the death rate are positively asso-ciated since higher frailty results in both higher risk of recurrence andhigher risk of death. The authors proposed a semiparametric penalizedlikelihood estimation method in which the model degree of freedom isused to specify the smoothing parameter. Their method yields unbiasedand efficient estimates. It is noteworthy to say that the work of Rondeauet al. is implemented in a very complete R package named frailtypack.
Huang and Liu (2007) considered a similar situation. Conditional on a Gamma frailty U , on a baseline vector of covariates associated with
survival ZD and on an episode-specific vector of covariates Z[k], they
assumed that the Y [k] for k = 1, 2, . are independent (but not identically
distributed) with Y [k] having conditional hazard function given by
λk(y U, Z[k]) = U λ (y) exp(β[k]′.Z[k])
and are independent of the death time with conditional hazard functiongiven by λD(t U, ZD) = U αλ0,D(t) exp(γ′.ZD).
The authors mentioned the fact that if covariate effects are believed to behomogeneous across gap times in some appropriate practical situation, a common β may be use instead of the β[k] for k = 1, 2, . in order to gainefficiency. However, the additional flexibility of this model with respectto that of Rondeau et al. (2007) induced by episode-specific baselinehazard functions and episode-specific covariate effects may be of limitedpractical use if data are sparse as k grows.
The standard assumption that the frailty U is fixed over time and independent of observed covariates is still strong. With this respect,
Du et al. (2011) proposed a more general model even though they do
not account for the possibility of associated death. Conditional on an
unobserved vector of random effects U and on two vectors of observed
covariates Z and e
Z, Du et al. (2011) assumed that the Y [k] for k = 1, 2, .
are independent with conditional hazard function satisfying log λ(y U, Z, e
Z) = η(y, Z) + e
The vector of covariates Z is expected to impact the gap time dis-
tribution while the vector of covariates e
Z is expected to impact the
random effects. The usual frailty model correspond to e Z′.U = U for
a scalar-valued random variable U that is both time-independent and
covariate-independent. The random effects multivariate distribution is
left completely unspecified which allows to incorporate time-varying
frailty. Moreover, the very general form of the hazard function gives
the possibility to investigate a general shape of the conditional hazard
function and extract useful information that might be missed by para-
metric or semiparametric models. Inference is carried out by iteratively
minimizing a penalized likelihood in which the smoothing parameter
selection is reported as potentially challenging. Extension of episode-
specific covariates Z[k] is claimed to be straightforward.
As as summary of this subsection, in the current state of the lit- erature, a balance has to be made between models relying on strongassumptions that are more or less hard to check and, on the other hand,more general flexible models in which one may have to face identifiabilityand efficiency problems.
As a second modeling strategy, it is possible to estimate meaningfuland identifiable conditional distributions related to the gap times in thepresence of non-informative censoring. Such models usually specify howthe probability (or hazard) function of subsequent recurrence depends Modeling and Inferential Thoughts for Consecutive Gap Times on the past event history which may not be a trivial task. Typically lessrobust to bad specification of subject-level correlation structure betweenevents, these models are useful for studying local process dynamics andpredicting recurrence experience at the subject level.
Amongst the primary attempt to estimate P[Y [2] ≤ y2 Y [1] ≤ y1] is the proposal by Lin et al. (1999) who investigated the estimate definedas the ratio of the estimate of the joint distribution function of (Y1, Y2)(see also next subsection) over the estimate of the marginal distributionof Y [1]. Their inferential procedure is based on the inverse probability ofcensoring weights which is a well-known and useful tool for adjusting theinduced dependent censoring when analyzing multiple gap times betweenrecurrent events. Adjusting for induced dependent censoring consistsof weighting risk set contributions by the inverse of the probability ofremaining uncensored. The estimate is obtained without any modelingassumptions regarding the dependence structure of the successive gaptimes. The standard errors are also derived.
Quite similarly, Schaubel and Cai (2004) proposed an estimator of the conditional survival function for the k-th gap time conditional onthe (k − 1)-th event occurring prior to some fixed time point. Theirwork shares with Lin et al. (1999) the fact that the estimate is obtainedwithout any modeling assumptions regarding the dependence structureof the successive gap times. However, instead of being based on a ratio ofestimate, Schaubel and Cai (2004) proposed an estimator that is deriveddirectly from a cumulative hazard function. From a technical viewpoint,their estimator is not subject to negative mass which is a problem thatmay arise with an estimate that depends on the joint distribution ofthe successive gap times. Another advantage of the proposed techniquesis the ease of computing standard errors which may be important topractitioners. A method for computing simultaneous confidence bandsis also provided.
Regression methods are also available to incorporate covariates into the analysis of conditional distributions.
Chang and Wang (1999) focused on semi-parametric regression for conditional gap times analysis using a Cox model incorporating time-dependent covariates and in which the number of past episodes servesas a stratification variable. Two types of time-dependent covariates areincluded. The first type of covariates has an effect which is expected toremain constant through the distinct episodes while the second kind of covariates effect is episode-specific. Setting k (y) = Z(u) =
: 0 ≤ u ≤ Y [j] + y λ[k](y Y [1], ., Y [k−1], Zk(y)) = P y ≤ Y [k] ≤ y + hY [k] ≥ y, Y [1], ., Y [k−1], Zk(y) , h→0+ h their model can be written as follows for y ≥ 0 λ[k](y Y [1], ., Y [k−1], Zk(y)) = (y) exp β′.Z 
Y [j] + y + γ′ Y [j] + y . k .Z2
Implicitly, the prior event history is summed up by the time-dependentcovariates. To estimate the parameter β, a profile likelihood approachbased on all of the data is adopted to handle the nuisance parameters γk.
In the data, because the number of subjects who experience at least krecurrent events decreases as k increases, a limitation in the estimationof the γk is the lack of sufficient data for consistent estimation when khas large values, as already mentioned elsewhere. However, the authorspoint out that, with appropriate conditions, the regression coefficient βcan be consistently estimated regardless of whether the parameters γkcan or cannot be.
Lawless et al. (2001) reviewed conditional regression models applied to shunt failure data. The following model λ[k](y Y [1], ., Y [k−1], Z[k]) = λ (y) exp β′
k .(Y [1], ., Y [k−1])′ + α′k .Z[k]
which gives a symmetric role to episode-specific covariates Z[k] and past
gap times was considered with an emphasis on the condition βk =
(0, ., 0, bk)′ so that their model incorporates first-order dependence.
The main difference with the model of Chang and Wang (1999) is
that the conditional hazard function now explicitly depends on previ-
ous event. Besides the fact that the functional relationship between the
gap times should be adequate, a potential drawback to the conditional
approach is that the parameters have to be interpreted conditionally to
previous event times.
Modeling and Inferential Thoughts for Consecutive Gap Times Schaubel and Cai (2004b) considered estimation via semi-parametric Cox regression models for conditional gap times hazard functions. Re-specting the identifiability issues, the authors focused on the followinggap-time-specific hazard functions P y ≤ Y [k] ≤ y + hY [k] ≥ y, Y [j] ≤ t h→0+ h for some pre-specified tk−1 chosen in the support of the total observa-
tion time distribution and for some external time-dependent covariates
Z[k](.). They assumed the proportional hazards formulation
λ[k](y; tk−1 Z[k](y)) = λ (y; t
(.) is an unspecified continuous function. A sensibility analysis is recommended for an appropriate choice of fixed time point tk−1. Infer-ence can be carried out without making assumptions about associationamong individual's gap times.
Clement and Strawderman (2009) proposed a method for estimating the parameters indexing the conditional means and variances of the gaptime distributions conditional on all the available explanatory covariateshistory as well as on past gap times. Precisely, their work deal with E Y [k] Y [1], ., Y [k−1], Z(u) : 0 ≤ u ≤
Y [j] = µk(θ) Var Y [k] Y [1], ., Y [k−1], Z(u) : 0 ≤ u ≤
Y [j] = σ2Vk(θ)2 (2) where µk(.) and Vk(.) are known scalar functions of the unknown pa-rameter θ. The scalar parameter σ2 > 0 is also to be estimated. Theproposed methodology is an adaptation of generalized estimating equa-tions for longitudinal data and permits the use of both time-fixed andtime-varying covariates, as well as transformations of the gap times.
Censoring is dealt with by imposing a parametric assumption on thecensored gap times. Simulations report the relative robustness to devia-tions from this assumption although this supposed adequacy is identifiedas a potential issue. It shall be emphasized that the parametric assump-tions in (1) and (2) bear on the two first moments of the conditional distribution and not on the conditional hazard function itself. The Rpackage condGEE implements this conditional GEE for recurrent eventgap times.
Note also that all these methods do not account the possibility of The main issue with conditional models lies in the more or less ques- tionable assumptions made to incorporate past events.
Models for Multivariate Functions
Several non-parametric statistical analysis have been proposed for jointinference on consecutive gap times through multivariate functions inthe absence of death but accounting for non-informative censoring. Thenonparametric approach is quite classical to this aim. Nonparametricstatistics have the benefit of avoiding too restrictive assumptions espe-cially regarding would-be independence or memoryless-type conditions.
It is a good method to understand basics and to produce descriptiveresults. It also allows a first investigation of effects of covariates shownby stratifying data into groups.
Such an approach was originally developed in Visser (1996). The author considered joint nonparametric estimation for two successive du-ration times in the presence of independent right-censoring restricted tothe setting where the gap times and the censoring variable are discrete.
His method can deal with situations where censoring may depend uponprevious gap times but relies on estimating the cumulative conditionalhazard of the second gap time given the first one and therefore discretecensoring time and gap times are mandatory.
Wang and Wells (1998) studied the same problem but for any arbi- trary distributions of the gap times and the censoring variable. Theyalso considered joint nonparametric estimation for two successive du-ration times.
They proposed an estimator for the bivariate survival function of (Y [1], Y [2]) by estimating the cumulative conditional hazardof Y [2] given Y [1] > y1. The estimator was shown to be consistent andasymptotically normal, but do not guarantee a non-negative weightingof the data. Moreover, no analytical variance expression is given due tothe complicated expression of the estimator.
Lin et al. (1999) proposed a nonparametric estimator for the joint distribution function of the gap times. Their estimator is based on theinverse probability censoring weighted method used with the Kaplan-Meier estimator.
To enable comparison with other proposals, let us Modeling and Inferential Thoughts for Consecutive Gap Times introduce the observable gap times T [1] = min(Y [1], C) and T [2] =min(Y [2], (C − Y [1])I(Y [1] ≤ C)), the observable total duration timeT = min(Y [1] + Y [2], C) and the observable indicator variable δ =I(Y [1] + Y [2] ≤ C) when C is the censoring variable independent of (Y [1], Y [2]). Let (T i, δi) for i = 1, ., n be i.i.d. replications of (T [1], T [2], T, δ). Lin et al. (1999)'s estimator of P[Y [1) ≤ y1, Y [2] ≤ y2]can be written as 1 ∑ I(T ∑ I(T ≤ y i=1 1 − b Gn is a suitable Kaplan-Meier estimate of the censoring distri- bution function G. However, their estimator is not always a properdistribution function in that it may have negative mass points thoughit converges to a proper distribution function as n goes to ∞. Meira-Machado and Moreira (2010) found out in simulation studies that thisestimator is almost unbiased but may have important variance.
Van der Laan et al. (2002) considered more general problems of estimation that can be exploited for successive gap times. They alsoused inverse probability of censoring weights techniques. The statisticalnovelty of their approach lies in the derivation of locally efficient one-stepestimator.
In a more general situation of dependent censoring including the present setting as special case, Van Keilegom (2004) derived a nonpara-metric estimator for the bivariate and marginal distribution functionsof two gap times. The proposal by Van Keilegom (2004) consists ofwriting the joint distribution function of (Y [1], Y [2]) as an average of theconditional distribution F2 1(y2 y) = P[Y [2] ≤ y2 Y [1) = y] ie as P[Y [1] ≤ y1, Y [2] ≤ y2] = where F1(y) = P[Y [1] ≤ y1]. The conditional Kaplan-Meier estimator ofBeran (1981) is used to estimate F2 1. This relies on a kernel smoothingaround Y [1] = y with the modification that only uncensored observa-tions of T [1] are allowed in the window. The practical choice of theaforementioned window may be a limiting factor to a more frequent useof this estimator even though this choice is reported as non crucial andif a bootstrap procedure is advocated for.
Alvarez and Meira-Machado (2008) proposed another non- parametric estimator of bivariate distribution function of two consecu-tive gap times. The estimator of de U˜ Alvarez and Meira-Machado (2008) is a weighted bivariate distribution function of the form Wi I(T where the weight Wi is the Kaplan-Meier weight attached to Ti whenestimating the marginal distribution of Y [1] + Y [2] from the observablerandom variables (Ti, δi). Their estimator is a proper distribution func-tion contrarily to the proposal of both Wang and Wells (1998) and Linet al. (1999). Simulations revealed that this new estimator is reasonablyunbiased and may achieve efficiency levels clearly above the previous pro-posals, which is promising. However, theoretical investigation is neededto get general conclusions. It is also noted that the method can easilybe extended to cope with more than two successive gap times. Meira-Machado and Moreira (2010) found out in simulation studies that thisestimator is almost unbiased but may still have important variance.
In general, the prize to pay for the absence of restrictive assump- tions is a lack of efficiency. Some methods however exists to deal withthis issue. Presmoothing techniques may be useful to gain efficiency.
The idea of presmoothing goes back at least to Dikta (1998), see alsoDikta (2000), (2001) and Dikta et al. (2005). Presmoothing consists ofreplacing the censoring indicator by a smooth fit of a binary regressionof the indicator on observable gap times. This replacement usually re-sults in estimators with improved variance. That is why de U˜ and Amorim (2011) applied the idea of presmoothing to the estimationof the bivariate distribution function of censored gap times. As in thepaper by de U˜ Alvarez and Meira-Machado (2008), the estimator of Alvarez and Amorim (2011) is a weighted bivariate distribution function of the form Wi I(T [1] ≤ y1, T [2] ≤ y2) Wi now uses a presmoothed version of the preceding Kaplan-Meier weight Wi. The consequence of this is that the estimatorde U˜ Alvarez and Amorim (2011) can attach positive mass to pair of gap times with censored second gap times which is not the case with theestimator of de U˜ Alvarez and Meira-Machado (2008). Note that in the limiting case of no presmoothing, the estimator de U˜ Amorim (2011) reduces to that of de U˜ Alvarez and Meira-Machado (2008). A simulation study by Meira-Machado and Moreira (2010) logi-cally concluded that the presmoothed estimator improves efficiency with Modeling and Inferential Thoughts for Consecutive Gap Times respect to the estimator of de U˜ Alvarez and Meira-Machado (2008) but its bias is greater.
Van Keilegom et al. (2011) considered a non-parametric location- scale model for the two first gap times assuming that the vector of gaptimes (Y [1], Y [2]) satisfies Y [2] = m(Y [1]) + σ(Y [1]) ε where the functions m and σ are smooth and ε is independent of Y [1].
This allows the transfer of tail information from lightly censored areasto heavily ones. Under this model, the authors proposed estimators ofP[Y [1] ≤ y1, Y [2] ≤ y2], P[Y [2] ≤ y2 Y [1] = y1] and other related quan-tities. In a related paper, Meira-Machado et al. (2011) discussed thepractical implementation and performance of the aforementioned esti-mators and proposed some modifications. In an extensive simulationstudy, the good performance of the method is shown. The main limita-tion of their work lies in the fact that the adequacy of the model to thedata needs to be tested. However, the authors mentioned that derivingsuch a test in the present setting is far from being straightforward.
The R package survivalBIV is much helpful to calculate the different estimates for the bivariate distribution function.
We already mentioned the possible use of presmoothing to improve efficiency. General and testable assumptions such as Koziol-Green modelalso termed as informative censoring (Koziol and Green (1976), Chengand Lin (1987)) or proportionality constraints (Dauxois and Kirmani(2003), Geffray and Guilloux (2011)) can also be used leading to moreefficient semi-parametric inference under not so much restrictive assump-tions. Adekpedjou et al. (2010) adopted this strategy to tackle theefficiency problem.
Two-sample tests have been briefly considered in the literature. Lin and Ying (2001) proposed several classes of two-sample nonparametricstatistics for comparing the gap time distributions based on the nonpara-metric estimator of the gap time distribution given by Lin et al. (1999).
These statistics are analogous to familiar censored data statistics, suchas weighted log-rank statistics.
Huang (2000) proposed a semi-parametric accelerated failure time model to compare two treatment groups in terms of their successivegap times. Let ∆ be the indicator function that takes value 1 whenthe subject is in the first treatment group and 0 when the subject is inthe second treatment group. Specifically, the author parametrized thegroup effect on gap times and survival time by a scale transformation assuming that the random variables exp(β1∆)Y [1], ., exp(βk ∆)Y [k0] follow an unspecified multivariate continuous distribution function thatis independent of ∆. A log-rank type statistic is then derived.
Joint regression models have also been considered.
Huang (2002) considered multivariate accelerated failure time models for which the variables log Y [k] for k = 1, ., k0 follow a multivariatelocation-scale distribution of the form: log Y [k] = β′k.Z[k] + ε[k]
where (ε[1], ., ε[k0]) have an unspecified joint distribution that does not
depend on the vector of episode-specific covariates Z[k]. Note that infer-
ence is robust to misspecification of the gap time association structure
at the expense of strong assumptions on the censoring mechanism. It
turns out that, in this paper, censoring is assumed to be independent of
both covariates and the recurrent event process. Moreover, it is implic-
itly considered that each subject can experience at most k0 events. This
model may consequently appear less suited when the numbers of events
vary substantially across subjects.
He and Lawless (2003) presented multivariate parametric regression models for proportional hazards specified either within a copula model orwithin a frailty model. The method employs flexible piecewise constantor spline specifications as baseline hazard functions in either models.
Because all the models considered are parametric, ordinary maximumlikelihood can be applied. The adequacy to the parametric assumptionsis crucial to get unbiased estimates which may be a drawback.
All these methods, however, assume that recurrent events are not terminated by death during the study. Some efforts have been madeto account for death in a joint analysis for two-sample comparison pur-poses. Chang (2000) proposed a semi-parametric accelerated failure timemodel to compare two treatment groups jointly in terms of their suc-cessive gap times and survival time. This model is similar to that ofHuang (2000) but accounts for death. Let ∆ be the indicator func-tion that takes value 1 when the subject is in the first treatment groupand 0 when the subject is in the second treatment group.
cally, the author parametrized the group effect on gap times and sur-vival time by a scale transformation assuming that the random variablesexp(α∆)D, exp(β1∆)Y [1], ., exp(βk ∆)Y [k0] follow an unspecified mul- tivariate continuous distribution function that is independent of ∆. Alog-rank type statistic is then derived.
As a brief summary, models for multivariate functions mostly belong to the realm of non-parametric statistics in the absence of covariates Modeling and Inferential Thoughts for Consecutive Gap Times information. Fewer papers are available for multivariate functions inthe regression framework or in the presence of death.
Nonparametric Estimation of Cause-Specific
Distributions
In the recurrent events with death framework, functions describing thestochastic dynamics in the tree of Figure 3 can be much useful. Theapproach of Li and Lagakos (1997) and Derzko and Leconte (2004)who treated death as a competing risk acting at each recurrence canbe adopted for that purpose. They modeled the terminal event as a de-pendent competing event for each recurrent event i.e. they treated thefailure time for each recurrence as the first occurrence of the recurringevent or terminating event whichever came first. Thus, for each recur-rence, the patient is submitted to two dependent competing risks (REand death) in the presence of independent right-censoring provided heor she survived the previous occurrences. These step-by-step competingrisks models do not specify the association structure between recurrentevents and death. The work is centered on crude functions since theseare the only identifiable quantities without any assumptions regardingthe association structure among the competing risks. Non-parametricinference under minimal assumption is investigated.
At risk subjects Figure 3: Competing risks at each recurrence in the presence of independentcensoring (RE = recurrent event).
We assume that the observed data consist of i.i.d.
(Y [0], ., Y [K], (D ∧ C) − Y [k], I(D ≤ C)) where D is the death time, C is the independent right-censoring. The number K ∈ N is ran-dom as in Wang and Chang (1999) and Pe˜ na et al. (2001), Y [0] is set as 0, if K ≥ 1, the Y [k] for k = 1, ., K are the observed gap times untila recurrent event while the last gap time ends either with a death or acensoring event.
With these remarks in view, the functions that can serve as useful descriptive devices are the following. We consider for y1, y2 ≥ 0: F [1(2)](y1) = P D ≤ y1, Y [1] > D , F [1(1),2(1)](y1, y2) = P Y [1] ≤ y1, Y [1] ≤ D, Y [2] ≤ y2, Y [2] ≤ D − Y [1] , F [1(1),2(2)](y1, y2) = P Y [1] ≤ y1, Y [1] ≤ D, D − Y [1] ≤ y2, Y [2] > D − Y [1] . This can be straightforwardly extended to further recurrences provided the data are not too sparse.
Let FD be the distribution function of D. Denote by C the non- negative random variable that stands for the independent right-censoringwith distribution function G. Let H be the distribution function definedby 1 − H = (1 − FD)(1 − G) and let τH = sup{x : H(x) < 1} be theright-endpoint of the distribution function H. The functions (3) to (5)can be consistently estimated and it can be shown that the correspond-ing estimators have an asymptotic Gaussian behavior on compact setssuch that the corresponding total observation time is inferior to τH . No-tice that F [1](y1) is estimable only if y1 < τH , that F [1(1),2](y1, y2) isestimable only if y1 + y2 < τH and so on. The objective of this section isto justify nonparametric estimation for the functions displayed in Equa-tions (3) to (5) without any assumption regarding either the dependencestructure among the multiple endpoints.
For ease of exposition, note that the observable random variables can be coded as follows.
• Let K +1 (with K ∈ N) be the total number of observed events for a given individual (including recurrent events, death and censoringevents).
For k = 1, . . , K + 1, let T be the random variable that stands for the gap time between the (k − 1)-th and the k-th event and set • For k = 1, . . , K + 1, the random variable Modeling and Inferential Thoughts for Consecutive Gap Times 0 if the k-th event is censored J [k] = 1 if the k-th event is a recurrent event 2 if the k-th event is a death indicates the nature of the k-th observed event.
We suppose that observations are taken on an i.i.d. sample of n individuals. For i = 1, . . , n, the data for the i-th individual consists ofKi + 1 couples where Ki is the number of observed (non-fatal) recurrent events. For k = 1, . . , Ki + 1, the k-th couple is given by (T which is distributed as (T , J [k]).
Estimation of the Censoring Distribution Function
An estimate of the censoring distribution function G will be used. Thissubsection deals with this preliminary step.
As noted in Section 1, the last observation for a given patient is either a censoring time or a death time. For a given patient, we do not observeboth the censoring event and the death but only the first event thatoccurs. Since the death and the censoring processes are independent,the censoring distribution function may be estimated by the Kaplan-Meier estimator based on the total observation time for each patient i.e.
on the data (T ) for i = 1, . . , n.
The Kaplan-Meier estimator of the censoring distribution function G is given for t ≥ 0 by: i ≤ t, Ji n(t) = 1 − I (Tℓ ≥ Ti) If there are ties between recurrent event times and censoring times, theKaplan-Meier estimator of G cannot be obtained by using the indicatorstatus equal to zero. In such cases, the R package prodlim provides auseful alternative to estimate the censoring distribution.
"Plug-in" Estimation of the Functions of Interest
To derive an estimator for the functions F [1(2)], F [1(1),2(1)] and F [1(1),2(2)],we introduce the following distribution functions for y1, y2 ≥ 0: H[1(1,2)](y1) = P T 1, J [1] = 2 H[1(1,1),2(1,j)](y1, y2) = P T , j = 1, 2. 2, J [1] = 1, J [2] = j For y ≥ 0, we obtain the following relation: H[1(1,2)](y) = P Y [1] ≤ y, Y [1] ≤ C, C[1] = 2 I (u ≤ y, u ≤ c) G(dc)F [1(2)](du) 1 − G−(u) F [1(2)](du) . with G− being the left-continuous modification of G.
F [1(2)](y) can be written in terms of the estimable functions G andH[1(1,2)]: F [1(2)](y) = u≤y 1 − G−(u) We can obtain in the same way for j = 1, 2 and y1, y2 ≥ 0 that I (u ≤ y1, v ≤ y2, c ≥ u + v) × G(dc)F [1(1),2(j)](du, dv) 1 − G−(u + v) F [1(1),2(j)](du, dv) F [1(1),2(j)](y1, y2) = 1 − G−(u + v) Consequently, we propose "plug-in" estimates of the functions F [1(2)]and F [1(1),2(j)] for j = 1, 2 by means of "plug-in" estimators denoted respectively by b for j = 1, 2. These estimators are obtained by replacing G by its Kaplan-Meier estimator defined in Sub-section 3.1 and H[1(1,2)] and H[1(1,1),2(1,j)] by their empirical counterpartswhich are defined respectively for y1, y2 ≥ 0 by: y1, J (y1, y2) = y1, T y2, J Modeling and Inferential Thoughts for Consecutive Gap Times Consequently, we let for y1, y2 ≥ 0 n (u) (y1, y2) = j = 1, 2 n (u + v) n is the left-continuous modification of b 1. For any σ < τH , the estimator b is strongly consistent on [0, σ] for F [1(2)]. 2. For any σ < τH , the estimators b are strongly consistent for F [1(1),2(j)], for j = 1, 2, on the set Tσ = {(y1, y2) : y1 + y2 < σ}. Remark 3.1.
If y1 is taken equal to ∞ in the definition of the esti- , one would have b (∞) equal to b is the last order statistic of the sample (T ) → F [1(2)](τ 1(2) ∧ τG) holds in probability where τ1(2) is the right-endpoint of F [1(2)] and where τG is the right-endpointof G. But F [1(2)](τ1(2) ∧ τG) may be strictly inferior to F [1(2)](∞). Thisis fulfilled in particular if τG < τ1(2) which is the case in a clinical trialfor example where this is to hope that some patients won't experience arecurrence by the end of study. This situation would lead to a biased es-timation of F [1(2)](∞). The same kind of restriction holds for the otherestimators mentioned here.
Assume that G is continuous. For any σ < τH , the − F [1(2)] and n F [1(1),2(j)] − F [1(1),2(j)] for j = 1, 2 converge jointly in distribution to zero-mean Gaussian pro-cesses in the Skorohod space of c ag functions on Tσ. Remark 3.2.
The condition that G is continuous is restrictive since it does not allow for a fixed time to follow-up. Further work would beneeded for such an extension.
To save place, the large sample arguments are purposefully sketchy.
Proof of Proposition 3.1. We decompose Fn
F [1(1),2(j)] for j = 1, 2 into (y1, y2) − b F [1(1),2(j)](y1, y2) = 1 − G−(u + v) n (u + v) (du, dv) − H[1(1,1),2(1,j)](du, dv)) 1 − G−(u + v) We carry out integration by parts on the second term in the aboveequality and get straightforwardly (y1, y2) − b t≤σ  1 − G(t) 1 − G(σ) (y1, y2) The fact that G(σ) < 1 together with Glivenko-Cantelli's theorem validwith and without independent right-censoring give the required almost sure convergence on Tσ. The proof is identical for b Proof of Proposition 3.1. First, we endow the space of cadlag func-
tions on Tσ with the appropriate topology. This can be obtained by
transporting the Skorohod topology of the space of cadlag functions on
[0, 1]2 build in Neuhaus (1971) since the spaces Tσ and [0, 1]2 are home-
omorphic.
The weak convergence result is then obtained by empirical processes techniques. It relies on appropriate decomposition of the processes − F [1(2)] and n Fn − F [1(1),2(j)] for j = 1, 2 that permits to apply the joint convergence of univariate and multivariateempirical processes based on the observed data. The functional delta-method as in Andersen et al.
(1993) is also of much use.
ingredient is the use of the existing results for the Kaplan-Meier processwhich makes the assumption that G is continuous necessary. Let us begin with the decomposition of − F [1(1),2(j)] for j = Modeling and Inferential Thoughts for Consecutive Gap Times 1, 2. The decomposition of − F [1(2)] is left to the reader.
(y1, y2) − F [1(1),2(j)](y1, y2) y2 Hn (du, dv) − H[1(1,1),2(1,j)](du, duv) 1 − G(u + v) n (u + v) − G(u + v) H[1(1,1),2(1,j)](du, dv) (1 − G(u + v))2 Gn(u + v) − G(u + v) (1 − G(u + v))2 n (u + v) − G(u + v) 1 − G(u + v) n (u + v) = I1(y1, y2) + I2(y1, y2) + I3(y1, y2) + I4(y1, y2) . Terms I3 and I4 are negligible uniformly on Tσ thanks to the functionaldelta-method and to the weak convergence of both the Kaplan-Meier process and the empirical process I2 needs to be further decomposed. Let H(0) be the censoring subdis-tribution function defined by H(0)(y) = P[T ≤ y, JK+1 = 0], let its empirical counterpart be defined by n (y) = I(Ti ≤ y, JKi+1 = 0) and let Hn be the empirical counterpart of H defined by Hn(y) = I(Ti ≤ y). Applying the methods and results of Cs¨ o (1996), we can state since G is assumed continuous that n(t) − G(t) d(Hn (s) − H(0))(s) 1 − G(t) 1 − H−(s) n (s) − H −(s) dH(0)(s)  = oP √ (1 − H−(s))2 Consequently, the process − F [1(1),2(j)] is asymptotically y2 Hn (du, dv) − H[1(1,1),2(1,j)](du, dv) 1 − G(u + v) (u+v)− d(Hn (s) − H(0))(s) 1 − H−(s) n (s) − H −(s) dH(0)(s)H[1(1,1),2(1,j)](du, dv) . (1 − H−(s))2 It remains to get the joint convergence of n(Hn − H(0)), n(Hn − −H[1(1,1),2(1,j)]) in order to apply once again the functional delta-method and get the result. To see this, set T R = and J R = J Ki+1 for i = 1, ., n such that K i + 1 ≥ 3. Set also T R = K+1 T with the sum being null if the summation index set is void and J R = J K+1. Then, decompose (Hn − H(0)) into n (y) − H (0)(y) ≤ y, J = 0) − P[T ≤ y, J[1] = 0] ≤ y, J = 1, J = 0) ≤ y, J[1] = 1, J[2] = 0] ≤ y, J = J = 1, J ≤ y, J[1] = J[2] = 1, J[R] = 0 and decompose (Hn − H) into Hn(y) − H(y) = ≤ y, J ̸= 1) − P[T ≤ y, J[1] ̸= 1] Modeling and Inferential Thoughts for Consecutive Gap Times ≤ y, J = 1, J ̸= 1) ≤ y, J[1] = 1, J[2] ̸= 1] ≤ y, J = J = 1, J ≤ y, J[1] = J[2] = 1, J[R] ̸= 1 Pollard's (1982) theorem valid for the empirical processes indexed by the VC-class of sets { 2 [0, y i] : y1 + y2 ≤ σ} concludes the argument.
Integration by parts permits to have the empirical processes in the in-tegrand rather than in the measure of integration. This technicality isleft to the reader. 2 The variance function of the limiting process is not mentioned. Ob- taining it through empirical processes techniques is quite cumbersome.
Moreover, a classical problem with competing risks is the variance ex-plosion which entails far too wide confidence bands, see e.g. Geffray(2009), making this calculus less interesting.
Martingale methods have been successfully used for survival analysispurposes from the mid 1970's. These developments go back to Aalen'swork, see e.g. Aalen (1978a), (1978b), Aalen et al. (1978) then moved totwo-sample tests, Cox's proportional hazards regression model, Markovtransition probabilities estimation among many other situations, see e.g.
Gill (1980), (1983), (1994), Andersen et al. (1993). It emerges that thecounting process and stochastic integral approach provides relativelysimple methods of inference in some situations where standard methodsof inference are too cumbersome or require too restrictive assumptions.
The martingale inference methods provide systematic methods for mo-ment calculation, establish asymptotic normality of empirical processesand gave rise to a variety of results in settings where a single time runs,in particular, it enables to extend asymptotic results up to the last or-der statistic of the observation instead of restricting their validity tocompact time-intervals.
The martingale methods of use in survival analysis could be gen- eralized to our setting where multiple times run. This would requiremulti-parameter counting processes and martingales.
time-continuous multi-parameter martingale and stochastic integral hasdeveloped extensively from the 1970's, see e.g.
Bickel and Wichura (1971), Zakai (1981), Zakai et al.
(1974), (1976), Merzbach (1988), Ivanoff (1996), with a special emphasis to set-indexed martingales inthe nineties, see e.g.
Ivanoff and Merzbach (2000) for an extensive study of this subject. These probability results opened a new perspec-tive for statisticians involved with multi-time periods inference prob-lems. The use of martingale methods for multi-parameter problems insurvival analysis was initiated by Pons (1986) in the setting of bivari-ate survival function estimation. In a tremendous paper, Ivanoff andMerzbach (2002) developed a general model for survival analysis wherecensored data are parametrized by sets instead of time points. Dis-appointingly, their work is of little help for our purpose. Our specificcensoring scheme complicate technical matters considerably and, in par-ticular, invalidate the direct use of their methods. However we anticipatethat some progress could be made in this direction.
We pointed out earlier the fact that the pure nonparametric approach may suffer from a lack of efficiency. Considering a presmoothed versionof the estimate of functions (3) to (5) as in Cao et al. (2004), (2005),Amorim et al. (2011), de U˜ na-Alvarez and Amorim (2011) could be an interesting remedy.
In the framework of Section 3, conditional analysis could be interest- ing since it allows dynamical prediction while incorporating a patient'shistory. The estimation of conditional probabilities such as P Y [2] ≤ y2, Y [2] ≤ D − Y [1] Y [1] = y1, Y [1] ≤ D P D ≤ y2, Y [2] > D − Y [1] Y [1] = y1, Y [1] ≤ D is currently under investigation via projection methods without any as-sumptions regarding dependence structure of successive gap times.
It is worth mentioning that investigators increasingly encounter data- sets in which some patients are expected to be cured. This is a seri-ous matter because a patient surviving the trial is considered censoredwhereas the patient is cured if he or she will never experiment the eventunder study. The difficulty comes from the fact that a cure can neverbe observed due to a finite monitoring time. To address this problem,cure rate models have been proposed and have received intensive atten-tion for their ability to account for the probability of a patient being Modeling and Inferential Thoughts for Consecutive Gap Times cured, see e.g. Yin and Ibrahim (2005). Recently, some progress havebeen made to incorporate the possibility of cure into recurrent eventsmodeling. Rondeau (2010) developed a cure frailty model to evaluatetime-dependent medical treatment effects on the times to recurrenceamong the uncured patients and on the cure probability. The probabil-ity of cure here may evolve with time and is defined as the probabilityof not developing further event after each event. Rondeau et al. (2011)compared several forms of cure rate model within a frailty model for therecurrent event part. To analyze recurrent events, it is first necessaryto define the cured proportion to be modeled. The first model considersthat immune patients are those who are not expected to experience theevents of interest over a sufficiently long time period. The other investi-gated models account for the possibility of cure after each event i.e. theprobability of cure may evolve with time. The focus is placed on timesto recurrence and death is accounted for.
Account for the possibility of cure should be dealt with in the frame- work of joint multivariate approach as well as th possibility of incorpo-rating covariates acting on the uncured population survival. Investiga-tion of both their practical and theoretical properties should be carefulreported for applied purposes.
Another point that is worth mentioning is the issue due to non- reliable cause of death. In this situation, relative survival models canbe of use, see e.g. Lambert et al. (2010). Subjects may die of thedisease they are diagnosed with but they may also die of something else.
Deaths due to another cause than strictly the disease under study canbe broadly classified into "totally independent death" i.e. death froma cause related neither to the disease under study nor to the treatmentand "possibly related death". The "totally independent death" usuallyconstitutes part of the independent right-censoring process and is notan issue.
But the class "possibly dependent" leads to difficulties in interpreting the results. Interest mostly lies in mortality strictly due tothe disease of interest and not to related causes. How to classify, forexample, deaths due to treatment complications? Consider a patientdiagnosed with lung cancer who dies following a myocardial infarction.
Do we classify this death as ‘due entirely to lung cancer' or ‘due entirelyto other causes' ? There may also exist problems with cause-specificdeath distribution due to inaccuracy of death certificates. An alternativeto cause-specific distribution estimation is then to model relative survivalor its converse which is termed as excess mortality. Suppose that, S∗(t)is the expected survival. Then the total survival S(t) can be written as the product of the relative survival R(t) and the expected survival S∗(t)i.e. S(t) = S∗(t)R(t). Relative survival is often preferred over cause-specific survival for the study of cancer patient survival. This issue isparticularly relevant here where possible applications are infarction orcancer recurrence and related death and is worth investigating.
A last interesting point to note is that in our setting some events may be rare. For instance, Cui et al. (2010) noted that the chance ofhaving two myocardial infarction events within 5 years was low amongall participants in the LIPID study. As a consequence, when analyzingthe pre-specified set, say (Y [1], Y [2]), the second gap times Y [2] won't beavailable for many patients leading to efficacy problems. The work ofBuyske et al. (2000) on two-sample log-rank statistics when the survivalevent is rare could be extended to the present setting.
Aalen, O. (1978a), Nonparametric estimation of partial transition prob- abilities in multiple decrement model. The Annals of Statistics, 6,
534–545.
Aalen, O. (1978b), Nonparametric inference for a family of counting processes. The Annals of Statistics, 6, 701–726.
Aalen, O. and Johansen, S. (1978), An empirical transition matrix for nonhomogeneous Markov chains based on censored observations.
Scandinavian Journal of Statistics, 5, 141–150.
Adekpedjou, A., Pe˜ na, E. A., and Quiton, J. (2010), Estimation and efficiency with recurrent event data under informative censoring.
Journal of Statistical Planning and Inference, 140, 597–615.
Amorim, A. P., de U˜ Alvarez, J., and Meira-Machado, L. (2011), Presmoothing the transition probabilities in the illness-death model.
Statistics and Probability Letters, Special Issue: Statistics in Bio-
logical and Medical Sciences, 81, 797–806.
Andersen, P., Borgan, O., Gill, R. D., and Keiding, N. (1993), Statis- tical Models Based on Counting Processes. New-York: Springer-Verlag.
Andersen, P. K. and Pohar Perme, M. (2008), Inference for outcome probabilities in multi-state models. Lifetime and data analysis,
14, 405–431.
Modeling and Inferential Thoughts for Consecutive Gap Times Beran, R. (1981), nonparametric regression with randomly censored survival data. Technical report, Univ. California, Berkeley.
Bickel, P. J. and Wichura, M. J. (1971), Convergence Criteria for Mul- tiparameter Stochastic Processes and Some Applications. Annals
of Mathematical Statistics, 42(5), 1656–1670.
Bowker, T., Clayton, T., and co authors (1996), A British Cardiac Society survey of the potential for secondary prevention of coro-
nary disease: ASPIRE (Action on Secondary Prevention through
Intervention to Reduce Events). Heart, 75, 334-342.
Buyske, S., Fagerstrom, R. and Ying, Z. (2000), A class of weighted log-rank tests for survival data when the event is rare. Journal of
the American Statistical Association, 95(449), 249–258.
acome, M. A. (2004), Presmoothed kernel density esti- mator for censored data. Journal of Nonparametric Statistics, 16,
289–309.
opez de Ullibarri, I., Janssen, P., and Veraverbeke, N. (2005), Presmoothed Kaplan-Meier and Nelson-Aalen estimators. Journal
of Nonparametric Statistics, 17, 31–56.
Chang, S.-H. (2000), A two-sample comparison for multiple ordered event data. Biometrics, 56, 183–189.
Chang, S.-H. (2004), Estimating marginal effects in accelerated fail- ure time models for serial sojourn times among repeated events.
Lifetime Data Analysis, 10, 175–190.
Chang, S. H. and Wang, M. C. (1999), Conditional regression analy- sis for recurrence time data. Journal of the American Statistical
Association, 94, 1221–1230.
Chen, Y. Q., Wang, M. C., Huang, Y. (2004), Semiparametric regres- sion analysis on longitudinal pattern of recurrent gap times. Bio-
statistics, 5(2), 277–290.
Cheng, P. E. and Lin, G. D. (1987), Maximum likelihood estimation of a survival function under the Koziol-Green proportional hazards
model. Statistics and Probability Letters, 5, 75–80.
Clement, D. Y. and Strawderman, R. L. (2009), Conditional GEE for recurrent event gap times. Biostatistics, 10(3), 451–467.
Cook, R. J. and Lawless, J. F. (2008), The Statistical Analysis of Re- current Events. New-York: Springer-Verlag.
o, S. (1996), Universal Gaussian approximations under random censorship. The Annals of Statistics, 24, 2744–2778.
Cui, J., Forbes, A., Kirby, A., Marschner, I., Simes, J., Hunt, D., West, M. and Tonkin, A. (2010), Semi-parametric risk prediction
models for recurrent cardiovascular events in the LIPID study.
BMC Medical Research Methodology, 10(27), 1471–2288.
Dauxois, J.-Y. and Kirmani, S. N. U. A. (2004), On testing the pro- portionality of two cumulative incidence functions in a competing
risks setup. Journal of Nonparametric Statistics, 16, 479–492.
Alvarez, J. and Amorim, A. P. (2011), A semiparametric esti- mator of the bivariate distribution function for censored gap times.
Biometrical Journal, 53, 113-127.
Alvarez, J. and Meira-Machado, L. (2008), A simple estima- tor of the bivariate distribution function for censored gap times.
Statistics and Probability Letters, 78, 2240-2445.
Dikta, G. (1998), On semiparametric random censorship models. Jour- nal of Statistical Planning and Inference, 66, 253–279.
Dikta, G. (2000), The strong law under semiparametric random cen- sorship models. Journal of Statistical Planning and Inference, 83,
1–10.
Dikta, G. (2001), Weak representation of the cumulative hazard func- tion under semiparametric random censorship models. Statistics,
35, 395–409.
Dikta, G., Ghorai, J., and Schmidt, C. (2005), The central limit theo- rem under semiparametric random censorship models. Journal of
Statistical Planning and Inference, 127, 23–51.
Du, P. (2009), Nonparametric modeling of the gap time in recurrent event data. Lifetime and Data Analysis, 15, 256–277.
Du, P., Jiang, Y. and Wang, Y. (2011), Smoothing spline ANOVA frailty models for recurrent event data. Biometrics, 67, 1330–1339.
Modeling and Inferential Thoughts for Consecutive Gap Times Duchateau, L. and Janssenn, P. (2008) , The Frailty Model. New York: alez, I., Busse, J. W., Heels-Ansdell, D., and et al. (2007), Problems with use of composite endpoints in cardiovascular trials:
systematic review of randomised controlled trials. BMJ Journal,
334(7597), 786–793.
Geffray, S. (2009), Strong approximations for dependent competing risks with independent censoring. TEST, 18(1), 76-95.
Geffray, S. and Guilloux, A. (2011), Maximum likelihood estimator for cumulative incidence functions under proportionality constraint.
Sankhya, series A, 73(2), 303–328.
Gent, S. and co-authors, (1996), A randomised, blinded trial of clopi- dogrel versus aspirin in patients at risk of ischaemic events (CA-
PRIE). The Lancet, 348, 1329–1339.
Gill, R. D. (1980), Censoring and Stochastic Integrals. Mathematical Centre Tracts 124, Mathematisch Centrum, Amsterdam.
Gill, R. D. (1983), Large sample behavior of the product-limit estimator in the whole line. The Annals of Statistics, 11, 49–58.
Gill, R. D. (1994), Glivenko-Cantelli for Kaplan-Meier. Mathematical Methods of Statistics, 3, 76–87.
Gill, R. D. and Keiding, N. (2010), Product-limit estimators of the gap time distribution of a renewal process under different sampling
patterns. Lifetime Data Analysis, 16(4), 571–579.
He, W. and Lawless, J. F. (2003), Flexible Maximum Likelihood Meth- ods for Bivariate Proportional Hazards Models, Biometrics, 59(4),
837–848.
Hougaard, P. (2000), Analysis of Multivariate Survival Data.
Huang, Y. (2000), Two-Sample Multistate Accelerated Sojourn Times Model. Journal of the American Statistical Association, 95(450),
619–627.
Huang, Y. (2002), Censored regression with the multistate accelerated sojourn times model. Journal of the Royal Statistical Society, se-
ries B, 64(1), 17–29.
Huang, Y. and Chen, Y. (2003), Marginal regression of gaps between recurrent events. Lifetime Data Analysis, 9, 293–303.
Huang, X. and Liu, L. (2007), A joint frailty model for survival and gap times between recurrent events. Biometrics, 63, 389–397.
Huang, C.-Y., Luo, X., and Follmann, D. A. (2010), A model checking method for the proportional hazards model with recurrent gap
time data. Biostatistics, 12(3), 535–547.
Huang, Y. and Wang, C.-Y. (2000), Cox regression with accurate co- variates unascertainable: a nonparametric-correction approach.
Journal of the American Statistical Association, 95, 1209–1219.
Ivanoff, G. (1996), Stopping times for multi-parameter martingales.
Statistics and Probability Letters, 28, 111–114.
Ivanoff, G. and Merzbach, E. (2002), Random censoring in set-indexed survival analysis. Annals of Applied Probability, 12(3), 944–971.
Ivanoff, G. and Merzbach, E. (2000) , Set-indexed martingales. Boca Raten, FL: CRC Press, Jokhodar, M., Jacobsen, S. J., Reeder, G. S., Weston, S. A., and Roger, V. L. (2004), Sudden death and recurrent ischemic events after my-
ocardial infarction in the community. American Journal of Epi-
demiology, 159(1), 1040–1046.
Kessing, L. V., Olsen, E. W., and Andersen, P. K. (1998), Recurrence in affective disorder: analyses with frailty models. American Journal
of Epidemiology, 149, 404–411.
Kleist, P. (2007), Composite endpoints for clinical trials: current per- spectives. International Journal of Pharmaceutical Medicine, 21,
187–198.
Koziol, J. A. and Green, S. B. (1976), A Cram´ er-Von Mises statistic for randomly censored data, Biometrika, 63, 465–474.
Kvist, K., Gerster, M., Andersen, P. K., and Kessing, L. V. (2007), Non-parametric estimation and model checking procedures for marginal
gap time distributions for recurrent events. Statistics in Medicine,
26, 5394–5410.
Modeling and Inferential Thoughts for Consecutive Gap Times Lambert, P. C., Dickman, P. W., Nelson, C. P., and Royston, P. (2010), Estimating the crude probability of death due to cancer and other
causes using relative survival models. Statistics in Medicine, 29,
885–895.
Lawless, J. F., Wigg, M. B., Tuli, S., Drake, J., and Lamberti-Pasculli, M. (2001), Analysis of Repeated Failures of Durations, with Ap-
plication to Shunt Failures for Patients with Paediatric Hydro-
cephalus. Journal of the Royal Statistical Society, series C, (Ap-
plied Statistics), 50(4), 449–465.
Li, Q. and Lagakos, S. (1997), Use of the Wei-Lin-Weissfeld method for the analysis of a recurring and a terminating event. Statistics
in Medicine, 16, 925–940.
Lin, D. Y., Sun, W., and Ying, Z. (1999), Nonparametric estimation of the gap time distribution for serial events with censored data.
Biometrika, 86, 59–70.
Lin, D. Y., and Ying, Z. (2001), Nonparametric Tests for the Gap Time Distributions of Serial Events Based on Censored Data, Biomet-
rics, 57, 369–375.
e, C. and Chevret, S. (1999), Estimation of the treatment effect in a clinical trial when recurrent events define the endpoint. Statistics
in Medicine, 18(14), 1821–1829.
Marschner, I.,Colquhoun, D., Simes, R., Glasziou, P., Harris, P., Singh, B., Fridlander, D., White, H., Thompson, P., and Tonkin, A.
(2001), Long-term risk stratification for survivors of acute coronary
syndrome. Results from the long-term intervention with Pravas-
tatin in ischemic disease (LIPID) study. Amer. J. Cardiol., 38(1),
56-63.
Meira-Machado, L., de U˜ Alvarez, J., Cadarso-Su´ arez, C., Andersen, P. K. (2009), Multi-state models for the analysis of time to event
data. Statistical Methods in Medical Research, 18, 195–222.
Meira-Machado, L. and Moreira, A. (2006), Estimation of hte bivariate distribution function for censored gap times. LIDA, 12, 1367-1374.
Meira-Machado, L., Roca-Pardi˜ nas, J., Van Keilegom, I., Cadarso- arez, C. (2011), Estimation of transition probabilities in a non- Markov model with successive survival times. preprint.
Merzbach, E. (1988), Point Processes in the Plane. Acta Applicandae Mathematicae, 12, 79–101.
Neuhaus, G. (1971), On weak convergence of stochastic processes with multidimensional time parameter. Annals of Mathematical Statis-
tics, 42(4), 1285–1295.
na, E., Strawderman, R., Hollander, M. (2001), Nonparametric es- timation with recurrent event data. Journal of the American Sta-
tistical Association, 96, 1299–1315.
Pfeffer, M. A., Sacks, F. M., Moy´ e, L. A., Brown, L., Rouleau, J.
L., Hartley, L. H., Rouleau, J., Grimm, R., Sestier, F., Wicke-
meyer, W., and et al. (1996), Cholesterol and Recurrent Events:
a secondary prevention trial for normolipidemic patients. CARE
Investigators. American Journal of Cardiology, 76(9), 98–106.
Pocock, S. J. (1997), Clinical trials with multiple outcomes: A sta- tistical perspective on their design, analysis, and interpretation,
Controlled Clinical Trials, 18, 530–545.
Pollard, O. (1982), A central limit theorem for empirical processes.
Journal of Australian Mathematical Society, series A, 33, 235–
248.
Pons, O. (1986), A test of independence between two censored survival times. Scandinavian Journal of Statistics, 13, 173–185.
Rondeau, V., Mathoulin-Pelissier, S., Jacqmin-Gadda, H., Brouste, C., and Soubeyran, P. (2007), Joint frailty models for recurring events
and death using maximum likelihood estimation: application on
cancer events. Biostatistics, 8(4), 708–721.
Rondeau, V. (2010), Statistical models for recurrent events and death: application on cancer events. Mathematical and Computer Mod-
elling, 52, 949–955.
Rondeau, V., Schaffner, E., Corbi ere, F., Gonzalez, J. R., Mathoulin- elissier, S. (2011) Cure frailty models for survival data: applica- tion to recurrences for breast cancer and to hospital readmissionsfor colorectal cancer. Statistical Methods in Medical Research,DOI: 10.1177/0962280210395521.
Modeling and Inferential Thoughts for Consecutive Gap Times Schaubel, D. E. and Cai, J. (2004), Non-parametric estimation of gap time survival functions for ordered multivariate failure time data.
Statistics in Medicine, 23, 1885–1900.
Schaubel, D. E. and Cai, J. (2004), Regression methods for gap time hazard functions of sequentially ordered multivariate failure time
data. Biometrika, 91(2) 291–303.
Scun, L., Park, D.-H., and Sun, J. (2006), The additive hazards model for recurrent gap times. Statistica Sinica, 16, 919–932.
Strawderman, R. L. (2005), The accelerated gap times model. Biometrika, The LIPID Study Group (1995), Design features and baseline charac- teristic of the LIPID (Long-Term Intervention with Pravastatin in
Ischemic Disease) study: a randomized trial in patients with pre-
vious acute myocardial and/or unstable angina pectoris. Amer. J.
Cardiol, 76, 474-479.
The LIPID Study Group (1998), Prevention of cardiovascular events and death with Pravastatin in patients with coronary heart disease
and a broad range on initial cholesterol level. N. Eng. J. Med.,
339(19),1349-1357.
Therneau, T. M. and Grambsch, P. M. (2000), Modeling Survival Data: Extending the Cox Model. New York: Springer-Verlag.
Van der Laan, M. J., Hubbard, A. E., and Robins, J. M. (2002), Locally efficient estimation of a multivariate survival function in longitu-
dinal studies. Journal of the American Statistical Association, 97,
494–507.
Van Keilegom, I. (2004), A note on the nonparametric estimation of the bivariate distribution under dependent censoring. Journal of
Nonparametric Statistics, 16, 659–670.
Van Keilegom, I., de U˜ Alvarez, J., and Meira-Machado, L. (2011), Nonparametric location-scale models for successive survival times
under dependent censoring. Journal of Statistical Planning and
Inference, 141, 1118–1131.
Visser, M. (1996), Nonparametric estimation of the bivariate survival function with application to vertically transmitted AIDS. Biometrika,
83, 507–518.
Wang, M. C. (1999), Gap time bias in incident and prevalent cohorts.
Statistica Sinica, 9, 999–1010.
Wang, M. C. and Chang, S.-H. (1999), Nonparametric estimation of a recurrent survival function. Journal of the American Statistical
Association, 94, 146–153.
Wang, W. and Wells, M. T. (1998), Nonparametric estimation of suc- cessive duration times under dependent censoring. Biometrika, 85,
561–572.
Wong, E. and Zakai, M. (1974), Martingales and Stochastic Integrals for Processes with a Multi-Dimensional Parameter. Z. Wahrschein-
lichkeitstheorie verw. Gebiete, 29, 109–122.
Wong, E. and Zakai, M. (1976), Weak Martingales and Stochastic In- tegrals in the Plane. The Annals of Probability, 4(4), 570–586.
Yin, G. S. and Ibrahim, J. G. (2005), Cure rate models: a unified approach. The Canadian Journal of Statistics, 33, 559–570.
Zakai, M. (1981), Some Classes of Two-Parameter Martingales. The Annals of Probability, 9(2), 255–265.

Source: http://jirss.irstat.ir/article-1-210-en.pdf