Improvement Screening for Ultra-High Dimensional Data with Censored Survival Outcomes and Varying Coefficients

1 Introduction

Survival analysis involves a collection of commonly used statistical methods for the analysis of failure time data such as biological death, mechanical failure, or credit default. In this context, death or failure is also referred to as an “event”. The time-to-event data is usually censored due to the termination of study. One practical objective of survival analysis is to identify the risk factors and quantify their effects on the failure outcome through hazard regression. Accurate risk prediction can be quite useful for patient treatment selection, prevention strategy or disease management in evidence-based medicine. Specifically, we may study how the conditional hazard function of survival time T depends on the observed p-dimensional covariate X = x

According to the probabilistic definition, the conditional hazard function may be interpreted as the instantaneous failure rate at time t given a particular covariate profile x. The most popular regression model is Cox proportional hazards (PH) model introduced by Cox [1, 2]

in which h0(t) is the baseline hazard. A linear model formulation ψ(x)=xTβ is usually adopted. See Klein and Moeschberger [3] and references therein for more detailed literature on Cox PH model.

In recent biomedical studies, especially in cancer research, profiling approaches have been extensively conducted, measuring genome-wide mRNA gene expression levels, DNA modifications (e.g., copy number variation – CNV, SNPs), epigenetic regulation (e.g., DNA methylation and histone modifications), post-transcriptional regulations (e.g., microRNA expression) and others. These studies naturally produce massive data and require special statistical methodology. A representative example is The Cancer Genome Atlas (TCGA) in the USA, which provides comprehensive genomic characterization for a cohort of cancer and normal samples. High-dimensional data offer a unique opportunity to more comprehensively describe the ethology and prognosis of important diseases. Correspondingly in statistical literature we have seen a surge of interest in ultra-high dimensionality survival data set where the covariate dimensionality p grows exponentially or non-polynomially fast with sample size n, i.e., log(p)=O(nα) for some α∈(0,1/2). For ultra-high dimensional data, the first step of data analysis is to rank the importance of covariates based on their marginal associations with the outcome variable and screen out unimportant covariates from the ordered list. Only after this necessary screening stage one can proceed to construct a parsimonious regression model applying more sophisticated model building techniques.

For variable screening procedure, Fan and Lv [4] proposed sure independence screening (SIS) based on marginal correlation ranking and also extended SIS to iterative version (ISIS) to handle difficult cases such as when some important predictors are marginally uncorrelated with the response. In order to deal with more complex real data, Fan et al. [5] extended SIS and ISIS to generalized linear models and robust regression; Fan and Song [6] and Fan et al. [7] extended SIS to nonparametric additive models. Moreover, variable screening methods have been quickly adapted in survival analysis with focus on Cox PH model [8×10], and [11] or other regression models [12, 13], and [14]. Recently, He et al. [15], Song et al. [16], and Li et al. [17] proposed non-model based screening methods.

Most statistical methods developed for failure time data assume that covariate effects on the logarithm of the hazard function are linear and the regression coefficients are constant parameters Fleming and Harrington [18] and Andersen et al. [19]. However, true covariate effects can be more complex than the log-linear effect. An important extension of the standard regression model with constant coefficients is the varying-coefficient model. In recent decades, the varying functional effect in survival analysis has been carefully studied by Murphy [20], Cai and Sun [21], Tian et al. [22], Cai et al. [23], and Cai et al. [24], among others. We will consider varying-coefficient Cox model in this paper and aim at providing a novel screening procedure to improve the model performance.

To assess the gain in model development, researchers proposed many improvement statistics in the literature. For example, Pencina et al. [25] proposed the integrated discrimination improvement (IDI) and the net reclassification improvement (NRI), which have drawn much attention in medical research. The IDI is the improvement of the integral of sensitivity and specificity over all possible cut-off values. Extensions of the IDI to survival data have been proposed by Uno et al. [26], Kerr et al. [27] and Chi and Zhou [28]. NRI attempts to quantify how well a new model correctly reclassifies subjects based on comparison between a baseline model and a new model with the added markers. See Pencina et al. [29], Uno et al. [30] and Shi et al. [31] for more discussion. We use such novel improvement statistics in this paper to assess the gain in screening models as well as some standard improvement statistics. We will then use the improvement measures as an effective way to screen ultra-high dimensional biomarkers.

The main contribution of this paper can be summarized as follows. A formal improvement screening methodology is proposed for survival analysis. In particular, the hazards effects of new biomarkers are modelled as varying coefficient functions of possible confounders in the Cox PH model. To the best of our knowledge, for the present setup, improvement screening has not been thoroughly discussed for varying-coefficient Cox PH model yet. Moreover, motivated by the growing importance of biomarker screening problem, in this paper we consider the ultra-high dimensional data setting, which presents a tremendous novel challenge, and, as far as we have reviewed, there is no previous work on this topic. The rest of the paper is organized as follows. In Section 2, we formulate the varying-coefficient Cox PH model and propose a detailed procedure for our improvement screening methodology. In Section 3, we evaluate the proposed screening procedures through simulation studies and illustrate our methods via application to the lung cancer data set in Section 4. Some discussion is provided in Section 5.

2 Methodology

Many survival analysis data set involves massive number of covariates, such as genetic and genomics studies. Usually low-dimensional confounders are also collected, such as demographical variables, known risk factors, environmental exposures or other variables [32]. Including relevant confounding variables in the model may produce better prediction results. Furthermore, comparison of the results across different studies is easier since most researches on the same response would adjust a similar set of confounders. In our paper, we consider the baseline model with only a confounding variable and add biomarkers to improve the model performance. We present several novel metrics to quantify the improvement from the new markers for prediction of the patients risk and target to screen ultra-high dimensional markers in varying-coefficient models.

3 Simulation

To examine the performance of the above improvement screening methods, we conduct simulation studies under practical settings. In all simulations, we consider the following true model:

where Xi=(Xi1,…,Xip),i=1,…,n are i.i.d. from multivariate normal distribution with mean 0p, and Cov(Xij,Xik)=ρ|j−k|. Both censoring time C and confounder U are generated from uniform distribution on (0,10). The correlation ρ is set to be 0.8. For the IDI and NRI approaches, we choose the first and third quartiles (corresponding to α=0.25 and α=0.50, respectively) of the observed time. We consider the ultra-high dimensional case with p=1000 and total sample size n=400 and 800. In our implementation, we adopt cubic spline with one inner knot without intercept. In the first two cases, we consider p∗=20 important variables among 1000 variables; in the latter four cases we increase the number of important variables to 30 and 50. Our purpose is to screen out the important markers from a large set of candidates. Case I Let g0(u)=1gk(u)=(3u−1)2 for k=1,5,9,13,17,gk(u)=3u3+0.5u for k=2,6,10,14,18,gk(u)=3sin(2πu) for k=3,7,11,15,19 and gk(u)=3exp(−4u) for k=4,8,12,16,20.Case II Let g0(u)=1,gk(u)=(3u−1)2 for k=1,5,9,13,17,gk(u)=2exp(−(3u−1)2)+exp(−4(u−3)2) for k=2,6,10,14,18,gk(u)=0.2sin(2πu)+0.2cos(2πu)+0.3sin(2πu)2+0.4cos(2πu)3+0.5sin(2πu)3 for k=3,7,11,15,19 and gk(u)=3exp(−4u) for k=4,8,12,16,20.Case III Let g0(u)=1,gk(u)=1.5u for k=1,7,13,19,25,gk(u)=2u for k=2,8,14,20,26,gk(u)=3exp(−4u) for k=3,9,15,21,27,gk(u)=2u3 for k=4,10,16,22,28,gk(u)=2exp(−(3u−1)2)+exp(−4(u−3)2) for k=5,11,17,23,29,gk(u)=(3u−1)2 for k=6,12,18,24,30.Case IV Let g0(u)=1,gk(u)=(3u−1)2 for k=1,⋯,5,gk(u)=2exp(−(3u−1)2)+exp(−4(u−3)2) for k=6,⋯,10,gk(u)=2u3 for k=11,⋯,15 and gk(u)=3exp(−4u) for k=16,⋯,20,gk(u)=1.5u3+0.5u for k=21,⋯,25 and gk(u)=2.5u for k=26,⋯,30.Case V Let g0(u)=1,gk(u)=1.5u for k=1,6,11,16,21,26,31,36,41,46,gk(u)=2u for k=2,7,12,17,22,27,32,37,42,47,gk(u)=3exp(−4u) for k=3,8,13,18,23,28,33,38,43,48,gk(u)=2u3 for k=4,9,14,19,24,29,34,39,44,49,gk(u)=2exp(−(3u−1)2)+exp(−4(u−3)2) for k=5,10,15,20,25,30,35,40,45,50.Case VI Let g0(u)=1,gk(u)=2u for k=1,6,11,16,21,26,31,36,41,46,gk(u)=3exp(−4u) for k=2,7,12,17,22,27,32,37,42,47,gk(u)=3u3+0.5u for k=3,8,13,18,23,28,33,38,43,48,gk(u)=2exp(−(3u−1)2)+exp(−4(u−3)2) for k=4,9,14,19,24,29,34,39,44,49,gk(u)=(3u−1)2 for k=5,10,15,20,25,30,35,40,45,50.

The selected functions are shown in Figure 1. We design the settings where all of the non-zero functions used to generate data are roughly equally related to the response. Similar settings were adopted by many earlier authors for numerical studies. They represent the practical study where only sparse sets of variables contribute to the response. We have examined other settings with different p values in extensive simulation studies and obtain quite similar findings as we report here.

Figure 1

The plots of functions gk(u) defined in simulation for each of the four cases.

The simulation results are based on 500 independent runs. Besides the four improvement screening methods, we also report the screening results by using the P-value from fitting marginal Cox regression model to each marker. P-value is usually obtained from a Wald test and is commonly used in practice to evaluate the significance of variables. This marginal screening method is widely applied in high-dimensional research.

The first criterion to evaluate the screening performance is the minimum model size which is the smallest number of covariates that we need to include in a model in order to ensure that all the active variables are selected. In Figure 2, we display the box plots of the minimum model size from the simulations. In most cases, the minimum model sizes concentrate around the number of truly important covariates, indicating satisfactory screening order.

Figure 2

The box plots of the minimum model size covers all important covariates.

In our implementation, as recommended by Fan and Lv [4], Xia et al. [12] and Cheng et al. [35], we keep the top ⌊nlogn⌋ variables in the ranked list and screen out the rest of the covariates. In Table 1, we report the coverage probability that this selected set of top ⌊nlogn⌋ covariates includes all truly important ones. We observe from Table 1 that all screeners perform quite well with the coverage probability close to one. As the number of important covariates increase in Cases III, IV, V and VI, it is slightly harder to select all of them by only keeping the top ⌊nlogn⌋ variables in the ranked list. In general performances of all methods improve as the sample size increases.

We note that there was little work on improvement screening available in the literature. Though C-index and l2-norm are commonly used in practice, screening based on their improvement has not been fully investigated. It is quite appealing to apply improvement screening based on these familiar statistics. We note that IDI and NRI focus on different aspects of accuracy improvement in diagnostic medicine and their results thus should be interpreted differently.

4 Example: lung cancer

Lung cancer is the leading cancer killer in both men and women in the US. About 1 out of 4 cancer deaths are from lung cancer. Each year, more people die of lung cancer than of colon, breast, and prostate cancers combined [41]. Assessing gene expressions in the lung cancer and detection of the relevant gene expressions are thus crucial for lung cancer prevention and treatment. A substantial number of studies have reported the development of gene expression-based prognostic signatures for lung cancer. The ultimate aim of such studies should be the development of well-validated clinically useful prognostic signatures that improve therapeutic decision making beyond current practice standards [42].

Our example was extracted from a large, training-testing, multi-site blinded validation study [43] based on 442 lung adenocarcinomas. The gene expression data in this example were generated by four different laboratories under a common protocol. The same data set has been analyzed by Li et al. [17] and Xie et al. [44]. In this analysis, we adopt age as the potential confounder, as it is strongly associated with end points of interest in occupational epidemiology [45]. After removing cases with missing information, a total of 439 subjects were included in downstream analysis. The median follow-up time is 46 months (range: 0.03 to 204 months) and the median age at diagnosis is 65 years (range: 33 to 87 years). The overall censoring rate is 46.47%. For each subject, the expressions of 22,283 genes are available. Our target is to identify a set of gene expression signatures that are highly predictive in the lung cancer.

We first fit a baseline Cox PH model including only the age covariate whose effect is modelled as an unknown function and then evaluate new models by adding the genes with varying coefficients as described in Section 2. We then quantify the incremental value for each gene by the four proposed improvement screening approaches. For illustration, we also discretized the event time using the first (α = 0.25), second (α = 0.50) and the third quartile (α = 0.75) points in IDI and NRI methods for the short-term, mid-term and long-term prediction, respectively.

In Table 2, we compare the screening results for the top 20 genes obtained by the proposed improvement screening methods and also the P-value based screening approach. We display the IDs of selected genes and the corresponding screening measures along with their 95% confidence intervals calculated by 500 perturbation-resampling. We highlight genes which are simultaneously selected by at least three screeners. Despite certain degree of overlapping, we note that these methods generally produce distinct list of top genes. Depending on the goal and emphasis of the study we may consider choosing appropriate screening statistics. For example, gene 750 achieves the highest IDI at α=0.75 and the third highest NRI value at α=0.75, yielding 14.3% and 35.2% increment from the baseline model, respectively. If the aim is to improve the prediction for long-term survival, we may consider this gene as an important biomarker. The top genes selected by C-statistics give high concordance for risk prediction. For example, gene 7,603 gives an increase of 15.5% in Harrell’s C, suggesting that more than 15.5% cases predicted by the new marker would yield risk prediction in the same direction as the failure time. The top genes selected by the l2-norm imply a greater contribution towards the hazard rates and all top genes selected by p-value are highly significant.

We also compare our proposed improvement screening approaches to improvement screening without the varying coefficients. In particular, we consider fitting the following two nested models:

where M0 is the baseline model only with a confounder variable U in a parametric form and Mk⋆ is the kth new model with the additional marker Xk on top of the baseline. The same improvement screening statistics are implemented and the results are shown in Table 3. In general IDI, NRI and the improvement in Harrell’s C values in Table 3 are smaller than those in Table 2. For example, at α=0.25, the first gene selected under IDI using the varying-coefficient attains an IDI 0.12 while that using the constant coefficient achieves the highest IDI 0.10. Such a comparison may suggest that the varying-coefficient model could lead to relatively stronger improvement for the hazards regression than the constant coefficient model. Using more sophisticated semiparametric varying-coefficient model may provide more accurate model specification for a high-dimensional analysis.

We next report the Pearson’s correlation matrix for these measurements in Table 4. The correlation indicates the degree of pairwise agreement between two quantitative improvement statistics. For example, the correlation between NRI and IDI is computed using the following sample statistic:

where NRIi and IDIi are the estimated NRI and IDI for the ith marker and NRˉI and IDˉI are the sample means of NRI and IDI, respectively. We observe that the degree of linear correlation between Harrell’s C-Statistics & p-value is relatively high with Pearson’s correlation 0.619. This may suggest a moderate agreement between these two methods. However, the degree of linear correlation between l2-norm and other measures are relatively low. Furthermore, there is a relative high agreement between IDI & NRI (or Harrell’s C-Statistics) and between NRI & Harrell’s C-Statistics for α=0.5 and α=0.75 but relatively low agreement for α=0.25. Since these methods may address different types of improvement, we thus expect such an empirical disagreement.

We display the estimated varying-coefficient functions gk(age) for the top three genes in the ranked list by methods IDI(α = 0.25), IDI(α = 0.75), NRI(α = 0.25) and NRI(α = 0.75) in Figure 3. All of the functions are quite different from a straight line. Source code for this data analysis is available at http://www.stat.nus.edu.sg/ stalj/.

Figure 3

Estimated function of gk(age) for the top three gene expressions in the ranked list by IDI(α = 0.25), IDI(α = 0.50), IDI(α = 0.75), NRI(α = 0.25), NRI(α = 0.50) and NRI(α = 0.75) using cubic spline with one inner knot and without intercept.

5 Discussion

We intend to provide a practical screening approach when the goal is to seek markers with the greatest improvement over existing risk prediction models. Besides the four approaches we examine in this paper, there are additional improvement screening procedures available such as the improvement in odds ratios, hazards ratios, P-values and other parametric and nonparametric measures. They may also be implemented in a similar way. More theoretical development is needed to justify the theoretical properties such as the sure screening properties and the control of false discovery rates.

We also want to highlight the importance of considering sophisticated nonparametric modelling techniques such as the spline estimates of varying coefficients. Many biological effects are nonlinear and a simple description using parametric methods may be misleading and non-informative. Nonparametric and semiparametric regression methods are well developed in theory and implementation and could serve as useful alternative choices for practitioners.

The method for varying coefficient function is limited to a scalar confounder U. Many applications may involve a set of confounders such as age, gender and other low-dimensional variables. It is usually difficult to smooth more than one index variables in the varying coefficient due to the curse of dimensionality. To address such a concern one may need to adopt certain functional structure for the multiple confounder effects. Model selection and model validation are necessary development for such a goal.