國
立
交
通
大
學
統計學研究所
碩
士
論
文
非混合治癒模型之統計推論
Statistical Inference for Cure Models
The Non-mixture Approach
研 究 生： 陳瀅竹
指導教授： 王維菁 教授
非混合治癒模型之統計推論
Statistical Inference for Cure Models
The Non-mixture Approach
研 究 生：陳瀅竹 Student：
Ying-jhu Chen
指導教授：王維菁 教授 Advisor：
Dr. Weijing Wang
國 立 交 通 大 學
統 計 學 研 究 所
碩 士 論 文
A Thesis
Submitted to Institute of Statistics College of Science
National Chiao Tung University in partial Fulfillment of the Requirements
for the Degree of Master in
Statistics
June,2009
Hsinchu, Taiwan, Republic of China
非混合治癒模型之統計推論
學生：陳瀅竹
指導教授：王維菁 教授
國立交通大學統計學研究所 碩士班
摘
要
在存活分析中，若在實驗期間能無限延長下，有部分的資料能夠
永久存活或者不會經歷到所感興趣的事件發生，必須考慮使用治癒模
型。本論文主要探討的是 Tsodikov (1998) 所提出的非混合治癒模型，
討論在數學模型的架構下代表的生物意義，藉模擬討論在有母數與無
母數的統計推論方法中對治癒率以及發病時間函數參數估計的影響。
最後在有解釋變數對模型的影響中，提供符合比例風險迴歸模型(PH)
與加速失敗時間模式(accelerated failure time model)的圖型診斷方法。
Statistical Inference for Cure Models
The Non-Mixture Approach
Student: Ying-jhu Chen
Advisor: Dr. Weijing Wang
Institute of Statistics National Chiao Tung University
Hsinchu, Taiwan
Abstract
Cure models are survival models which allows for the existence of cure despite of
long-term follow-up. In this thesis, we consider statistic inference for cure models
based on the non-mixture approach proposed by Tsodikov(1998). This model has some interesting biological interpretations. Parametric and semi-parametric methods will be investigated. We propose a model diagnostic tool to verify the form of regression models.
誌
謝
兩年的碩士生涯，在此論文完成同時也響起了尾聲。
很幸運的在兩年前能夠進入學術界神聖的殿堂，衷心感謝交大統
計所老師這些日子的教導。尤其要感謝我的恩師王維菁教授，您總是
細心的帶領著我，無論是在課業或論文遇到困難時，引導我如何剖析
問題並一步一步去解決。您給了我一個不一樣的視野，使我看問題的
核心更加準確。而對於我的生活及未來規劃，也會適時地給予意見與
協助。
感謝研究所班上二十一位同學們。這段全天候打拼的日子，沒有
你們的鼓勵還有支持，我不會成長了那麼多。一同討論問題、腦力激
盪找原因，過程中偶而參插著歡樂的笑聲，我永遠會記著這每一刻。
由衷的感謝你們，讓我在交大統研所的生活更多采多姿。
特別要感謝東華大學應用數學系的曾玉玲教授。是您，我才會踏
入統計這個領域。我的統計生涯，因為有您，開始有了起點。在這畢
業之際，對於未來，我會繼續努力。
最後，謹將此論文獻給我最愛的家人還有王維菁教授，謝謝你們
無私的奉獻，給我這個機會支持我，讓我順利的完成碩士學歷。
感謝這一路上曾經幫助過我的人，你們都是我生命中無法取代的
貴人。
陳瀅竹 謹誌于
國立交通大學統計學研究所
中華民國九十八年六月
Contents
1 Introduction 1
1.1 Background . . . 1
1.2 The Mixture Approach . . . 2
1.3 The Non-Mixture Approach . . . 2
1.4 Comparison of the Two Models . . . 3
1.5 Estimation under Censoring . . . 3
1.6 Outline of the Thesis . . . 4
2 Model Properties 5 2.1 Biological Interpretation . . . 5 2.2 Mathematical Properties . . . 6 3 Parametric Inference 9 3.1 Likelihood Function . . . 9 3.2 Two Examples . . . 10
3.3 Generation of a Parametric Distribution . . . 10
3.3.1 Algorithm I . . . 11 3.3.2 Algorithm II . . . 11 3.4 Simulation Results . . . 12 4 Nonparametric Inference 13 4.1 Likelihood Function . . . 13 4.2 Score Functions . . . 14
4.2.1 An algorithm based on the method of Lagrange Multiplier . . . 14
4.2.2 An algorithm based on change of variables . . . 16
4.3 Simulation Results . . . 18
5 Regression Analysis 19 5.1 Classification of Covariate Effects . . . 19
5.3 Short Term Effect: Accelerated Failure Time Model . . . 23
6 Conclusion 25
Reference 26
Appendix: Additional Figures 27
List of Tables
A-1 Relationship between C∼Unif[0,k] and Pr(δ = 0) for Gamma Distributions . 37 A-2 Relationship between C∼Unif[0,k] and Pr(δ = 0) for Weibull Distributions . 37 A-3 Relationship between C∼Unif[0,k] and Pr(δ = 0) for Log-normal Distributions 37 A-4 Maximized likelihood estimators for Gamma distributions with p = 0.3 and
Pr(δ = 0) = 0.3 . . . 38 A-5 Maximized likelihood estimators for Weibull distributions with p = 0.3 and
Pr(δ = 0) = 0.3 . . . 38 A-6 Maximized likelihood estimators for Log-normal distributions with p = 0.3
and Pr(δ = 0) = 0.3 . . . 39 A-7 Maximized likelihood estimators for Gamma distributions with p = 0.3 and
Pr(δ = 0) = 0.4 . . . 40 A-8 Maximized likelihood estimators for Weibull distributions with p = 0.3 and
Pr(δ = 0) = 0.4 . . . 40 A-9 Maximized likelihood estimators for Log-normal distributions with p = 0.3
and Pr(δ = 0) = 0.4 . . . 41 A-10 NPMLE of θ and p for Gamma distributions with p = 0.3 and Pr(δ = 0) = 0.3. 42 A-11 NPMLE of θ and p for Weibull distributions with p = 0.3 and Pr(δ = 0) = 0.3. 42 A-12 NPMLE of θ and p for Log-normal distributions with p = 0.3 and Pr(δ = 0) =
0.3. . . 42 A-13 NPMLE of θ and p for Gamma distributions with p = 0.3 and Pr(δ = 0) = 0.4. 43 A-14 NPMLE of θ and p for Weibull distributions with p = 0.3 and Pr(δ = 0) = 0.4. 43 A-15 NPMLE of θ and p for Log-normal distributions with p = 0.3 and Pr(δ = 0) =
List of Figures
1-1 K-M Survival function for patients with Hodgkin’s disease (Tsodikov, 2003) . 1
2-1 The Number of Recurrence tumor . . . 6
2-2 The failure time T given tumor number N . . . 8
4-1 Estimated survival funtions under model mis-specification . . . 18
5-1 The effects of covariate . . . 20
5-2 K-M curves and diagnostic plots when the PH assumption holds. . . 22
5-3 K-M curves and diagnostic plots when the PH assumption is violated. . . 22
5-4 K-M curves and diagnostic plots when the AFT assumption holds. . . 24
5-5 K-M curves and diagnostic plots when the AFT assumption is violated. . . . 24
A-1 Survival density and hazard functions for selected parametric distributions. . 27
A-2 Estimated survival functions when ˜F is correctly specified as Gamma(1,1) . 28 A-3 Estimated survival functions when ˜F is correctly specified as Gamma(4,1.5) 29 A-4 Estimated survival functions when ˜F is correctly specified as Gamma(6,2) . 30 A-5 Estimated survival functions when ˜F is correctly specified as Weibull(1,3) . . 31 A-6 Estimated survival functions when ˜F is correctly specified as Weibull(2,10) . 32 A-7 Estimated survival functions when ˜F is correctly specified as Weibull(4,15) . 33 A-8 Estimated survival functions when ˜F is correctly specified as log-normal(1,1) 34 A-9 Estimated survival functions when ˜F is correctly specified as log-normal(2,0.4) 35 A-10 Estimated survival functions when ˜F is correctly specified as log-normal(2.5,0.25) 36
Chapter 1
Introduction
1.1
Background
Let T be the failure time of interest which is subject to right censoring by C. Denote S(t) = Pr(T > t), h(t) = lim∆→0Pr(T ∈[t,t+∆)|T ≥t)∆ and H(t) =
Rt
x=0h(x)dx as the survival
hazard and cumulative hazard functions respectively. Observed data can be denoted as {(Xi, δi)(i = 1, . . . , n)}, where Xi = Ti ∧ Ci and δi = I(Ti ≤ Ci) and (Ti, Ci) are random
replications of (T, C). Let t(1) < t(2) < · · · < t(D) be distinct ordered failure times. Under
the assumption that Ti and Ci are independent, S(t) can be estimated by the Kaplan-Meier
estimator ˆ S(t) = Y k:t(k)<t 1 −dk rk ,
where dk is the number of failures and rk is the number at risk at t(k).
Figure 1-1 depicts the Kaplan-Meier curve of an example in the paper of Tsodikov(2003). The curve estimates the survival function for patients with Hodgkin’s disease treated by radiotherapy. The event of interest was defined as death due to the disease. An interesting
1064 Journal of the American Statistical Association, December 2003
(a) (b)
Figure 1. (a) Relapse-Free Survival for Patients With Hodgkin’s Disease Treated by Radiotherapy. Data from the International Database on Hodgkin’s Disease. (b) Breast cancer specic survival in local stage by age. Data from the Surveillance Epidemiology and End Results Program.
where M is a binary random variable taking on the values of 0 and 1 with probability p and 1¡ p, respectively, with
pD PrfX D 1g;
and G0.t/ is de ned as the survival function for the time to
failure conditional upon ultimate failure, that is,
G0.t/D PrfX ¸ tjX < 1g: (1.3)
When designing regression counterparts of model (1.2), it is common practice to use the logistic regression model for in-corporating covariates into the probability p; and proportional hazards regression for modeling the effect of the covariates on the conditional survival function G0.
An alternative, but equally general, representation of an im-proper survival time distribution can be obtained by assuming that the cumulative hazard 3.t/DR0t¸.t/ dt has a nite posi-tive limit, say µ , as t tends to1. In this case one can write
G.t/D e¡µF .t/; µ > 0; t ¸ 0; (1.4)
where F .t/D 3.t/=µ is the c.d.f. of some nonnegativerandom variable such that F .0/D 0. In what follows we will call the model given by (1.4) the bounded cumulative hazard (BCH) model.
Within the nonparametric framework it makes no difference whether representation (1.2) or (1.4) is used as a basis for the estimation of p, but the situation is not the same when G0.t/ is
parametrically speci ed. Using de nition (1.3), one can repre-sent the survival function (1.4) in the form of formula (1.2), but both p and G0become functions of the common parameter µ :
pD e¡µ; G0.t/D e
¡µF .t/¡ e¡µ
1¡ e¡µ :
A similar confounding occurs when one represents (1.2) in the form of (1.4). A summary of useful formulas showing the rela-tionship between the two models was given in Chen, Ibrahim, and Sinha (1999).
Although the models given by (1.2) and (1.4) are just two different ways of rescaling the survival function G.t/, it took
will be discussed later), which is the main subject matter of the present article. The rst move in this direction was due to Hay-bittle (1959, 1965). The author proceeded from the observation that in some clinical data on cancer survival, an actuarial es-timate of the hazard function tends to decrease exponentially with time. If the same property holds for the true hazard, ex-pression (1.4) assumes the form:
G.t/D expf¡µ.1 ¡ e¡³t/g; ³ > 0: (1.5) A comprehensive treatment of this model was given by Cantor and Shuster (1992) who used the following parameterization:
G.t/D exp » ¯ ³.1¡ e ³ t/ ¼ ; ³ 6D 0; (1.6)
where ¯ > 0, but ³ may take either sign. If ³ > 0 the survival function (1.6) corresponds to the proper Gompertz distribution, but if ³ < 0 the distribution is improper with the surviving function equaling e¯=³. For this reason the distribution given
by (1.6) can be called a generalized Gompertz distribution. In more recent work Cantor (2001) used representation (1.6) to de-termine the projected variance of estimated survival probabili-ties in clinical trials. A generalization of the model (1.6) was proposed by Cantor (1997) who suggested approximating a log hazard by a polynomial to obtain an estimate of the cure rate based on formula (1.1).
Clearly, the Gompertz-like model given by formula (1.5) is a special case of formula (1.4) with the function F .t/ spec-i ed by an exponential c.d.f. with parameter ³ . The repre-sentation (1.4) was rst introduced in an article of Yakovlev et al. (1993) and discussed later as an alternative to the mix-ture model by Yakovlev (1994). Interestingly enough, Yakovlev et al. (1993) proceeded from purely biological considerations; the idea of imposing a constraint on the behavior of the haz-ard function was introduced later in Tsodikov, Loef er, and Yakovlev (1998b). In fact, the authors proposed a simple mech-anistically motivated model of tumor recurrence yielding an improper survival time distribution. Under this model the prob-ability of tumor cure is de ned as the probprob-ability of no
clono-Figure 1-1: K-M Survival function for patients with Hodgkin’s disease (Tsodikov, 2003) feature of the plot is that the tail does not descend to zero. One explanation is that the study period is not long enough to observe large failure points. A different interpretation
is that some patients would never develop the event of interest despite of long-term follow-up. In other words, these patients can be viewed as ”cured” and hence free from the event of interest. In many practical applications, such a plateau in the Kaplan-Meier curves is commonly seen. Estimating the proportion of cure may have important medical implication. Cure models are survival models that allow for cured individuals (Boag, 1949). Here we introduce two types of models which include the possibility of cure in survival analysis.
1.2
The Mixture Approach
It is assumed that the population can be divided into two sub-populations, namely the susceptible group and the cured group. For people in the cured group (with proportion p), the failure time is set to be ∞. Subjects in the susceptible group (with proportion 1 − p) will ultimately experience the event of interest. Denote the corresponding survival function as S0(t) which is a proper survival function satisfying limt→∞S0(t) = 0. Hence the survival
function of the population can be written as the following mixture form,
S(t) = p + (1 − p)S0(t). (1)
The corresponding hazard function is given by
h(t) = (1 − p)f (t) S(t) ,
where f (t) is the density of T . This model was developed by Berkson and Gage (1952) and later studied extensively by a number of statisticians.
1.3
The Non-Mixture Approach
In the thesis, we will focus on the second type of cure fraction model. Yakovlev et al. (1993) proposed the so-called ”non-mixture cure model” or ”bounded cumulative hazard” (BCH) model by defining an asymptote for the cumulative hazard function. Specifically it is assumed that the cumulative hazard function can be written as
H(t) = Z t
0
where 0 ≤ ˜F (t) ≤ 1 satisfies properties of a cumulative distribution function. It follows that lim t→∞H(t) = θ < ∞. Since S(t) = exp (−H(t)) = exp−θ ˜F(t), (2) it follows that lim t→∞S(t) = exp(−θ) = p > 0.
1.4
Comparison of the Two Models
The survival function in (2) can be re-expressed as the mathematical form of model (1). Specifically define p = e−θ and then S0(t) = e−θ ˜F (t)− e−θ 1 − e−θ . (3)
Notice, that in (3), p and S0(t) are both functions of θ. If the cure fraction p is the only
interest, it makes no difference which model formulation is chosen under a nonparametric framework. However, an important feature of the mixture model in (1) is that p and S0(t)
are modeled separately so that S0(t) does not contain information of p. Therefore, the two
approaches become different when additional model assumptions are imposed on the cure rate and the latency distributions.
1.5
Estimation under Censoring
In practice, censoring happens when subjects are lost to follow-up. Observed data are of the form, {(Xi, δi), i = 1, 2, · · · , n}. Our main purpose is to estimate θ and ˜F (·). It is easy
to see that Pr(δ = 0) ≥ exp(−θ) and Pr(δ = 0) = exp(−θ) if C = ∞. Estimation of θ will be affected by both the distribution of C and ˜F (t). In particular, those outliers are those with extreme values of T which are likely to be censored.
1.6
Outline of the Thesis
The thesis considers statistical inference for the non-mixture model in (2). In Chapter 2, we provide more discussions on this model which include useful biological interpretations. Parametric and non-parametric inference methods are studied in Chapter 3 and Chapter 4 respectively. Regression analysis that models the effect of covariates on p and ˜F is discussed in Chapter 5. We also propose diagnostic plots to verify the proportional hazard (PH) and accelerated failure time (AFT) assumptions. Concluding remarks are given in Chapter 6.
Chapter 2
Model Properties
In this chapter, we will examine properties of the non-mixture model: S(t) = exp−θ ˜F (t).
In Section 2.1, we introduce a model with interesting biological interpretations which would yield the above representation. More detailed mathematical properties are examined in Section 2.2.
2.1
Biological Interpretation
Consider the following biological mechanism about tumor development. Assume that an individual in the population has N tumor cells where N = 0, 1, 2 · · · . Each tumor cell will develop the event of the metastasis with the failure time ˜T . Let ˜Tj be the random time for the
jthtumor cell to produce the detectable metastatic disease. Cancer occurs under a competing
risk framework such that the time to develop cancer is defined as T = min{ ˜Tj, 1 ≤ j ≤ N }
if N ≥ 1. If a person has no tumor cell with N = 0, he/she is called ”immune” or ”cured”. We will set Pr( ˜T0 = ∞) = 1.
Assume that N follows a Poisson distribution with mean θ. Also assume that, Given N , the random variables ˜Tj (j = 1, . . . , N ) are independent and identically distributed with a
common distribution function ˜F (t) = Pr(Tj ≤ t) that does not depend on N . It follows that
S(t), the survival function for T = min{ ˜Tj, 0 ≤ j ≤ N }, can be written as
S(t) = Pr(no metastatic cancer by time t) = Pr(N = 0) + ∞ X k=1 Pr(T > t, N = k) = Pr(N = 0) + ∞ X k=1 Pr(T > t)Pr(N = k) = exp (−θ) + ∞ X k=1 ( ˜S(t))kexp(−θ)θ k k! = exp (−θ + θ ˜S(t)) = exp (−θ ˜F (t)),
which is exactly the model considered here. There are two distinct characteristics of tumor growth; namely the initial number of tumor cells and the rate of their progression. Thus parameters in the model have a clear biological meaning. Such a model is reasonable to describe cancer relapse, time to inflection and so forth.
Notice that lim
t→∞S(t) = Pr(N = 0) = exp(−θ). As θ → ∞, the cure rate tends to
zero; as θ → 0, the cure rate reduces to one. Usually we have 0 < p = exp(−θ) < 1. Since H(t) = θ ˜F (t), the survival function is S(t) = e−θ ˜F (t) and the hazard function becomes h(t) = θ ˜f (t), where that ˜f is the density of Tj.
2.2
Mathematical Properties
Based on the biological framework, we discuss how the shape of S(t) is affected by N and ˜F (t). Notice that the cure rate equals Pr(N = 0) = exp(−θ) = p which measures the heaviness of the tail of S(t). The failure time T is determined by both N and the distribution of Tj for j = 1, . . . , N . Larger value of N (i.e. more tumor cells) tends to be associated with
smaller value of T = min{ ˜Tj, 1 ≤ j ≤ N } (i.e. early onset of the event). In other words, if a
person has more tumor cells, the time of cancer relapse tends to happen earlier.
N Pr(N=n) 0 2 4 6 8 10 0.0 0.2 0.4 0.6 0.8 1.0 N~Poisson( 2.3 ) N Pr(N=n) 0 2 4 6 8 0.0 0.2 0.4 0.6 0.8 1.0 N~Poisson( 1.2 ) N Pr(N=n) 0 1 2 3 4 5 0.0 0.2 0.4 0.6 0.8 1.0 N~Poisson( 0.69 ) N Pr(N=n) 0 1 2 3 4 0.0 0.2 0.4 0.6 0.8 1.0 N~Poisson( 0.36 )
Figure 2-1: The Number of Recurrence tumor
Figure 2-1 plots the probability function of N for selected values of θ. The shaded region corresponds to the cure rate, Pr(N = 0). Figure 2-2 contains several boxplots of T |N
for N = 0, 1, 2, · · · where ˜Tj follows the Gamma, Weibull and Log-normal distributions
respectively. As expected, the median of T |N shows a decreasing pattern as a function of N . The pattern of outliers in each plot deserves special attention. Outliers of T happen more frequently for those with N = 1. Especially in Figure 2-2, there are has a lot of outliers for T |N = 1 with Tj ∼ Weibull(1,3). If the follow-up period is not very long, large outliers
are more likely to be censored. It becomes more difficult to estimate p = exp(−θ) if T has a long tail. We will investigate this conjecture via simulations. Among the three Weibull distributions, Var(T |N ) for N ≥ 1 is largest for Weibull(4,15). We will also examine how this characteristic affects parametric inference.
In Figure A-1 of the Appendix, we plot the survival function of T , and the density and the hazard functions of Tj for selected parametric distributions.
1 2 3 4 5 01234567 N T Gamma ( 1 , 1 ) 1 2 3 4 5 6 246 N T Gamma ( 4 , 1.5 ) 1 2 3 4 5 2468 N T Gamma ( 6 , 2 ) (a) Gamma 1 2 3 4 5 0 5 10 15 N T Weibull ( 1 , 3 ) 1 2 3 4 5 6 0 5 10 15 20 N T Weibull ( 2 , 10 ) 1 2 3 4 5 7 5 1 01 52 0 2 5 N T Weibull ( 4 , 15 ) (b) Weibull 1 2 3 4 5 6 0 1 02 03 04 05 0 N T Log-normal ( 1 , 1 ) 1 2 3 4 5 6 5 1 01 5 2 02 5 3 0 N T Log-normal ( 2 , 0.4 ) 1 2 3 4 5 6 5 1 01 5 2 02 5 N T Log-normal ( 2.5 , 0.25 ) (c) Log-normal
Chapter 3
Parametric Inference
Parametric analysis can provide more clear description of the data if the model assumption is correct. In Section 3.1, we derive the likelihood and score functions under a parametric framework. Specific examples are given in Section 3.2. Section 3.3 presents numerical analysis that compares the bias and standard deviation of the maximum likelihood estimates under several parametric settings.
3.1
Likelihood Function
Recall that T denotes the failure time with the survival function S(t) = Pr(T > t). The censoring variable is denoted by C with the survival function G(t) = Pr(C > t) and the density function g(t). Assume that T and C are independent. In presence of right censoring, the observed variables are (X, δ), where X = T ∧ C and δ = I(T ≤ C). The sample contains independently and identically distributed observations of (X, δ), denoted as (Xi, δi)(i = 1, 2, · · · , n).
The likelihood function has the form
n Y i=1 [f (xi)G(xi)]δi × [S(xi)g(xi)]1−δi ∝ n Y i=1 [f (xi)δiS(xi)1−δi].
Under the non-mixture model, S(t) can be express as
S(t) = exp(−θ ˜F (t)), (4)
and the density function can be expressed as f (t) = −dS(t)
dt = θ ˜f (t) exp(−θ ˜F (t)), (5) where ˜f (t) = ˜F0(t). Accordingly the parametric log-likelihood function can be written as
n
X
i=1
{δi[log(θ) + log ˜f (xi)] − θ ˜F (xi)}. (6)
Suppose that ˜F (·) is modeled as a parametric function ˜Fη(·) so that the log-likelihood is
a function of (θ, η), say `(θ, η). The maximum likelihood estimator of (θ, η) can be obtained by solving the equations ∂θ∂`(θ, η) = 0 and ∂η∂`(θ, η) = 0.
3.2
Two Examples
First suppose that ˜f (·) follows an exponential distribution with parameter α. The log-likelihood function in (6) becomes
`(θ, α) ∝
n
X
i=1
{δilog(θαe−αxi) − θ(1 − e−αxi)}.
The score functions are θ = Pn i=1δi Pn
i=1(1−e−αxi)
, Pn i=1{δi( 1 α − xi) − θxie −αxi} = 0.
The solution ˆα and ˆθ can be obtained by numerical methods.
In the second example, we assume that ˜f (·) follows a Weibull distribution with the param-eters k and λ, where k is the shape parameter and λ is the scale parameter. The decreasing rate of hazard is controlled by k. If k > 1, the hazard he hazard has a peak in the λ; if k ≤ 1, it is decreasing. Note that the shape of hazard function has right tail when coefficient of skewness is larger than zero. The log-likelihood function in (6) becomes
n X i=1 δi
log θ + log k + (k − 1) log xi− (k − 1) log λ −
xi λ k − θ1 − e(xiλ) k . The score functions are
θ = Pn i=1δi Pn i=11−e( xi λ) k, Pn i=1 n δi h 1 k + log xi− log λ − xi λ k log xi λ i + θ xi λ k log xi λe( xi λ) ko = 0, Pn i=1 n δi k−1 λ + kλ −k−1 − θkλ−k−1e(xiλ)ko= 0.
The maximum likelihood estimates, ˆk, ˆλ and ˆθ, can be obtained by numerical methods.
3.3
Generation of a Parametric Distribution
To generate data from a parametric model of the form S(t) = exp(−θ ˜F (t)), we propose two algorithms. One is based on mathematical properties of the random variables and the other utilizes properties of the biological model as discussed earlier.
3.3.1 Algorithm I
This algorithm is constructed based on the biological model of tumor development. • Step 1: Generate the number of tumor cell N from a Poisson distribution with mean θ
such that θ = − log(p) , where p is the cure fraction. • Step 2: If N = 0, set T to be a large number;
if N > 0, generate ˜Tj from ˜F (·) independently for j = 1, . . . , N .
• Step 3: Set T = min{ ˜T1, ˜T2, · · · , ˜TN}.
• Step 4: Generate C from U nif (0, K).
• Step 5: Set X = min{T, C} and δ = I(T ≤ C). 3.3.2 Algorithm II
Without the biological setting, we has to handel the challenge that T is not a proper random variable since positive mass is put on T = ∞. The second algorithm uses the idea of mixture models to generate random variables. Specifically under the mixture model, we have
S(t) = (1 − p)S0(t) + p,
where S0(t) is the survival function of a proper random variable, say T0. It follows that
S0(T0) ∼ U (0, 1) so that T0 = S0−1(U ) where U ∼ U (0, 1). Hence we can set T = T0 with
probability 1 − p and T as a very large number with probability p.
The relationship between T0 and ( ˜F (·), θ) is given in equation (3). Thus we have
S0(T0) =
e−θ ˜F (T0)− e−θ
1 − e−θ . (7)
Suppose that the parametric form of ˜F (·) is specified. It follows that ˜
F (T0) =
− log{U (1 − e−θ) + e−θ}
θ .
The second algorithm is described below. • Step 1: Generate U ∼ U(0, 1).
• Step 2: Set T0 = ˜F−1 −log(U {1 − e −θ} + e−θ θ . • Step 3: Generate V ∼ U(0, 1).
If V > p = exp(−θ), set T = T0;
if V < p, set T equal to a very large number. • Step 4: Generate C from a uniform distribution. • Step 5: Set X = T ∧ C and δ = I(T ≤ C).
3.4
Simulation Results
In the simulations, we evaluate the MLE performance for selected parametric distribu-tions. One set if simulations did not include censoring by setting C = ∞ with p = 0.3. The other setting was designed such that p = 0.3 and Pr(δ = 0) = 0.4. This means that 10% subjects who are censored but actually susceptible. Since C follows a U(0,k), we need to find the value of k such that Pr(T ≤ C) can achieve the target goal. Table ?? lists the value of k which yields the target value Pr(δ = 0) for different parametric distributions.
Data were generated using the algorithm derived from model. Two sample size with n = 100 and n = 300 were considered. For each estimator, the average bias and standard deviation of parameter estimates are reported base on 500 replications.
Table A-4 ∼ A-6 summarize the result of MLE under C = ∞ when Tj follows Gamma,
Weibulll and log-normal distributions respectively. Note that p = exp(−θ) = 0.3. Table A-7 ∼ A-9 contain the results with same p = 0.3 but Pr(δ = 0) = 0.4. That is 10% subjects who have at least one tumors but are censored due to limited observation period. The estimate of p is poor for Weibull(1,3) which has more outliers than the other two Weibull distributions. The estimator of λ is poor for Weibull(1,3) and Weibull(4,15). Parameter λ controls the spread of the distribution which involves the tail behavior. Weibull(4,15) has large tail values which are likely to be censored.
Chapter 4
Nonparametric Inference
In this chapter, we discuss estimation of θ and ˜F (·) when the distribution of the latter is not specified. The paper of Tsodikov (2001) considered this setting but did not provide the details. Here we give detailed derivations. The likelihood function and score functions are discussed first. Two algorithms for solving the score equations are proposed.
4.1
Likelihood Function
Based on the original data notations, the nonparametric likelihood function can be written as
n
Y
i=1
[Pr(T = xi, C > xi)]δi[Pr(C = xi, T > xi)]1−δi .
After taking logarithm and assuming that Ti and Ci are independent, the log-likelihood
function can be written as
n
X
i=1
{δilog[Pr(T = xi)] + (1 − δi) log[Pr(T > xi)]} .
In order to estimate ˜F (·) nonparametrically, we assume that is takes jumps only in observed failure points. Now, we express the log-likelihood function in terms of ordered observations. Let t(1) < t(2) < . . . < t(D) be observed failure points. Set t(0) = 0 and t(D+1) = ∞. Denote
mj =
Pn
i=1I(Xi = t(j), δi = 1) which measures the number of failure points at time t(j) and
and nj =
Pn
i=1I(t(j)≤ Xi < t(j+1), δi = 0) which measures the number of censored points in
the interval [t(j), t(j+1)). Under the model assumption, it follows that
log L(θ, ˜F ) ∝
D
X
i=1
milog[S(t(i−)) − S(t(i))] + nilog[S(t(i))]
= D X i=1 ( mi " −θ i−1 X k=1 4 ˜Fk+ log h 1 − exp −θ4 ˜Fi i # + ni " −θ i X k=1 4 ˜Fk #) = −θ D X i=2 i−1 X k=1 mi4 ˜Fk+ D X i=1 milog h 1 − exp−θ4 ˜Fi i − θ D X i=1 i X k=1 ni4 ˜Fk
where 4 ˜Fi = ˜F (t(i))− ˜F (t(i−1)) subject to the restriction that
PD
i=14 ˜Fi = 1. Now we discuss
4.2
Score Functions
In this section, we present two ways to solve the score functions. In Section 4.2.1, since we want to maximize log-likelihood function subject to some constraints, the lagrange multiplier method is natural choice. The other way to solve the variables is presented in Section 4.2.2. We change the form of variables in the log-likelihood function and use the condition to solve the variables.
4.2.1 An algorithm based on the method of Lagrange Multiplier
Applying the Lagrange multiplier method, we can handle the constraint of the parameters in the maximization procedure. Adding the Lagrange multiplier λ, we consider the function
h(θ, ˜F , λ) = log L(θ, ˜F ) + λ(1 −
D
X
i=1
4 ˜Fi).
Taking the derivatives with respect to the unknown parameters, the score equations can be written as ∂h(θ, ˜F , λ) ∂4 ˜FD = 0 ⇔ λ = ∂ log L(θ, ˜F ) ∂4 ˜FD = mDθ exp(−θ4 ˜FD) 1 − exp(−θ4 ˜FD) − θnD (8) and ∂h(θ, ˜F , λ) ∂4 ˜Fj = 0, 1 ≤ j ≤ D − 1. It follows that for each j = 1, 2, · · · , D − 1, we have
− θ D X i=j+1 mi− θ D X i=j ni+ mjθ exp(−θ4 ˜Fj) 1 − exp(−θ4 ˜Fj) − λ = 0. (9) By taking ∂h(θ, ˜F , λ) ∂θ = 0, we obtain − D X i=2 i−1 X k=1 mi4 ˜Fk− D X i=1 i X k=1 ni4 ˜Fk+ D X i=1 mi 4 ˜Fiexp −θ4 ˜Fi 1 − exp−θ4 ˜Fi = 0 (10) and also D X i=1 4 ˜Fi = 1. (11)
• Step 1: Take summation (9)×4 ˜Fi from i =1 to D-1 − θ D−1 X j=1 D X i=j+1 mi− θ D−1 X j=1 D X i=j ni+ D−1 X j=1 4 ˜Fj mjθ exp −θ4 ˜Fj 1 − exp −θ4 ˜Fj − D−1 X j=1 λ = 0 (12) Then multiply equation((10) by θ
−θ D X i=2 i−1 X k=1 mi4 ˜Fk− θ D X i=1 i X k=1 ni4 ˜Fk+ θ D X i=1 mi 4 ˜Fiexp −θ4 ˜Fi 1 − exp−θ4 ˜Fi = 0 and subtract from (12). We obtain
4 ˜FD mDθ exp −θ4 ˜FD 1 − exp−θ4 ˜FD − θnD4 ˜FD+ D−1 X j=1 4 ˜Fjλ = 0.
and equation (11). Accordingly we have
λ = θnD4 ˜FD− 4 ˜FD mDθ exp(−θ4 ˜FD) 1−exp(−θ4 ˜FD) 1 − 4 ˜FD . (13)
• Step 2: Comparing (13) and (8), we can get θnD4 ˜FD− 4 ˜FD mDθ exp(−θ4 ˜FD) 1−exp(−θ4 ˜FD) 1 − 4 ˜FD = mDθ exp(−θ4 ˜FD) 1 − exp(−θ4 ˜FD) − θnD. It follows that mDθ exp(−θ4 ˜FD) 1 − exp(−θ4 ˜FD) − θnD = 0.
• Step 3: Using the estimation of 4 ˜FD in (9), we obtain
mDθ exp(−θ4 ˜FD) 1 − exp(−θ4 ˜FD) − mjθ exp(−θ4 ˜Fj) 1 − exp(−θ4 ˜Fj) + D X i=j+1 mi+ D−1 X i=j ni = 0,
which is used to solve 4 ˜Fj for j = 1, 2, · · · , D − 1.
• Step 4: Finally, obtain the estimate of θ based on equation (10): − D X i=2 i−1 X k=1 mi4 ˜Fk− D X i=1 i X k=1 ni4 ˜Fk+ D X i=1 mi 4 ˜Fiexp −θ4 ˜Fi 1 − exp −θ4 ˜Fi = 0.
The above procedure is summarized below. • Step 1: Obtain θ4 ˜FD from
mDθ exp(−θ4 ˜FD)
1 − exp(−θ4 ˜FD)
− θnD = 0
• Step 2: Solve θ4 ˜Fj for j = 1, 2, · · · , D − 1 based on
mjθ exp(−θ4 ˜Fj) 1 − exp(−θ4 ˜Fj) − θ D X i=j+1 mi− θ D X i=j ni = 0 • Step 3: Since PD
i=14 ˜Fi = 1, then set
PD
i=1θ4 ˜Fi = θ and 4 ˜Fi = θ4 ˜Fi
θ , 1 ≤ i ≤ D − 1
4.2.2 An algorithm based on change of variables
Denote θk = θ4 ˜Fk for k = 1, 2, · · · , D. The log-likelihood function as a function of θk
can be expressed as the following simpler form: log L(θ, ˜F ) = − D X i=2 i−1 X k=1 miθk+ D X i=1 milog [1 − exp (−θi)] − D X i=1 i X k=1 niθk.
Taking the derivatives with respect to θj for D first and then j = 1, . . . , D − 1, we obtain
the following score equations: ∂L(θ, ˜F ) ∂θD = 0 ⇔ mDexp(−θD) 1 − exp(−θD) − nD = 0; and for 1 ≤ j ≤ D − 1 ∂L(θ, ˜F ) ∂θj = 0 ⇔ mjexp(−θj) 1 − exp(−θj) − D X i=j+1 mi− D X i=j ni = 0.
Notice that θj does not depend on other θi for i 6= j.
Now we summarize the numerical algorithm. • Step 1: Solve the equation
mDexp(−θD)
1 − exp(−θD)
− nD = 0
• Step 2: For 1 ≤ j ≤ D − 1, the estimator of θj solves mjexp(−θj) 1 − exp(−θj) − D X i=j+1 mi− D X i=j ni = 0,
• Step 3: Under the constraint, PD
i=14 ˜Fi = 1, we get
θ1+ θ2+ · · · + θD = θ4 ˜F1+ θ4 ˜F2+ · · · + θ4 ˜FD = θ
4.3
Simulation Results
Here we evaluate the non-parametric approach. Two sample sizes with n = 100 and n = 300 are evaluated. The censoring variable C by simulations was generated from a uniform distribution in the the interval [0, k], where k is determined based on Table ??. For θ and p, the average bias and standard deviation are reported based on 500 replications. The result are summarized in Tables A-10 to A-15.
The estimated survival functions using the parametric and non-parametric methods are shown in Figure A-2 to Figure A-10 with n = 100 when the parametric assumption is correctly specified. We see that the parametric estimator is closer to the true function and both parametric and non-parametric estimator are poorer with the increasing of curing rate. In Figure 4-1 below shows the result when the parametric form ˜F (t) is mis-specified. The NPMLE approach is still accurate while the parametric approach fails.
0.1 0.2 0.3 0.4 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE
Chapter 5
Regression Analysis
In practice, data also contain covariate information. Let Z be the covariate so that observed data can be written as {(Xi, δi, Zi)(i = 1, · · · , n)}. The non-mixture cured model
in presence of covariate can be written as
S(t|Z) = exp(−θ(Z) ˜F (t|Z)). (14)
5.1
Classification of Covariate Effects
In (14), covariate Z affects both θ and ˜F (t). Here we consider some special cases which are simplified conditions of the general situation. The long-term effect stands for its influence on the tail exp(−θ(Z)) = p(Z). The short-term effect is often associated with the timing of tumor occurrence for those with N ≥ 1. Hence the model in (14) combines both effects.
Example 1: Short-term effect without long-term effect
Here we assume θ(Z) is independent of Z. For example, if two treatment groups have the same long-term survival rates, and one of them is characterized by more rapidly developing tumors, we could see a significant short-term effect, as shown in Figure 5-1(a). In different groups, the cure rate θ is the same but ˜F (t|Z) depends on Z. Equation (14) can be simplified as
S(t|Z) = exp(−θ ˜F (t|Z)).
Survival curves for the two groups would diverge initially, and converge again as time passes by.
Example 2: Long-term effect without short-term effect
In contrast, we assume that Z affects the cure rate but not ˜F (t), equation (14) can be rewritten as
S(t|Z) = exp(−θ(z) ˜F (t)). Figure 5-1(b) depicts the situation for Z = 0, 1.
0 1 2 3 4 5 0.0 0.2 0.4 0.6 0.8 1.0 t St Short−term effect
(a) short-term effect
0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.0 0.2 0.4 0.6 0.8 1.0 t St Long−term effect (b) Long-term effect
Figure 5-1: The effects of covariate
5.2
Short-term Effect: Proportional Hazards Model
We now illustrate two forms of ˜F (t|Z) which are the most popular choices of regression models in survival analysis. We will also propose model checking procedures to verify the validity of the model assumption. To simplify the presentation, we focus on the two-sample case with Z = 0, 1.
Assume the survival function ˜S(t) follows the form of a PH model such that ˜
S(t|Z = 1) = ˜S(t|Z = 0)k, (15) where k is a pre-determined constant. Notice that
S(t|Z) = exp −θ(z) ˜F (t|Z) = exp θ(z) ˜S(t|Z) − 1 . (16) Thus we have log S(t|Z) θ(z) + 1 = ˜S(t|Z). (17) It follows that log S(t|Z = 1) θ(1) + 1 = ˜S(t|Z = 1) = n ˜S(t|Z = 0)ok = log S(t|Z = 0) θ(0) + 1 k ,
Accordingly log log S(t|Z = 1) θ(1) + 1 = k log log S(t|Z = 0) θ(0) + 1 . (18)
Notice that equation (18) now involves only estimable quantities. The function S(t|z) can be estimated by the Kaplan-Meier estimator
ˆ S(t|Z = z) = Y 0<u≤t 1 − Pn i=1I(Xi = u, δi = 1, Zi = z) Pn i=1I(Xi ≥ u, Zi = z) . Denote ˆθ(z) be the nonparametric estimate of θ(z). Define
Xi = log log ˆS(ti|z = 1) ˆ θ(1) + 1 ! and Yi = log log ˆS(ti|z = 0) ˆ θ(0) + 1 ! ,
where t1 < t2 < · · · < tD are ordered failure points. If the PH assumption holds, equation
(18) indicates that Y follow a linear relationship passing through the origin. We present some plots of ˆS(t|Z) for z = 0, 1 and the corresponding diagnostic plot of Xi versus Yi for
i = 1, 2, · · · , D based on 1000 simulated observations.
The diagnostic plots in Figure 5-2 reveal clear linear pattern for most points. Figure 5-3 presents the plots when the PH assumption is violated. There is a curved relationship between Xi and Yi.
0 1 2 3 4 5 0.0 0.2 0.4 0.6 0.8 1.0 Time Survival
Gamma( 1 , 1 ),cure rate= 0.3 ,censore= 0.4 Gamma( 1 , 2 ),cure rate= 0.5 ,censore= 0.58
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −4 −3 −2 −1 0 −4 −3 −2 −1 0 PH model x y (a) Case I 0 2 4 6 8 10 0.0 0.2 0.4 0.6 0.8 1.0 Time Survival
Weilbull( 1 , 3 ),cure rate= 0.3 ,censore= 0.35 Weilbull( 1 , 2 ),cure rate= 0.5 ,censore= 0.53
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −4 −3 −2 −1 0 −4 −3 −2 −1 0 PH model x y (b) Case II
Figure 5-2: K-M curves and diagnostic plots when the PH assumption holds.
0 5 10 15 20 0.0 0.2 0.4 0.6 0.8 1.0 Time Survival
Log−normal( 0.3 , 5 ),cure rate= 0.3 ,censore= 0.36 Log−normal( 1 , 0.1 ),cure rate= 0.4 ,censore= 0.53
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −5 −4 −3 −2 −1 0 −1.0 −0.8 −0.6 −0.4 −0.2 0.0 PH model x y (a) Case I 0 1 2 3 4 5 0.0 0.2 0.4 0.6 0.8 1.0 Time Survival
Log−normal( 1 , 1 ),cure rate= 0.3 ,censore= 0.62 Log−normal( 2.5 , 0.25 ),cure rate= 0.4 ,censore= 0.79
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −3.0 −2.0 −1.0 0.0 −2.5 −2.0 −1.5 −1.0 −0.5 0.0 PH model x y (b) Case II
5.3
Short Term Effect: Accelerated Failure Time Model
Under the AFT model, the survival function ˆS(t|z) for z = 0, 1 follows the relationship ˜
S(t|Z = 1) = ˜S(kt|Z = 0) for k being a prespecified constant. It follows that
S(t|z = 1) = expnθ(1)h ˜S(t|z = 1) − 1io = expnθ(1)h ˜S(kt|z = 0) − 1io. Accordingly we have log S(t|Z = 1) θ(1) + 1 = ˜S(t|Z = 1) = ˜S(kt|Z = 0) = log S(kt|Z = 0) θ(0) + 1. (19) We want to find a clear relationship from the above equation.
Define p1 < · · · < pM as some constants locating in (0, 1). Then we solve Xi satisfying
log S(xi|Z = 1)
θ(1) = pi, i = 1, 2, · · · , M . Thus, for each pi we have
Xi =SZ=1−1ˆ
n
exppiθ(1)ˆ
o .
Similar steps can be derived based on the right-hand side of equation (19). Set Yi =SZ=0−1ˆ
n
exppiθ(0)ˆ
o
for i = 1, 2, · · · , M . When the AFT model holds, (Xi, Yi) (i = 1, 2, · · · , M ) will follow a
straight line through the origin. Plots following the AFT model are presented in Figure 5-4, in which n=1000.
The diagnostic plots in Figure 5-4 reveal clear linear pattern for most points. Figure 5-5 presents the plots when the AFT assumption is violated. There is a curved relationship between Xi and Yi.
0 1 2 3 4 5 0.0 0.2 0.4 0.6 0.8 1.0 Time Survival
Gamma( 1 , 3 ),cure rate= 0.3 ,censore= 0.38 Gamma( 1 , 3 ),cure rate= 0.5 ,censore= 0.57
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0 1 2 3 0 1 2 3 4 AFT x1 x2 (a) Case I 0 1 2 3 4 5 0.0 0.2 0.4 0.6 0.8 1.0 Time Survival
Gamma( 2 , 3 ),cure rate= 0.3 ,censore= 0.41 Weilbull( 2 , 5 ),cure rate= 0.5 ,censore= 0.64
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0 1 2 3 4 0.0 0.5 1.0 1.5 2.0 2.5 3.0 AFT x1 x2 (b) Case II
Figure 5-4: K-M curves and diagnostic plots when the AFT assumption holds.
0 1 2 3 4 5 0.0 0.2 0.4 0.6 0.8 1.0 Time Survival
Weilbull( 0.5 , 1 ),cure rate= 0.4 ,censore= 0.47 Log−normal( 1 , 5 ),cure rate= 0.3 ,censore= 0.62
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0 1 2 3 4 0 1 2 3 4 AFT x1 x2 (a) Case I 0 2 4 6 8 10 0.0 0.2 0.4 0.6 0.8 1.0 Time Survival
Gamma( 5 , 2 ),cure rate= 0.4 ,censore= 0.47 Log−normal( 1 , 5 ),cure rate= 0.3 ,censore= 0.48
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0 2 4 6 0.0 0.5 1.0 1.5 2.0 2.5 3.0 AFT x1 x2 (b) Case II
Chapter 6
Conclusion
In the thesis, we study the non-mixture approach for analyzing survival data in presence of cure. This formulation has an interesting biologial interpretation. In parametric analysis, we find that outliers in T which often occur when N = 1 will affect estimation of the cure rate. For nonparametric inference, we propose two algorithms to solve the score funcions of nonparametric MLE. One is the classical Lagrange multiplier method and the other is by change of variables. Two regression models are considered under the simplified two-sample setting. One is the proportional hazard model and the other is the accelerated failure time model. We propose diagnostic plots which can verify the form of regression effect.
References
[1] J. Berkson and R. P. Gage (1952). Survival Curve for Cancer Patients Following Treatment. J. Amer. Statist. Assoc. 47, 501-515.
[2] Boag, JW. (1949). Maximum Likelihood Estimates of the Proportion of Patients Cured by Cancer Therapy. J. Royal Statist. Soc., Series B 11:15-44.
[3] Tsodikov, A. D. (1998). A Proportional Hazards Model Taking Account of Long-Term Survivors. Biometrics 54, 1508-1516.
[4] Tsodikov, A. D. (2001). Estimating of Survival Based on Proportional Hazards When Cure is a Possibility. Mathematical and Computer modelling. 33, 1227-1236.
[5] Tsodikov, A. D., Ibrahim, J. G. and Yakovlev A. Y. (2003). Estimating Cure Rates From Survival Data : An Alternative to Two-Component Mixture Models . E. J. Am. Statist. Assoc. 98, 1063-1078.
Appendix: Additional Figures
0 5 10 15 20 0.0 0.4 0.8 survival function t ST |Z (( t )) Gamma( 1 , 1 ) Gamma( 4 , 1.5 ) Gamma( 6 , 2 ) 0 5 10 15 20 0.0 0.1 0.2 0.3 0.4 density function t fT|Z (( t )) Gamma( 1 , 1 ) Gamma( 4 , 1.5 ) Gamma( 6 , 2 ) 0 5 10 15 20 0.0 0.1 0.2 0.3 0.4 hazard function t hT |Z (( t )) Gamma( 1 , 1 ) Gamma( 4 , 1.5 ) Gamma( 6 , 2 ) 0 5 10 15 20 0.0 0.4 0.8 survival function t ST |Z (( t )) Weibull( 1 , 3 ) Weibull( 2 , 10 ) Weibull( 4 , 15 ) 0 5 10 15 20 0.0 0.1 0.2 0.3 0.4 density function t fT|Z (( t )) Weibull( 1 , 3 ) Weibull( 2 , 10 ) Weibull( 4 , 15 ) 0 5 10 15 20 0.0 0.1 0.2 0.3 0.4 hazard function t hT |Z (( t )) Weibull( 1 , 3 ) Weibull( 2 , 10 ) Weibull( 4 , 15 ) 0 5 10 15 20 0.0 0.4 0.8 survival function t ST |Z (( t )) Log−normal( 1 , 1 ) Log−normal( 2 , 0.4 ) Log−normal( 2.5 , 0.25 ) 0 5 10 15 20 0.0 0.1 0.2 0.3 0.4 density function t fT|Z (( t )) Log−normal( 1 , 1 ) Log−normal( 2 , 0.4 ) Log−normal( 2.5 , 0.25 ) 0 5 10 15 20 0.0 0.1 0.2 0.3 0.4 hazard function t hT |Z (( t )) Log−normal( 1 , 1 ) Log−normal( 2 , 0.4 ) Log−normal( 2.5 , 0.25 )0 1 2 3 4 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.33 0.0 0.5 1.0 1.5 2.0 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.45 0.0 0.5 1.0 1.5 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.49
1 2 3 4 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.34 1 2 3 4 5 6 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.48 1 2 3 4 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.64
1 2 3 4 5 6 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.3 1 2 3 4 5 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.42 1 2 3 4 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.51
0 2 4 6 8 10 12 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.33 0 1 2 3 4 5 6 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.38 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.53
5 10 15 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.36 5 10 15 20 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.47 5 10 15 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.6
4 6 8 10 12 14 16 18 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.33 5 10 15 20 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.47 6 8 10 12 14 16 18 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.63
0 2 4 6 8 10 12 14 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.34 0 1 2 3 4 5 6 7 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.43 0 1 2 3 4 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.51
2 4 6 8 10 12 14 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.3 2 4 6 8 10 12 14 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.38 4 6 8 10 12 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.5
10 15 20 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.39 8 10 12 14 16 18 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.49 8 10 12 14 16 18 0.0 0 .2 0.4 0 .6 0.8 1 .0 Time S(t) True MLE NPMLE p= 0.3 , censoring rate = 0.57
Figure A-10: Estimated survival functions when F˜ is correctly specified as log-normal(2.5,0.25)
