Nonparametric estimation of the sojourn time distributions for a multipath model

全文

(1)J. R. Statist. Soc. B (2003) 65, Part 4, pp. 921–935. Nonparametric estimation of the sojourn time distributions for a multipath model Weijing Wang National Chiao-Tung University, Hsin-Chu, Taiwan [Received May 2000. Final revision April 2003] Summary. We use a multipath (multistate) model to describe data with multiple end points. Statistical inference based on the intermediate end point is challenging because of the problems of nonidentifiability and dependent censoring. We study nonparametric estimation for the path probability and the sojourn time distributions between the states. The methodology proposed can be applied to analyse cure models which account for the competing risk of death. Asymptotic properties of the estimators proposed are derived. Simulation shows that the methods proposed have good finite sample performance. The methodology is applied to two data sets. Keywords: Competing risks; Cure models; Dependent censoring; Identifiability; Illness–death models; Long-term survivors; Multiple end points; Semicompeting risks; Susceptibility. 1.. Introduction. 1.1. Background Many biomedical applications involve the analysis of multiple end points. Consider a study of bone marrow transplants for leukaemia patients. Some patients will experience recurrence of the malignancy before death but others may die without relapse. A second example is related to the study of a life-threatening acute condition (Betensky and Schoenfeld, 2001). Only a proportion of patients will leave the hospital alive and later die from other causes, whereas the remaining patients, who are not cured, will die in the hospital. These phenomena can be described by the multipath or multistate model that is depicted in Fig. 1. The terminal event, death, leads to an absorbing state whereas the non-terminal event, indicating some progression status, may be bypassed. Let .X, Y , C/ be the times to an intermediate end point, a terminal end point and an external censoring end point respectively. Censoring may be due to termination of the study period or the loss of patients to follow-up. Let {.Xi , Yi , Ci /, i = 1, . . . , n} be identically and independently distributed replications of .X, Y , C/. We shall assume that C is independent of both X and Y , whereas X and Y may be correlated. Statistical methods for analysing a terminal end point and a non-terminal end point are very different. Let S1 .t/ = Pr.X > t/. When death is of primary interest, the Kaplan–Meier estimator is the standard method for estimating the survival function of Y . Progression information can also be incorporated to improve estimation (Gray, 1994). When the prediction of individual survival times is of interest, we need to specify how the intermediate end point affects subsequent survival. Popular choices of model include Markov and semi-Markov models and the proportional hazards model, in which the progression information is treated as a time-dependent covariate. Andersen Address for correspondence: Weijing Wang, Institute of Statistics, National Chiao-Tung University, Hsin-Chu, Taiwan, Republic of China. E-mail: wjwang@stat.nctu.edu.tw  2003 Royal Statistical Society. 1369–7412/03/65921.

(2) 922. W. Wang State 1: healthy. XX≤Y. YY<X. State 2: progression. Y-X  X ≤ Y State 3: death. Fig. 1.. Multipath model. et al. (1993) and Hougaard (2000) contain many examples and good summaries of related results. When interest focuses on the progression of disease, statistical issues become more complicated and sometimes controversial (Prentice et al., 1978). Death precludes the biological development of progression whereas censoring terminates observation of the other two events. The major challenge comes from the fact that no progression information is available after death. For those who do not experience progression during their lifetime, X is not well defined and therefore its distribution is non-identifiable without further assumptions. Furthermore, a possible dependence between X and Y complicates a statistical analysis. One popular approach treats death as a dependent competing risk for disease progression. Most methods in the context of competing risks require an explicit specification of the dependent censoring structure; see Zheng and Klein (1995). However, these assumptions are not testable. Furthermore, the implicit assumption underlying many models of the existence of a hypothetical progression distribution after death is debatable. Rather than relying on artificial assumptions, an alternative is to consider only nonparametrically estimable functions, such as the descriptive measures that were proposed by Pepe and Mori (1993). These quantities, however, may have no direct interpretation for the problem of interest. Another branch of research, under the context of ‘long-term survivors’ or ‘cure models’, assumes that only a proportion of subjects will experience the event of interest. Maller and Zhou (1996) provide a useful reference on the subject. This concept may be applied to model the distribution of progression. Let B be an indicator of susceptibility for progression, with B = 1 representing susceptibility and B = 0 representing immunity or cure. Assuming that X = ∞ if B = 0, Pr.X > t/ can be written as the mixture form S1 .t/ = Pr.X > t|B = 1/pb + 1 − pb ,. .1/. where pb = Pr.B = 1/. When pb < 1, X has an improper distribution. In developing inference procedures for cure models, the definition of cure plays a crucial role. For most methods in the literature, cure is not explicitly defined. Such a model is imposed when a fraction of observations do not experience the event of interest despite long-term follow-up (Li et al., 2001). The majority of references on this rely on the assumption of independent censoring, which does not account for the effect of death. For example, Maller and Zhou (1992) considered estimating S1 .t/ by the Kaplan–Meier estimator. They claimed that, when the follow-up time is sufficiently long to observe all the susceptible individuals, the tail probability of the Kaplan–Meier estimator can be used to estimate 1 − pb and hence Pr.X > t|B = 1/ can be estimated via equation (1). Maller and Zhou (1994) further proposed a test for sufficient follow-up on the basis of the distance between the largest and the largest uncensored observations. Although their analysis.

(3) Sojourn Time Distributions. 923. looks theoretically appealing, practitioners feel concerned about identifiability of the cure fraction based on nonparametric estimation of the tail probability. Alternatively regression models, which utilize covariate information, have been studied by Farewell (1982), Kuk and Chen (1992) and Taylor (1995), among others. Identifiability of these models has been examined by Li et al. (2001). When cures become observable, identifiability is no longer an issue for nonparametric analysis. In Laska and Meisner (1992), cure is observed if a subject does not develop the event of interest within a prespecified time span. Betensky and Schoenfeld (2001) proposed a cure model with random cure times. Using the second example, they defined hospital discharge as the cure event and treated death as a competing risk. The major limitation of their method is that the death and cure times are assumed to be independent of each other. 1.2. Semicompeting risks data ˜ i , Y˜ i , δix , δiy /, i = 1, . . . , n}, where X ˜ i = Xi ∧Yi ∧Ci , In this paper we consider data of the form {.X y ˜Yi = Yi ∧ Ci , δix = I.Xi Yi ∧ Ci / and δi = I.Yi Ci /. Fine et al. (2001) called such a structure ‘semicompeting risks data’ since death is a competing risk for progression but not vice versa. Semiparametric estimation procedures for studying the dependence between X and Y when X Y have been proposed by Day et al. (1997), Fine et al. (2001) and Wang (2003). Twosample comparisons on the basis of X have been studied by Lin et al. (1996) and Chang (2000). Without imposing additional assumptions, this problem is non-identifiable owing to dependent censoring by Y . Lin et al. (1996) assumed that the distributions of X for the two groups, on a logarithmic scale, follow a bivariate location–shift model, and Chang (2000) assumed a bivariate accelerated failure time model. Although neither specified the relationship between X and Y , they implicitly assumed that the dependence structures are the same for the two groups. 1.3. Main ideas We consider nonparametric estimation of the following descriptive measures for the multipath model. Let Sij .t/ be the survival function of the sojourn time between state i and state j for i, j = 1, 2, 3 and i < j. Specifically S12 .t/ = Pr.X > t|X Y/, S13 .t/ = Pr.Y > t|X > Y/ and S23 .t/ = Pr.Y − X > t|X Y/. Also define S123 .t/ = Pr.Y > t|X Y/, which is the survival function for those who will develop progression within their lifespan. Now we discuss applications of these quantities. Reconsider the problem of a two-sample comparison on X. Let Z be the group indicator, with Z = 1 representing the treatment group and Z = 0 representing the control group. Treatment could affect both the intensity of X and that of Y from state 1. The difference between Pr.X Y |Z = 1/ and Pr.X Y |Z = 0/ reveals the treatment effect on the incidence of progression under the competing risk of death. Treatment could also affect transitions between state 1 and state 2 and/or between state 1 and state 3 (as well as between states 2 and 3). The differences of Pr.X > t|X Y , Z = j/ for j = 0, 1 reflect how the treatment changes the latency distribution of the time to progression if it does occur eventually. These effects may have different scientific meanings and hence should be studied separately. Because progression cannot be observed after death, it is natural to assign X = ∞ if X > Y . Letting ∆ = I.X Y/, we can write S1 .t/ = Pr.X > t|∆ = 1/p + 1 − p,. .2/. where p = Pr.∆ = 1/. Note that expressions (1) and (2) cannot be distinguished on the basis ˜ Y˜ , δ x , δ y /, if the cure event B is not observable. It of the information on observed variables, .X, should be noted that the susceptibility model (1) is not sensible if there is any reasonable chance.

(4) 924. W. Wang. of a competing risk of Y < X. In the second example, where discharge from hospital represents the cure event, it is reasonable to set ∆ = B. Therefore the approach that is proposed here can be directly compared with that proposed by Betensky and Schoenfeld (2001). Our method has the advantage that the times to death and cure can be correlated. For the first example of leukaemia relapse, it is questionable to assign ∆ = B. If someone dies of an unrelated cause just after entering state 1, it does not make sense to assume that he or she has been cured, or that this person has no chance of progression. Therefore we may want to exclude those who die early from the cured population even though they may not experience recurrence. In such a case, we may define the cure event as 1 − B = I.∆ = 0, Y > ξ/, where ξ is a prespecified constant, and the cured rate becomes .1 − p/ S13 .ξ/. 2.. Methodology. We propose nonparametric estimators for the path probability and the duration distributions between the states. If the path indicator ∆ = I.X Y/ is known for each subject, nonparametric estimators of these quantities can be easily constructed. When censoring is present, if δ x = 1, ∆ = I.X Y/ = 1; and, if δ x = 0 and δ y = 1, ∆ = 0. The major difficulty of estimation comes from double-censored observations with δ x = δ y = 0, for which the value of ∆ is unknown, but X > C and Y > C, and C is observable. For these observations, we define the conditional path probability given the observed value of C: 1 p.c/ = Pr.X Y |X > C, Y > C, C = c/ = Pr.X ∈ [v, v + dv/, Y v/, .3/ H.c/ v>c where H.t/ = Pr.X > t, Y > t/ is the survival function of X ∧ Y = min.X, Y/. We shall propose an estimator of p.c/ and then use it in estimating p, S12 .t/ and S13 .t/. 2.1. Estimating an unknown path probability The objective is to express p.c/ as a function of estimable quantities. To simplify the derivation, we assume that Y and C are continuous random variables and X is continuous for X Y . Our first important result is the identity ˜ ∈ [v, v + dv/, Y˜ v, δ x = 1/ 1 Pr.X ˜ p.c/ = + Pr.T < X Y/ , .4/ G.v/ H.c/ v>c where G.v/ = Pr.C > v/ and T˜ = sup{t : Pr.X > t, Y > t/ Pr.C > t/ > 0}. Note that [0, T˜ ] is ˜ All the components of p.c/ are estimable nonparametrically except the last the support of X. term in the braces. Specifically the survival function of X ∧ Y , H.t/, can be estimated by the Kaplan–Meier estimator   n.  ˜i = 1/  ˜   I. X = v, δ i   i=1 ˆ H.t/ = 1− , n.  vt    ˜ i v/   I.X i=1. where δ˜ = I.X ∧ Y C/ = δ x + δ y − δ x δ y . There are two versions of the Kaplan–Meier estiˆ 1 .t/ as the estimator based on {.X˜ i , 1 − δ˜i /, i = 1, . . . , n} mator for estimating G.t/. Denote G y ˆ and G2 .t/ as the estimator based on {.Y˜ i , 1 − δi /, i = 1, . . . , n}. Now Pr.X˜ u, Y˜ v, δ x = 1/ can be estimated by the empirical function.

(5) Sojourn Time Distributions n. i=1. 925. ˜ i u, Y˜ i v, δix = 1/=n: I.X. ˆ X˜ .n/ /, We can write Pr.T˜ < X Y/ = p.T˜ / H.T˜ /. The quantity H.T˜ / can be estimated by H. ˜ i .i = 1, . . . , n/. When the support of X ∧ Y lies within [0, T˜ ], where X˜ .n/ is the largest value of X H.T˜ / = 0 and Pr.T˜ < X Y/ = 0. When the support lies outside [0, T˜ ], p.T˜ / is not identifiable and additional assumptions are required. Assuming that p.T˜ / = p, we can derive an explicit ˜ .n/ / be an estimate of p.T˜ /. estimator for p whose formula will be given later. For now, let p. ˜ X The estimator proposed for p.c/ is n I.X ˜ i = v, Y˜ i v, δix = 1/. 1 ˆ X˜ .n/ / , p.c/ ˆ = .5/ + p. ˜ X˜ .n/ / H. ˆ ˆ H.c/ n G.v/ v>c i=1 ˆ ˆ j .t/ .j = 1, 2/. where G.v/ can be either one of G Similar arguments can be applied to the function q.c/ = 1 − p.c/. Specifically we obtain 1 Pr.X˜ v,Y˜ ∈ [v, v + dv/, δ x = 0, δ y = 1/ ˜ ˜ q.c/ = + q.T / H.T / .6/ H.c/ c<v G.v/ and 1 qˆ .c/ = ˆ H.c/. n I.X ˜ i v, Y˜ i = v, δix = 0, δiy = 1/. ˆ X˜ .n/ / , + q˜ .X˜ .n/ / H. ˆ nG.v/ v>c i=1. . .7/. ˜ .n/ /. In computing p.c/ ˜ .n/ / = 1 − p. ˜ X ˆ and qˆ .c/, we follow the convention that 0=0 = 0. where q˜ .X Theoretically we only need to estimate either p.c/ or q.c/ since they sum to 1. For finite samples, the sum is usually not equal to 1. In most of our simulations, p.c/ ˆ + qˆ .c/ ≈ 1. In formulae (4) and (6), each observation is weighted by ‘the inverse probability of censoring’, namely 1=G.v/, to adjust for the censoring bias. Similar ideas have been used in other estimation problems. Note ˆ that, as the value of v increases, G.v/ approaches closer to 0 and becomes more variable. We shall see in simulations that the performance of p.c/ ˆ and qˆ .c/ grows worse as the value of c increases. ˆ 2 .t/ is a better estimator of G.t/ than G ˆ 1 .t/ is because Pr.C Y/ Pr.C X ∧ Y/, Although G ˆ 2 .t/ would produce more efficient estimators of p.c/ ˆ and qˆ .c/. it is not guaranteed that using G In deriving asymptotic properties of all the estimators proposed, we assume that the support of X∧Y lies within [0, T˜ ], which is the condition of sufficient follow-up. In Appendix A, we prove that, for each c ∈ [0, T˜ ], sup0cT˜ | p.c/ ˆ − p.c/| → 0 with probability 1 and n1=2 { p.c/ ˆ − p.c/} converges to a mean 0 normal random variable. Properties of qˆ .c/ can be established by using similar arguments. The asymptotic variance of p.c/ ˆ can be estimated by the moment-type estimator (17), but this is quite complicated and its validity also depends on the condition of sufficient follow-up. Therefore we recommend the bootstrap estimator (18) for variance estimation. 2.2. Estimation of the sojourn time distributions Now we derive the relationships between p.c/ and other quantities of interest. By elementary properties of conditional expectation, it follows that Pr.X Y/ = Pr.δ x = 1/ + Pr.X Y , δ x = δ y = 0/ ¯ A .c/, = Pr.δ x = 1/ + p.c/ dG ¯ A .c/, Pr.X > Y/ = Pr.δ x = 0, δ y = 1/ + q.c/ dG.

(6) 926. W. Wang. ¯ A .c/ = Pr.C c, X > C, Y > C/. We can write where G ˜ ∈ [u, u + du/, δ x = 1/ Pr.X S12 .t/ = 1− , ˜ u, X Y/ Pr.X ut ˜ v, Y˜ ∈ [v, v + dv/, δ x = 0, δ y = 1/ Pr.X 1− , S13 .t/ = Pr.Y˜ v, X > Y/ vt where ˜ u, X Y/ = Pr.X ˜ u, δ x = 1/ + Pr.X. c u. Pr.Y˜ v, X > Y/ = Pr.Y˜ v, δ x = 0, δ y = 1/ +. ¯ A .c/, p.c/ dG. c v. ¯ A .c/: q.c/ dG. ¯ A .c/ can be estimated by G ¯ An.c/ = Σn I.X˜ i c, δ x = δ y = 0/=n. Using the The quantity G i i i=1 plug-in principle, we obtain the following nonparametric estimators: n n. 1 y x x ˜ pˆ = I.δ = 1/ + I.δi = δi = 0/ p. ˆ Xi / , n i=1 i i=1 n n. 1 y y ˜ i/ : qˆ = I.δix = 0, δi = 1/ + I.δix = δi = 0/ qˆ .X n i=1 i=1 Letting p. ˜ X˜ .n/ / = pˆ and qˇ .˜x.n/ / = qˆ , explicit estimators of p and q are given by n. pˆ = p. ˇ X˜ .n/ / =. i=1. I.δix = 1/ +. i∈An. ˜ i /=H. ˜ i/ ˆ 1 .X ˆ X L. , ˜ .n/ / ˆ X n − NA H. n. y ˆ 2 .˜xi /=H.˜ ˆ xi / L I.δix = 0, δi = 1/ +. qˆ = qˇ .˜x.n/ / =. i∈An. i=1. .8/. ,. ˆ x.n/ / n − NA H.˜. .9/. y. where An = {i : δix = δi = 0}, NA = ˆ 1 .c/ = L ˆ 2 .c/ = L. . . n. v>c j=1 n. v>c j=1. i∈An. ˆ xi /, 1=H.˜. ˆ ˜ j = v, Y˜ j v, δix = 1/=n G.v/, I.X y. ˆ ˜ j v, δjx = 0, Y˜ j = v, δj = 1/=n G.v/: I.X. The sojourn time survival functions can be estimated as   n. x = 1/   ˜   I. X = u, δ i i   i=1 ˆ S 12 .t/ = 1− n , n. ut    ˜ i u, δix = 1/ + ˜ i u, δix = δiy = 0/ p. ˜ i/    I.X I.X ˆ X i=1. i=1. .10/.

(7) Sojourn Time Distributions.   n. x = 0, δ y = 1/   ˜   I. Y = v, δ i i i   i=1 1− n , Sˆ13 .t/ = n.  y y x x vt    ˜ ˜ ˜  I.Yi v, δi = 0, δi = 1/ + I.Yi v, δi = δi = 0/ q. ˆ Xi /  i=1. 927. .11/. i=1.   n. y     I.Y˜ i = v, δix = 1, δi = 1/   i=1 Sˆ123 .t/ = 1− n : n.  y vt     I.Y˜ i v, δix = 1/ + I.Y˜ i v, δix = δi = 0/ p. ˆ Y˜ i /  i=1. .12/. i=1. Asymptotic properties of these estimators also depend on the condition of sufficient follow-up. Consistency and asymptotic normality of pˆ and Sˆ12 .t/ are proved in Appendices B and C. Similar arguments can be applied to derive the properties of q, ˆ Sˆ13 .t/ and Sˆ123 .t/. We also suggest use of the bootstrap for variance estimation. Nonparametric estimation of S23 .t/ is a different problem since it does not involve the missing path information. Only subjects who have experienced state 2 (i.e. δ x = 1) will contain useful information for estimating the survival function of Y −X. The challenge comes from the problem of dependent censoring. Specifically, for those with δ x = 1, the larger the value of X, the higher the chance that Y − X will be censored. Using the evolution of acquired immune deficiency syndrome as an example, Wang and Wells (1998) and Lin et al. (1999) applied the weighting technique mentioned earlier to handle the effect of dependent censoring. Both methods can be directly applied to estimate S23 .t/. 2.3. Estimation under an independent censoring model Many methods in the literature on cure models consider data of the form {.Xˇ i , δˇix /, i = 1, . . . , n}, where Xˇ i = Xi ∧ Ci and δˇix = I.Xi Ci /. This data structure ignores the competing risk of death, which is inevitable, however. Confusion arises when death without progression happens before the end of the study, i.e. Y C ∧ X. When this happens, neither Xˇ nor δˇx is identifiable since X and C are not observable after death. Nevertheless, if the only form of censoring is due to the end of the study, the value of C may be known even if Y < C. Further assuming that X = ∞ for X > Y , we can set Xˇ = C and δˇx = 0. Then S1 .t/ can be estimated by the Kaplan–Meier estimator,   n. x = 1/   ˇ ˇ   I. X = u, δ i i   i=1 1− , .13/ Sˇ 1 .t/ = n.  ut      I.Xˇ i u/ i=1. and p and S12 .t/ can be estimated by pˇ = 1 − Sˇ1 .Xˇ .n/ / and Sˇ 12 .t/ = {Sˇ 1 .t/ − 1 + p}= ˇ pˇ respectively, where Xˇ .n/ is the largest observed value of Xˇ i .i = 1, . . . , n/. We shall compare these estimators with the estimators proposed for S12 .t/ and p via simulations and data analysis. Sometimes practitioners estimate S1 .t/ by   n. x = 1/   ˜   I. X = u, δ i i   1 − i=1 n S˜1 .t/ = : .14/.  ut    ˜   I.Xi u/ i=1.

(8) 928. W. Wang. Now we discuss these two estimators of S1 .t/. In equation (8), the product–limit estimator is based on estimating the hazard for XÅ , where XÅ = X provided that X Y , given that censoring has not occurred. If X Y implies that X < ∞, Sˇ 1 .t/ is a valid estimator of S1 .t/ for t ∈ [0, T˜ ]. Beyond this range, Sˇ 1 .t/ will overestimate S1 .t/. In equation (9), the product–limit estimator is based on estimating the hazard for XÅ given that neither death nor censoring has occurred. It can be shown that S˜ 1 .t/ Sˇ 1 .t/ and S˜ 1 .t/ will underestimate S1 .t/. 3.. Numerical results. 3.1. Simulation A series of Monte Carlo simulations was conducted to examine the finite sample performance of the estimators proposed. The algorithm in Prentice and Cai (1992) was used to simulate identically and independently distributed replicates {.Xi0 , Yi /, i = 1, . . . , n} from the Clayton model: Pr.X0 > x, Y > y/ = {S Å .x/1−θ + S .y/1−θ − 1}1=.1−θ/ , θ > 1, 2. 1. where S1Å .x/ = Pr.X0 > x/ = exp.−βx/ and S2 .y/ = Pr.Y > y/ = exp.−y/. The parameter θ controls the degree of association between X0 and Y and is related to Kendall’s τ , τ = .θ − 1/=.θ + 1/. To construct the progression time for susceptible individuals, let X = I.X0 Y/X0 + I.X0 > Y/M, where M is a very large number such that Pr.Y > M/ ≈ 0. The value of β controls the level of Pr.X Y/ = p. Censoring variables were generated from U.0, 6/. Two types of data, namely semicompeting risks data and the data of independent censoring, were constructed. The censoring proportions are controlled by the values of τ and β. For example, as τ = 0:5 and β = 1, Pr.δ x = 0/ ≈ 60%, Pr.δ y = 0/ ≈ 18% and Pr.δ x = δ y = 0/ ≈ 15%. In Table 1, we evaluated the finite sample performances of the estimators proposed for p.c/ at some preselected values of c for n = 100 and 200. In simulations that are not presented here, we ˆ found that the choice of G.t/ has only a little effect on the resulting estimators. The estimator ˆ ˆ 2 . In Table 2, we compare p.c/ ˆ using G1 .t/ tends to be a little less variable than that using G several nonparametric estimators of p and q. Define pˆ j and qˆj as the estimators proposed for ˆ ˆ j .t/ .j = 1, 2/ respectively. We also evaluate the na¨ıve estimator p and q using G.t/ =G n. p˜ =. n. {I.δix i=1. i=1. I.δix = 1/ y. = 1/ + I.δix = 0, δi = 1/}. ,. the ‘empirical’ estimator of p, p¯ = Σni=1 I.∆i = 1/ and pˇ = Sˇ 1 . Xˇ .n/ /. The difference between the estimators proposed for p and p¯ reflects the effect of additional estimation of p.c/ ˆ for double-censored observations. Table 2 shows that pˆ 1 is a little more efficient than pˆ 2 but the performances of qˆ1 and qˆ2 are close. The na¨ıve estimator p˜ appears to be very unstable. It is quite good when p is close to 0.5 but can be very unreliable if p is far from 0.5. In all the cases, pˇ has a larger variation than the estimators proposed. In Table 3, we present the summary statistics for the estimators of S12 .t/ and S13 .t/ at preselected grid points for n = 200. For comparison, we also evaluated the estimators that use the true value of ∆i for all i = 1, . . . , n. Specifically define   n. x = 1/   ˜   I. X = u, δ i i   1 − i=1 , .15/ S¯ 12 .t/ = n.  ut    ˜  I.Xi u, ∆i = 1/  i=1.

(9) Sojourn Time Distributions Table 1. Summary statistics of p.c/ ˆ for .X,Y / Clayton(τ D 0:5) with S1 .x/ D exp.βx/ and S2 .y/ D exp.y/† c. p(c). 0.06. 0.38. 0.11. 0.38. 0.17. 0.37. 0.22. 0.37. 0.28. 0.36. 0.33. 0.36. 0.40. 0.35. 0.44. 0.35. 0.50. 0.34. Results for β = 5/6 and the following values of n: n = 100. n = 200. 3.1 (5.9) 3.9 (6.1) 5.4 (6.4) 5.2 (6.7) 5.1 (6.9) 5.7 (7.2) 7.1 (7.5) 7.8 (7.9) 8.5 (8.2). 1.6 (4.0) 1.2 (4.2) 1.5 (4.4) 0.9 (4.6) 1.3 (4.8) 1.0 (5.0) 1.7 (5.2) 2.6 (5.5) 2.9 (5.7). p(c). 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50. Results for β = 1 and the following values of n: n = 100. n = 200. 3.7 (5.9) 3.4 (6.2) 3.1 (6.6) 3.3 (7.0) 2.8 (7.3) 2.9 (7.7) 3.0 (8.1) 3.7 (8.4) 2.5 (9.0). −0.8 (4.3) −0.5 (4.5) −0.6 (4.8) −0.9 (5.0) −1.0 (5.3) −0.0 (5.5) −0.7 (5.7) −0.6 (6.0) −0.9 (6.2). p(c). 0.64 0.65 0.66 0.66 0.67 0.68 0.69 0.69 0.70. Results for β = 5/4 and the following values of n: n = 100. n = 200. 0.8 (5.8) 0.5 (6.3) 0.5 (6.7) 1.0 (7.0) 0.4 (7.4) 1.3 (7.7) 1.3 (8.3) 0.7 (8.7) 0.9 (8.9). 0.7 (4.2) 0.3 (4.4) 0.7 (4.6) −0.4 (4.9) −0.0 (5.1) −0.0 (5.4) −0.2 (5.7) −0.6 (6.0) −0.9 (6.2). †In each cell, the first number (times 10−3 ) is the average bias and the number in parentheses (times 10−2 ) is the standard deviation of the estimator.. Table 2. Summary statistics for various nonparametric estimators of p and q with n D 200† τ. β. p. p¯. pˆ 1. pˆ 2. qˆ1. qˆ2. p˜. pˇ. 0. 5/6. 0.46. 0. 1. 0:5. 0. 5/4. 0:56. 0.5. 5/6. 0.39. 0.5. 1. 0.50. 0.5. 5/4. 0.63. 0.8. 5/6. 0.25. 0.8. 1. 0.5. 0.8. 5/4. 0.78. 0.30 (3.61) −2.29 (3.64) −1.66 (3.43) −0.31 (3.56) 0.13 (3.51) 0.72 (3.57) −0.88 (2.93) 2.0 (3.64) 0.49 (2.91). 1.41 (3.74) −2.28 (3.72) 0.00 (3.58) 1.81 (3.82) −0.60 (3.81) −0.26 (3.59) 1.34 (3.16) 2.95 (4.04) −0.18 (3.05). 1.19 (3.74) −2.77 (3.73) −0.61 (3.58) 1.78 (3.82) −0.47 (3.86) −0.22 (3.65) 1.34 (3.16) 2.86 (4.06) −0.57 (3.16). −1.73 (3.74) 2.32 (3.73) 0.04 (3.58) −1.74 (3.83) 0.66 (3.80) 0.32 (3.59) −1.21 (3.16) −2.88 (4.04) 0.24 (3.05). −1.73 (3.75) 1.91 (3.73) −0.49 (3.56) −1.96 (3.86) 0.57 (3.81) 0.17 (3.60) −1.18 (3.17) −2.98 (4.03) 0.23 (3.05). 1.56 (3.70) −2.36 (3.71) −0.19 (3.55) 12.8 (3.83) −0.45 (3.78) −9.79 (3.60) 24.5 (3.37) 2.17 (3.84) −18.6 (3.27). 1.78 (3.88) −3.0 (3.73) 0.02 (3.58) −0.18 (3.93) −1.19 (4.04) 2.07 (3.94) −0.25 (3.17) 2.38 (4.28) 0.41 (3.70). †In each cell, the number (times 10−3 ) is the average bias and the number in parentheses (times 10−2 ) is the standard deviation of the estimator based on 500 simulation runs.. 929.

(10) 930. W. Wang Table 3. Summary statistics for the estimators of sojourn time survival functions with n D 200, τ D 0:5 and β D 1† S12 (t) and S13 (t). S¯ 12 (t). Sˆ 12 (t). Sˇ 12 (t). S¯ 13 (t). Sˆ 13 (t). 0.105. 0.825. 0.223. 0.686. 0.357. 0.570. 0.511. 0.469. 0.693. 0.378. 0.916. 0.295. 1.204. 0.217. 1.609. 0.143. 2.302. 0.071. −2.4 (4.0) −2.8 (4.6) −2.1 (5.1) 1.6 (5.1) 0.3 (5.0) 2.9 (4.9) 1.5 (4.3) −1.1 (3.8) −0.6 (3.0). −2.7 (4.0) −3.5 (4.7) −3.1 (5.3) 0.4 (5.4) −1.0 (5.4) 1.5 (5.3) 0.4 (4.8) −2.0 (4.2) −1.0 (3.4). −3.2 (4.1) −4.3 (4.8) −4.3 (5.4) −1.2 (5.6) −3.1 (5.7) −0.6 (5.6) −1.9 (5.2) −4.2 (4.7) −2.9 (4.0). −0.5 (3.6) −1.2 (4.5) −1.0 (4.9) −1.5 (5.0) −0.9 (4.9) −0.3 (4.7) 0.3 (4.3) 1.3 (3.8) 0.6 (3.0). −0.5 (3.7) −1.1 (4.6) −0.8 (4.9) −1.2 (5.1) −0.6 (5.1) −0.2 (5.1) 0.2 (4.7) 1.1 (4.1) 0.0 (3.2). t. †In each cell, the first number (times 10−3 ) is the average bias and the number in parentheses (times 10−2 ) is the standard deviation of the corresponding estimator..   n. x = 0, δ y = 1/   ˜   I. Y = v, δ i i i   i=1 ¯S13 .t/ = 1− : n.  vt      I.Y˜ i v, ∆i = 0/. .16/. i=1. Sˇ 12 .t/ is also more variable than Sˆ12 .t/. We also ran simulations under different levels of Pr.δ x = δ y = 0/. In general the variance of the estimators proposed increases as Pr.δ x = δ y = 0/ increases. 3.2. Real data examples The first data set was obtained from Table D.1 of Klein and Moeschberger (1997). There were 137 leukaemia patients receiving bone marrow transplants. Among these patients, 81 died with relapse of leukaemia, only two died without relapse and the remaining 54 patients were doubly censored. We define state 1, transplantation, state 2, relapse, and state 3, death. The na¨ıve estima˜ .n/ / = H.2640/ ˆ X ˆ tor of p gives p˜ = 81=83 = 0:976. Because H. = 0:336 > 0, we applied formulae ˆ =G ˆ 1 .t/. Because (8) and (9) and then obtained pˆ = 0:963 and qˆ = qˇ .2640/ = 0:023 using G.t/ pˆ ≈ 1, this data set may not be a typical example of the multipath model that is considered here. It turns out that Sˆ1 .t/ and Sˆ12 .t/ are very close. Since there were only two patients with δ x = 0 and δ y = 1, S˜ 1 .t/ and Sˆ1 .t/ are very close. The second data set is the Stanford heart transplantation data (Crowley and Hu, 1977). Define state 1, acceptance, state 2, transplantation, and state 3, death. Among 103 participants, 69 received transplants, 30 died without transplantation and only four observations were double censored. The na¨ıve estimator gives p˜ = 69=99 = 0:697. The data set provided the date of acceptance, date of transplant, date last seen and an indicator for the status of dead or alive..

(11) 931. 0.6 0.4 0.0. 0.2. survival probability. 0.8. 1.0. Sojourn Time Distributions. 0. 100. 200. 300. 400. 500. time. Fig. 2. Estimated survival probability of transplantation time based on the heart transplant data: – – – , Sˇ 1 .t/; - - - - - - - , Sˆ1 .t/; - - -, S˜1 .t/; , Sˆ12 .t/. ˜ .n/ / = H.1400/ ˆ X ˆ We computed H. = 0:02, which implies that the missing path proportion is tiny. There were two cases of p.c/ ˆ whose values exceed 1 but the values of q.c/ ˆ looked more reasonable. Therefore we use formula (7) for q.c/ ˆ and set p.c/ ˆ = 1 − q.c/. ˆ For the four observations with δ x = δ y = 0, the estimated path probabilities are q.1401/ ˆ = 0:304, q.428/ ˆ = 0:304, q.31/ ˆ = 0:386 and q.11/ ˆ = 0:361 which give pˆ = 0:71 and qˆ = 0:304. Note that there were 30 patients with δ x = 0 and δ y = 1. For these observations, we set Xˇ = C, which is the difference between the date of acceptance and the end of the study in April 1974, and δˇx = 0. The four estimated curves of Sˆ12 .t/, Sˆ1 .t/, S˜ 1 .t/ and Sˇ 1 .t/ are plotted in Fig. 2. We found that Sˇ 1 .t/ has a plateau at qˇ = 0:318, which is close to qˆ = 0:304, and hence Sˆ1 .t/ and Sˇ 1 .t/ are close. It turns out that the incorrect estimator S˜ 1 .t/ has a much lighter tail than Sˇ 1 .t/ and hence is misleadKM KM ing. In Fig. 3, we present the three estimated curves, Sˆ13 .t/, Sˆ2 .t/ and Sˆ123 .t/, where Sˆ2 .t/ denotes the Kaplan–Meier estimator of Pr.Y > t/. It is well known in the medical literature that the phenomenon Sˆ123 .t/ > Sˆ13 .t/ might be attributed to selection bias instead of the transplant effect. 4.. Concluding remarks. We used the framework of a multipath model to describe data with multiple end points. The methodology proposed was applied to cure models. It is well known that, when no cures are observable, the cured fraction is not identifiable (Li et al., 2001). Our analysis shows that the problem of non-identifiability is partly due to the competing risk of death. An important fact is that the commonly used data structure, {. Xˇ i , δˇix /, i = 1, . . . , n}, is by itself not identifiable if Pr.δ x = 0, δ y = 1/ > 0. Therefore the multipath model, which accounts for the competing risk of death, is a more natural formulation for studying the cure population. In the literature of cure models, ‘cure’ does not have a single definition. We used two examples to illustrate possible relationships between ∆ and B. When ∆ = B, both the method proposed and the Kaplan–Meier estimator can estimate the cured fraction and the latency distribution for the susceptible individuals. Our approach has several advantages over the Kaplan–Meier estimator. Although both meth-.

(12) W. Wang. 0.6 0.4 0.0. 0.2. survival probability. 0.8. 1.0. 932. 0. 500. 1000. 1500. time. Fig. 3. Estimated survival probability of death time based on the heart transplant data: - - - - - - - , Sˆ 123 .t/; KM - - -, Sˆ 2 .t/; , Sˆ 13 .t/. ods rely on the condition of sufficient follow-up, our approach is more robust if this assumption is violated. Specifically, under our setting, as long as H.T˜ / is close to 0 or the assumed value of p.T˜ / is reasonable, the bias of p.c/ ˆ is still minimal. The method proposed allows the dependence between death and progression and hence is more general than the method proposed by Betensky and Schoenfeld (2001). Furthermore our approach also provides estimators of S13 .t/ and S123 .t/, which may be useful in other contexts. For example when 1 − B = I.∆ = 0, Y > ξ/ indicates the cure status, S13 .ξ/.1 − p/ measures the cure probability. The technique that was used in estimating p and the sojourn time distributions is similar to the E-step in the EM algorithm. The missing values of the path indicator for double-censored observations are estimated by their conditional expected values, which are then used to estimate sufficient statistics for the quantities of interest. Acknowledgements The author is grateful to the Joint Editor, the referees and Dr Chen-Hsin Chen, Dr Tony Chen and Dr Adam Ding for their helpful comments on the work. The author also thanks J. Hsieh for performing the data analysis. The research for this paper was funded by National Science Council grant 86-2115-M-001-033-T. Appendix A: Asymptotic properties of p(c) ˆ Under sufficient follow-up, p.c/ = L.c/=H.c/, where Ku .du/=G.u/ L.c/ = c<u<T˜. and ˜ u, Y˜ v, δ x = 1/: Kv .u/ = Pr.X ˆ ˆ We can write p.c/ ˆ = L.c/= H.c/, where.

(13) ˆ L.c/ =. u>c. Sojourn Time Distributions. 933. ˆ ˆ u .du/=G.u/ K. and ˆ v .u/ = K. n i=1. ˜ i u, δix = 1, Y˜ i v/=n: I.X. We now prove strong consistency of p.c/. ˆ Since H.·/ is non-increasing, 1=H.c/ 1=H.T˜ / = M1 < ∞ for all 0 c < T˜ . Similarly 1=G.c/ 1=G.T˜ / = M2 < ∞ for all 0 c < T˜ . By strong consistency of the Kaplan–Meier estimators and applying the triangle inequality, one can show that, with probability 1, ˆ ˆ sup | p.c/ ˆ − p.c/| M1 sup |L.c/ − L.c/| + M1 sup |H.c/ − H.c/|:. 0c<T˜. 0c<T˜. 0c<T˜. By properties of the limit supremum function, strong consistency of the Kaplan–Meier estimators and applying integration by parts, one can show that, for any " > 0, ˆ ω/ − p.c/| " = 0, Pr lim sup ω : sup | p.c, 0cT˜. where ω is an element in the probability space. ˆ − p.c/}. We can write Now we show asymptotic normality of n1=2 { p.c/ 1 1=2 1 1=2 1 1=2 a ˆ ˆ n { p.c/ n {L.c/ n {H.c/ ˆ − p.c/} = − L.c/} − − H.c/}: p.c/ L.c/ H.c/ ˆ It can be shown that n1=2 {L.c/ − L.c/}=L.c/ =a an + bn , where ˆ v .v/ − Kv .v/} n1=2 {K 1 1=2 1 ˆ c .c/ − Kc .c/} + an = − n {K G.dv/ , L.c/ G.c/ G2 .v/ v>c ˆ n1=2 {G.v/ − G.v/} 1 Kv .dv/: bn = − L.c/ v>c G2 .v/ ˆ ˆ v .v/ − Kc .v/} and n1=2 {G.v/ − G.v/} to mean 0 Gaussian processes, By weak convergence of n1=2 {K asymptotic normality of bn and hence of n1=2 { p.c/ ˆ − p.c/} can be established. Each component of n1=2 { p.c/ ˆ − p.c/} can be further expressed explicitly as the sum of mean 0 random ˜ i , Y˜ i , δix , δiy , c/, where variables. Specifically we can write an = n−1=2 Σni=1 A.X 1 1 ˜ i c, Y˜ i c, δix = 1/ − Pr.X ˜ c, Y˜ c, δ x = 1/} ˜ i , Y˜ i , δix , δiy , c/ = − {I.X A.X L.c/ G.c/ 1 ˜ i v, δix = 1, Y˜ i v/ − Pr.X ˜ v, Y˜ v, δ x = 1/} dG.v/ : {I. X + 2 v>c G .v/ ˆ 2 .t/ = G.t/, ˆ If we use G it follows that. v n Mic .du/ 1 1 1=2 ˆ a , n {G.v/ − G.v/} = − 1=2 n i=1 0 Pr.Y˜ u/ G.v/ u. y where Mic .u/ = 0 I.Y˜ i s, δi = 0/ − su I.Y˜ i s/ Λc .ds/ and Λc .ds/ = Pr.C ∈ [s, s + ds/|C s/, and y hence bn = n−1=2 Σni=1 B.Y˜ i , δi , c/, where 1 dMic .du/ 1 y B.Y˜ i , δi , c/ = Kv .dv/: L.c/ v>c G.v/ 0<uv Pr.Y˜ u/. Also n 1 1 1=2 a ˆ − H.c/} = − 1=2 n {H.c/ H.c/ n i=1. where. 0. c. xy n 1 Mi .du/ ˜ i , δ˜i , c/, = 1=2 C.X ˜ u/ n i=1 Pr.X.

(14) 934. W. Wang ˜ i , δ˜i , c/ = C.X . xy. Mi .u/ =. 0. u. . c. 0. ˜ u/, Mi .du/=Pr.X xy. ˜ i s, δ˜i = 0/ − I.X. su. ˜ i s/Λxy .ds/ I.X. and Λxy .ds/ = Pr.X ∧ Y ∈ [s, s + ds/|X s, Y s/: We have shown that n 1 1 1=2 a ˜ i , Y˜ i , δix , δiy , c/ + B.Y˜ i , δiy , c/ + C.X ˜ i , δ˜i , c/}, n { p.c/ ˆ − p.c/} = 1=2 {A.X p.c/ n i=1. ˜ Y˜ , δ x , δ y / which converges in distribution to a mean 0 normal random variable with variance Σc = E[{A.X, ˜ 2 ]. Therefore a moment-type estimator of the asymptotic variance of p.c/ ˜ δ/} + B.Y˜ , δ y / + C.X, ˆ is pˆ 2 .c/. n i=1. ˆ X˜ i , Y˜ i , δix , δiy , c/ + B. ˆ X˜ i , δ˜i , c/}2 =n2 , ˆ Y˜ i , δiy , c/ + C. { A.. .17/. ˆ Bˆ and Cˆ are obtained by plugging in nonparametric estimates of the unknown components in where A, A, B and C respectively. The variance estimator (17) depends on the assumption of sufficientÅfollow-up. Å Å The bootstrap method provides a more flexible alternative for variance estimation. Let {.X˜i ,Y˜i , δix , δiy*/, i = 1, . . . , n} be a bootstrapped sample. The resampling procedure is repeated B times. Let pˆ bÅ .c/ be the estimator proposed for p.c/ based on the bth bootstrapped sample. The asymptotic variance of p.c/ ˆ can be estimated by B b=1. 2 { pˆ bÅ .c/ − p.c/} ˆ =.B − 1/:. .18/. Appendix B: Asymptotic properties of pˆ It follows that . . . n ¯ An .c/ − p.c/ dG ¯ A .c/. : ˆ dG | pˆ − p| I.δix = 1/=n − Pr.δ x = 1/ + p.c/ i=1. c. c. By the law of large numbers, The second term is bounded the first term converges to 0 in probability.. ¯ An.c/ and r2n = | p.c/{dG ¯ An.c/ − dG ¯ A .c/}|. By uniform by r1n + r2n , where r1n = c | p.c/ ˆ − p.c/| dG c and strong consistency of p.c/ ˆ and the bounded convergence theorem, r1n = op .1/. Applying integration ¯ An .c/, r2n = op .1/. It follows that by parts to r2n and by strong consistency of the empirical estimator G n1=2 . pˆ − p/ = p1n + p2n + p3n + op .1/, where p1n = n−1=2 p2n = n−1=2. n i=1. . n i=1. {I.δix = 1/ − Pr.δ x = 1/} ,. ˆ X ˆ X ˜ i , Y˜ i , δix , δiy , c/ + B. ˜ i , δ˜i , c/} dGA .c/, ˆ Y˜ i , δiy , c/ + C. p.c/{ A. p3n =. ¯ An .c/ − dG ¯ A .c/}√n: p.c/{dG. Asymptotic normality of p1n , p2n and p3n can be established by using standard techniques. The bootstrap method also provides a practical solution for variance estimation..

(15) Sojourn Time Distributions. 935. Appendix C: Asymptotic properties of Sˆ 12 (t) ˜ u, δ x = 1/ and R12 .u/ = We can write S12 .t/ = Πut {1 − F10 .du/=R12 .u/}, where F10 .u/ = Pr.X ˜ u, X Y/, which can be estimated by Pr.X n 1 ¯ An .c/ ˜ i u, δix = 1/ + I.X p.c/ ˆ dG Rˆ 12 .u/ = n i=1 cu ˜ i u, δix = 1/=n. Asymptotic normality of Sˆ 12 .t/ can be established by showing and Fˆ 10 .u/ = Σni=1 I.X that it is a smoothed function of Rˆ 12 .u/ and Fˆ 10 .u/, each of which is consistent and asymptotically normal.. References Andersen, P. K., Borgan, Ø., Gill, R. D. and Keiding, N. (1993) Statistical Models based on Counting Processes. New York: Springer. Betensky, R. A. and Schoenfeld, D. A. (2001) Nonparametric estimation in a cure model with random cure times. Biometrics, 57, 282–286. Chang, S. H. (2000) A two-sample comparison for multiple ordered event data. Biometrics, 56, 183–189. Crowley, J. and Hu, M. (1977) Covariance analysis of heart transplant survival data. J. Am. Statist. Ass., 72, 27–36. Day, R., Bryant, J. and Lefkopoulou, M. (1997) Adaptation of bivariate frailty models for prediction, with application to biological markers as prognostic indicators. Biometrika, 84, 45–56. Farewell, V. T. (1982) The use of mixture models for the analysis of survival data with long-term survivors. Biometrics, 38, 1041–1046. Fine, J. P., Jiang, H. and Chappell, R. (2001) On semi-competing risks data. Biometrika, 88, 907–919. Gray, R. (1994) A kernel method for incorporating information on disease progression in the analysis of survival. Biometrika, 81, 527–539. Hougaard, P. (2000) Analysis of Multivariate Survival Data. New York: Springer. Klein, J. P. and Moeschberger, M. L. (1997) Survival Analysis: Techniques for Censored and Truncated Data. New York: Springer. Kuk, A. Y. C. and Chen, C. (1992) A mixture model combining logistic regression with proportional hazards regressions. Biometrika, 79, 531–541. Laska, E. M. and Meisner, M. J. (1992) Nonparametric estimation and testing in a cure model. Biometrics, 48, 1223–1234. Li, C.-S., Taylor, J. M. G. and Sy, J. P. (2001) Identifiability of cure models. Statist. Probab. Lett., 54, 389–395. Lin, D. Y., Robins, J. M. and Wei, L. J. (1996) Comparing two failure time distributions in the presence of dependent censoring. Biometrika, 83, 381–393. Lin, D. Y., Sun, W. and Ying, Z. (1999) Nonparametric estimation of the gap time distributions for serial events with censored data. Biometrika, 86, 59–70. Maller, R. A. and Zhou, S. (1992) Estimating the proportion of immunes in a censored sample. Biometrika, 79, 731–739. Maller, R. A. and Zhou, S. (1994) Testing for sufficient follow-up and outliers in survival data. J. Am. Statist. Ass., 89, 1499–1506. Maller, R. A. and Zhou, S. (1996) Survival Analysis with Long-term Survivors. New York: Wiley. Pepe, M. S. and Mori, M. (1993) Kaplan–Meier, marginal or conditional probability curves in summarizing competing risks failure time data? Statist. Med., 12, 737–751. Prentice, R. L. and Cai, J. (1992) Covariance and survival function estimation using censored multivariate failure time data. Biometrika, 79, 495–512. Prentice, R. L., Kalbfleisch, J. D., Peterson, A. V., Flournoy, N., Farewell, V. T. and Breslow, N. E. (1978) The analysis of failure times in the presence of competing risks. Biometrics, 34, 541–554. Taylor, J. M. G. (1995) Semi-parametric estimation in failure time mixture models. Biometrics, 51, 899–907. Wang, W. (2003) Estimating the association parameter for copula models under dependent censoring. J. R. Statist. Soc. B, 65, 257–273. Wang, W. and Wells, M. T. (1998) Nonparametric estimation of successive duration times under dependent censoring. Biometrika, 85, 561–572. Zheng, M. and Klein, J. P. (1995) Estimates of marginal survival for dependent competing risks based on an assumed copula. Biometrika, 82, 127–138..

(16)