Analytical Strategies for Failure Time Data with a Cured Fraction

Some failure time data come from a population that consists of some subjects who are susceptible to and others who are non-susceptible to the event of interest. The data typically have heavy censoring at the end of the follow-up period, and a traditional survival analysis would not always be appropriate, yet it is commonly seen in literatures. For such kind of data, we carry out simulation studies to compare the performances of the Cox’s PH model with the proportional hazards mixture cure (PHMC) model and the accelerated failure model (AFT model) with the AFT mixture cure (AFTMC) model respectively. Then we apply the models to the datasets of Lung Cancer and Eastern Cooperative Oncology Group (ECOG) phase III clinical trial E1684. The conclusions are as follows. The PHMC model and the AFTMC model do not have obvious advantages for time-to-event data without a cured fraction. In this case, it is recommended to use the Cox’s PH model or AFT model for analysis. If some subjects are non-susceptible to the event of interest in the data, it is recommended to use the PHMC model or AFTMC model for analysis, however, which may need a sufficient sample size.


BACKGROUND
The Cox's PH model and AFT model are two most commonly used models for analyzing survival data [1,2]. The underlying assumption of these two models is that each patient will eventually experience the event of interest (death, relapse or recurrence of cancer) if the follow-up time is sufficiently long [3,4,5,6,7]. However, with the development of medical science, advanced treatments can substantially improve survival rates for many diseases that are previously poor prognosis and even fatal, which means that a substantial proportion of subjects will never experience the event of interest [4,5,6]. An estimated Kaplan-Meier survival curve that presents a steady and long plateau at the tail and levels off at a value greater than zero with heavy censoring can be used as empirical evidence for the existence of cured individuals [3,5,8]. In fact, for most diseases, we are more concerned with the improvement of the cured rate rather than prolongation of time to the event of interest. For such time-to-event data, the cured rate can't be estimated with traditional methods such as Cox's PH model or AFT model, even some inappropriate conclusions may be drawn. 3 The cure model, which has been proposed by Boag (1949) and Berkson (1952) [9,10] can be used to analyze such data. The mixture cure model assumes that the study population is a mixture of susceptible subjects, that are uncured or will experience the event of interest, and non-susceptible subjects, that are cured or will never experience the event of interest [11,12]. At present, parametric approaches have been proposed for the mixture cure model by many scholars [12,13,14]. Since it is usually difficult to testify a parametric assumption, there has been increasing enthusiasm in the semiparametric mixture cure model [3,5,7,8,15]. An attractive characteristic of the mixture model is that it includes two parts which can be interpreted separately, in particular, the estimation of the proportion of the cured patients and of the failure time distribution of uncured patients.
However, for the data with a cured fraction, a great deal of literatures indicated that non-statistical professional researchers would still choose the traditional Cox's PH model or AFT model for analysis. In view of this situation, this study intended to compare the performances of the Cox's PH model with the proportional hazard mixture cure model (PHMC) and the AFT model with the accelerated failure time mixture cure model (AFTMC) respectively. In the meantime, two real data were analyzed to illustrate that for survival data with a cured fraction, the mixture cure model is a more effective method compared with traditional survival analysis methods.
The rest of the paper is organized as follows. In section 2, the mixture cure model and its estimation method are described. In section 3, we carry out simulation studies to compare the performances of the Cox's PH model with the PHMC model and the AFT model with the AFTMC model. In section 4, simulation results are illustrated using Lung Cancer Data and Eastern Cooperative Oncology Group (ECOG) Data. In section 5, we will give some discussions.

Let
be the indicator variable, where an individual will eventually ( = 1) or never ( = 0) experience the event of interest. When = 1, let denote the time to the occurrence of an event of interest and ( | = 1) be the survival function of . The mixture cure model can be denoted as: Note that ( ) → 1 − if → ∞. = ( = 1) represents the probability of a subject being 4 uncured and is also referred to as "incidence". ( | = 1) is referred to as "latency". Assume that a logistic regression is employed to estimate the cured rate 1 − ( ) Where 0 is a constant, = ( 1 , … ) is a vector of unknown parameters and = ( 1 , … , ) is a vector of covariates. If Cox's PH model is considered as "latency", which is the mixture cure model is called the proportional hazard mixture cure model (PHMC model). If AFT model is considered as "latency", which is the mixture cure model is called the accelerated failure time mixture cure model (AFTMC model).
The Semiparametric Estimation Method of PHMC model and AFTMC model can be seen the literatures [3,5,15,16]. R package "smcure" was used to fit the PHMC model and AFTMC model.

The Cox's PH model and the PHMC model
3.1.1 Simulation Settings. We generate the data according to the structure of the PHMC model.
The data should satisfy proportional hazards assumption, and the methods of tests and diagnostics are based on weighted residuals [17]. The incidence part of the PHMC model is assumed that conditional probability of being cured follows a logistic model: where 1~( 0,1) is sampled from a standard normal distribution and x 2~( 0.5) follow the binomial distribution. We assume 1 = 2 = l (2) and 0 is selected such that the conditional probability of being cured are: 0, 0.1, 0.2, 0.3, 0.4 and 0.5. In the latency part of the PHMC model, we generate 1~( 0,1) that follow a standard normal distribution and z 2~( 0.5) that obey the binomial distribution. Simultaneously, we consider 1 = 2 = l (2) corresponding to 1 and 2 respectively. According to the method proposed by Bender [18], we generate the survival time. We assume the baseline hazard function is from the Weibull distribution and both the scale and shape parameters are set to be 1. The censored time is from the exponential distribution, the parameter λ controls the censored rate, and the value of λ is 5 adjusted so that the difference between the average censored rate and cured rate for each set of data are: 0.1, 0.2 and 0.3. The sample size is set to be100 and 200, and the number of simulations is 1000.
The evaluation criteria include Estimated biases, Mean square error (MSE), Confidence interval capture rate, and Concordance probability K index [19,20].  The MSE of the two models are shown in Figure 2. When the cured rate is 0, the MSE is similar in two models. With the sample size being100, when both the cured rate and the censored rate are high (e.g. cured rate=0.5 and censored rate=0.8), the MSE of the PHMC model will become remarkably large, especially for binary covariate. In this case, the PHMC model loses its accuracy. The K indexes of the two models are presented in Figure 4. When the cured rate is 0, the K indexes of the two models are almost equal. When the cured rate ranges from 0.1 to 0.5, the K indexes of the PHMC model are always greater than the Cox's PH model, indicating that when survival data has a substantial proportion of subjects being cured, the overall fit effect of the PHMC model is better. When the corresponding cured rate (0.1 to 0.5) is fixed, the K indexes of the PHMC model are greater than that of the Cox's PH model. With the increase of censored rate, the K indexes of both the PHMC model and the Cox's PH model will increase, implying that the smaller the cured proportion in censored individuals is, the better the model fit effect is both in the two models.    The MSE of the two models are exhibited in Figure 6. When the cured rate is 0, the MSE of the two models are near 0. When the cured rate is from 0. The K indexes of the two models are shown in Figure 8. When the cured rate is 0, the K indexes

Lung Cancer Data
Survival in patients with lung cancer is from the North Central Cancer Treatment Group [24]. This study was aimed to explore factors that may affect patients' prognosis. The data contained 228 patients. Covariates included Survival time, censored status, age, sex, ECOG score, Karnofsky performance score rated by doctor, Karnofsky performance score rated by patient, calorie consumed in meals, and weight loss in past six months. K-M survival curves of male and female patients are presented in Figure 5, which shows that after approximately 1000 days of follow-up, the survival rate is almost 0, implying that there is no cured individual. With LASSO regression [25,26], the covariates selected were sex, ECOG performance score, and 14 Karnofsky performance score rated by patients. The proportional hazards test exhibits the data meets the proportional hazards assumption (P=0.165). We used Cox's PH model and PHMC model to analyze it respectively. The results are exhibited in Table 1. The results of the incidence part of the PHMC model show that three covariates are not statistically significant (P values were greater than 0.9) and aren't displayed in Table1. As shown from (0=control group, 1=treatment group). The K-M survival curve of the treatment group and the control group is presented in Figure 6, which show a steady and long plateau at the tail, implying that there may be cured individuals in melanoma patients. Proportional hazards test shows the proportional hazards assumption is met (P=0.404), then Cox's PH model and PHMC model were used for analysis. The results are presented in Table 2 and Table 3. In Table 2, the results of the incidence part in PHMC model exhibit that TRT (OR=0.555, P=0.053) may be a factor that affects the cured rate of melanoma patients, suggesting that high dose interferon alpha-2b (IFN) may increase the cured rate of melanoma patients.
Both models in Table 3  treatments. The methods of traditional survival analysis are more sensitive to find out treatments that delay the occurrence of events (higher power) but are less sensitive to treatments that increase the cured rate, especially when the length of follow-up is limited, giving doctors a false estimate of the best treatment. For those studies with long follow-up and the existence of cured individuals, the power of the traditional survival analysis methods will reduce. Therefore, unless all patients will eventually experience the event of interest or there are no long-term survivors, the primary measure to evaluate patient's benefit should be the patient's cured rate rather than the risk ratio. Compared with traditional Cox or AFT models, which only focus on the time of event and analysis results may mislead doctors or patients, the cure model can provide more specific and accurate information on the long-term survival benefit of the treatment.
In this article, the difference between the average censored rate and the average cured rate is set from 0.1 to 0.3, and scenarios of greater censored rate are not involved in this study. Additionally, appropriate sample size for cure model has not yet been studied. Furthermore, we just focus on mixture cure models in this article and non-mixture cure models are not involved in our comparison.
A competing risks derivation of a mixture cure model has been proposed by JB. GreenHouse as early as in 1984 [28]. Semiparametric accelerated failure time cure rate mixture models with competing risks has been put forward recently [29]. Thus it is worth further research to compare the competing risks model with its mixture model.