High-dimensional variable selection for Cox’s proportional hazards model

First let us review the definition of sure screening property.

For each covariate $X_{m}$ ($1\leq m\leq p$), define its marginal utility as the maximum of the partial likelihood of the single covariate:

Here $x_{im}$ is the $m$-th element of $\mbox{\bf x}_{i}$, namely $\mbox{\bf x}_{i}=(x_{i1},x_{i2},\cdots,x_{ip})^{T}$. Intuitively speaking, the larger the marginal utility is, the more information the corresponding covariate contains the information about the survival outcome. Once we have obtained all marginal utilities $u_{m}$ for $m=1,2,\cdots,p$, we rank all covariates according to their corresponding marginal utilities from the largest to the smallest and select the $d$ top ranked covariates. Denote by $\mathcal{I}$ the index set of these $d$ covariates that have been selected.

The index set $\mathcal{I}$ is expected to cover the true index set $\mathcal{M}^{*}$ with a high probability, especially when we use a relative large $d$. This is formally shown by Fan and Lv (2008) for the linear model with Gaussian noise and Gaussian covariates and significantly expanded by Fan and Song (2010) to generalized linear models with non-Gaussian covariates. The parameter $d$ is usually chosen large enough to ensure the sure screening property. However the estimated index set $\mathcal{I}$ may also include a lot of unimportant covariates. To improve performance, the penalization based variable selection approach can be applied to the selected subset of the variables $\{X_{j},j\in\mathcal{I}\}$ to further delete unimportant variables. Mathematically, we then solve the following penalized partial likelihood problem:

where $\mbox{\bf x}_{\mathcal{I},i}$ denotes a sub-vector of $\mbox{\bf x}_{i}$ with indices in $\mathcal{I}$ and similarly for $\mbox{\boldmath$\beta$}_{\mathcal{I}}$. It will lead to sparse regression parameter estimate $\hat{}\mbox{\boldmath$\beta$}_{\mathcal{I}}$. Denote the index set of nonzero components of $\hat{}\mbox{\boldmath$\beta$}_{\mathcal{I}}$ by $\hat{\mathcal{M}}$, which will serve as our final estimate of $\mathcal{M}^{*}$.

Fan and Lv (2008) pointed out that SIS can fail badly for some challenging scenarios such as the case that there exist jointly related but marginally unrelated covariates or jointly uncorrelated covariates having higher marginal correlation with the response than some important predictors. To deal with such difficult scenarios, iterative SIS (ISIS) has been proposed. Comparing to SIS which is based on marginal information only, ISIS tries to make more use of joint covariates’ information.

The iterative SIS begins with using SIS to select an index set $\hat{\mathcal{I}}_{1}$, upon which a penalization based variable selection step is applied to get regression parameter estimate $\hat{}\mbox{\boldmath$\beta$}_{\hat{\mathcal{I}}_{1}}$. A refined estimate of the true index set is obtained and denoted by $\hat{\mathcal{M}}_{1}$, the index set corresponding to nonzero elements of $\hat{}\mbox{\boldmath$\beta$}_{\hat{\mathcal{I}}_{1}}$.

As in Fan, Samworth and Wu (2009), we next define the conditional utility of each covariate $m$ that is not in $\hat{\mathcal{M}}_{1}$ as follows:

This conditional utility measures the additional contribution of the $m$th covariate given that all covariates with indices in $\hat{\mathcal{M}}_{1}$ have been included in the model.

Once the conditional utilities have been defined for each covariate that is not in $\hat{\mathcal{M}}_{1}$, we rank them from the largest to the smallest and select these covariates with top rankings. Denote the index set of these selected covariates by $\hat{\mathcal{I}}_{2}$. With $\hat{\mathcal{I}}_{2}$ having been identified, we minimize

with respect to $\mbox{\boldmath$\beta$}_{\hat{\mathcal{M}}_{1}\cup\hat{\mathcal{I}}_{2}}$ to get sparse estimate $\hat{}\mbox{\boldmath$\beta$}_{\hat{\mathcal{M}}_{1}\cup\hat{\mathcal{I}}_{2}}$. Denote the index set corresponding to nonzero components of $\hat{}\mbox{\boldmath$\beta$}_{\hat{\mathcal{M}}_{1}\cup\hat{\mathcal{I}}_{2}}$ to be $\hat{\mathcal{M}}_{2}$, which is our updated estimate of the true index set $\mathcal{M}^{*}$. Note that this step can delete some variables $\{X_{j}\in\hat{\mathcal{M}}_{1}\}$ that were previously selected. This idea was proposed in Fan et al. (2009) and is an improvement of the idea in Fan and Lv (2008).

The above iteration can be repeated until some convergence criterion is reached. We adopt the criterion of either having identified $d$ covariates or $\hat{\mathcal{M}}_{j}=\hat{\mathcal{M}}_{j-1}$ for some $j$.

Fan, Samworth and Wu (2009) noted that the idea of sample spliting can also be used to reduce the false selection rate. Without loss of generality, we assume that the sample size $n$ is even. We randomly split the sample into two halves. Then apply SIS or ISIS separately to the data in each partition to obtain two estimates $\hat{\mathcal{I}}^{(1)}$ and $\hat{\mathcal{I}}^{(2)}$ of the true index set $\mathcal{M}^{*}$. Both these two estimates could have high FSRs because they are based on a simple and crude screening method. Yet each of them should include all important covariates with high probabilities. Namely, important covariates should appear in both sets with probability tending to one asymptotically. Define a new estimate by intersection $\hat{\mathcal{I}}=\hat{\mathcal{I}}^{(1)}\cap\hat{\mathcal{I}}^{(2)}$. The new estimate $\hat{\mathcal{I}}$ should include all important covaraites with high probability as well due to properties of each individual estimate. However by construction, the number of unimportant covariates in the new estimate $\hat{\mathcal{I}}$ is much smaller. The reason is that, in order for an unimportant covariate to appear in $\hat{\mathcal{I}}$, it has to be included in both $\hat{\mathcal{I}}^{(1)}$ and $\hat{\mathcal{I}}^{(2)}$ randomly.

For the new variant method based on random splitting, Fan, Samworth and Wu (2009) obtained some non-asymptotic probability bound for the event that $r$ unimportant covaraites are included in the intersection $\hat{\mathcal{I}}$ for any natural number $r$ under some exchangeability condition on all unimportant covariates. The probability bound is decreasing in the dimensionality, showing a “blessing of dimensionality”. Please consult Fan, Samworth and Wu (2009) for more details. We want to remark that their theoretical bound is applicable to our setting as well while studying time-to-event data because theoretical bound is based on splitting the sample into two halves and only requires the independence between these two halves.

While defining new variants, we may use the same $d$ as used in the original SIS and ISIS. However it will lead to a very aggressive screening. We call the corresponding variant the first variant of (I)SIS. Alternatively, in each step we may choose larger $\hat{\mathcal{I}}^{(1)}$ and $\hat{\mathcal{I}}^{(2)}$ to ensure that their intersection $\hat{\mathcal{I}}^{(1)}\cap\hat{\mathcal{I}}^{(2)}$ has $d$ covariates, which is called the second variant. The second variant ensures that there are at least $d$ covariates included before applying penalization in each step and is thus much less aggressive. Numerical examples will be used to explore their performance and prefer to the first variant.

In this section, we conduct simulation studies to show the power of the (I)SIS and its variants by comparing them with LASSO (Tibshirani, 1997) in the Cox’s proportional hazards model. Here the regularization parameter for LASSO is tuned via five fold cross validation. Most of the settings are adapted from Fan and Lv (2008) and Fan, Samworth and Wu (2009). Four different configurations are considered with $n=300$ and $p=400$. And two of them are revisited with a different pair of sample size $n=400$ and dimensionality $p=1000$. Covariates in different settings are generated as follows.

• Case

  1: $X_{1},\cdots,X_{p}$ are independent and identically distributed $N(0,1)$ random variables.

• Case

  2: $X_{1},\cdots,X_{p}$ are multivariate Gaussian, marginally $N(0,1)$, and with serial correlation $\mbox{corr}(X_{i},X_{j})=\rho$ if $i\neq j$. Here we take $\rho=0.5$.

• Case

  3: $X_{1},\cdots,X_{p}$ are multivariate Gaussian, marginally $N(0,1)$, and with correlation structure $\mbox{corr}(X_{i},X_{4})=1/\sqrt{2}$ for all $i\neq 4$ and $\mbox{corr}(X_{i},X_{j})=1/2$ if $i$ and $j$ are distinct elements of $\{1,\cdots,p\}\backslash\{4\}$.

• Case

  4: $X_{1},\cdots,X_{p}$ are multivariate Gaussian, marginally $N(0,1)$, and with correlation structure $\mbox{corr}(X_{i},X_{5})=0$ for all $i\neq 5$, $\mbox{corr}(X_{i},X_{4})=1/\sqrt{2}$ for all $i\notin\{4,5\}$, and $\mbox{corr}(X_{i},X_{j})=1/2$ if $i$ and $j$ are distinct elements of $\{1,\cdots,p\}\backslash\{4,5\}$.

• Case

  5: Same as Case 2 except $n=400$ and $p=1000$.

• Case

  6: Same as Case 4 except $n=400$ and $p=1000$.

Here, Case 1 with independent predictors is the most straightforward for variable selection. In Cases 2-6, however, we have serial correlation such that $\mbox{corr}(X_{i},X_{j})$ does not decay as $|i-j|$ increases. We will see later that for Cases 3, 4 and 6, the true coefficients are specially chosen such that the response is marginally independent but jointly dependent of $X_{4}$. We therefore expect variable selection in these situations to be much more challenging, especially for the non-iterated versions of SIS. Notice that in the asymptotic theory of SIS in Fan and Lv (2008), this type of dependence is ruled out by their Condition (4).

In our implementation, we choose $d=\lfloor\frac{n}{4\log n}\rfloor$ for both the vanilla version of SIS (Van-SIS) and the second variant (Var2-SIS). For the first variant (Var1-SIS), however, we use $d=\lfloor\frac{n}{\log n}\rfloor$ (note that since the selected variables for the first variant are in the intersection of two sets of size $d$, we typically end up with far fewer than $d$ variables selected by this method). For any type of SIS or ISIS, we apply SCAD with these selected predictors to get a final estimate of the regression coefficients at the end of the screening step. Whenever necessary, the BIC is used to select the best tuning parameter in the regularization framework.

In all setting, the censoring time is generated from exponential distribution with mean $10$. This corresponds to choosing the baseline hazard function $h_{0}(t)=0.1$ for $t\geq 0$. The true regression coefficients and censoring rate in each of the six cases are as follows:

• Case

  1: $\beta_{1}=-1.6328,\beta_{2}=1.3988,\beta_{3}=-1.6497,\beta_{4}=1.6353,\beta_{5}=-1.4209,\beta_{6}=1.7022$, and $\beta_{j}=0$ for $j>6$. The corresponding censoring rate is 33%.

• Case

  2: The coefficients are the same as Case 1. The corresponding censoring rate is 27%.

• Case

  3: $\beta_{1}=4,\beta_{2}=4,\beta_{3}=4,\beta_{4}=-6\sqrt{2}$, and $\beta_{j}=0$ for $j>4$. The corresponding censoring rate is 30%.

• Case

  4: $\beta_{1}=4,\beta_{2}=4,\beta_{3}=4,\beta_{4}=-6\sqrt{2},\beta_{5}=4/3$ and $\beta_{j}=0$ for $j>5$. The corresponding censoring rate is 31%.

• Case

  5: $\beta_{1}=-1.5140,\beta_{2}=1.2799,\beta_{3}=-1.5307,\beta_{4}=1.5164,\beta_{5}=-1.3020,\beta_{6}=1.5833$, and $\beta_{j}=0$ for $j>6$. The corresponding censoring rate is 23%.

• Case

  6: The coefficients are the same as Case 4. The corresponding censoring rate is 36%.

In Cases 1, 2 and 5 the coefficients were chosen randomly, and were generated as $(4\log n/\sqrt{n}+|Z|/4)U$ with $Z\sim N(0,1)$ and $U=1$ with probability 0.5 and -1 with probability 0.5, independent of $Z$. For Cases 3, 4, and 6, the choices ensure that even though $\beta_{4}\neq 0$, we have that $X_{4}$ and $Y$ are marginally independent. The fact that $X_{4}$ is marginally independent of the response is designed to make it difficult for the common independent learning to select this variable. In Cases 4 and 6, we add another important variable $X_{5}$ with a small coefficient to make it even more difficult.

We report our simulation results based on 100 Monte Carlo repetitions for each setting in Tables 1-7. To present our simulation results, we use several different performance measures. In the rows labeled $\|\mbox{\boldmath$\beta$}-\hat{}\mbox{\boldmath$\beta$}\|_{1}$ and $\|\mbox{\boldmath$\beta$}-\hat{}\mbox{\boldmath$\beta$}\|_{2}^{2}$, we report the median $L_{1}$ and squared $L_{2}$ estimation errors $\|\mbox{\boldmath$\beta$}-\hat{}\mbox{\boldmath$\beta$}\|_{1}=\sum_{j=1}^{p}|\hat{\beta}_{j}-\beta_{j}|$ and $\|\mbox{\boldmath$\beta$}-\hat{}\mbox{\boldmath$\beta$}\|_{2}^{2}=\sum_{j=1}^{p}|\hat{\beta}_{j}-\beta_{j}|^{2}$, respectively, where the median is over the 100 repetitions. In the row with label $P_{1}$, we report the proportion of the 100 repetitions that the (I)SIS procedure under consideration includes all of the important variables in the model, while the row with label $P_{2}$ reports the corresponding proportion of times that the final variables selected, after further application of the SCAD penalty, include all of the important ones. We also report the median model size (MMS) of the final model among 100 repetitions in the row labeled MMS.

We report results of Cases 1 and 2 in Table 1. Recall that the covariates in Case 1 are all independent. In this case, the Van-SIS performs reasonably well. Yet, it does not perform well for the dependent case, Case 2. Note the only difference between Case 1 and Case 2 is the covariance structure of the covariates. For both cases, vanilla-ISIS and its second variant perform very well. It is worth noticing that the ISIS improves significantly over SIS, when covariates are dependent, in terms of both the probability of including all the true variables and in reducing the estimation

error. This comparison indicates that the ISIS performs much better when there is serious correlation among covariates.

While implementing the LASSO penalized Cox’s proportional hazards model, we adapted the Fortran source code in the R package “glmpath.” Recall that the objective function in the LASSO penalized Cox’s proportional hazards model is convex and nonlinear. What the Fortran code does is to call a MINOS subroutine to solve the corresponding nonlinear convex optimization problem. Here MINOS is an optimization software developed by Systems Optimization Laboratory at Stanford University. This nonlinear convex optimization problem is much more complicated than a general quadratic programming problem. Thus generally it takes much longer time to solve, especially so when the dimensionality is high as confirmed by Table 3. However the algorithm we used does converge as the objective function is strictly convex.

Table 1 shows that LASSO has the sure screening property as the ISIS, however, the median model size is ten times as large as that of ISIS. As a consequence, it also has larger estimation errors in terms of $\|\mbox{\boldmath$\beta$}-\hat{}\mbox{\boldmath$\beta$}\|_{1}$ and $\|\mbox{\boldmath$\beta$}-\hat{}\mbox{\boldmath$\beta$}\|_{2}^{2}$. The fact that the median absolute deviation error is much larger than the median square error indicates that the LASSO selects many small nonzero coefficients for those unimportant variables. This is also verified by the fact that LASSO has a very large median model size. The explanation is the bias issue noted by Fan and Li (2001). In order for LASSO to have a small bias for nonzero coefficients, a smaller $\lambda$ should be chosen. Yet, a small $\lambda$ recruits many small coefficients for unimportant variables. For Case 2, the LASSO has a similar performance as in Case 1 in that it includes all the important variables but has a much larger model size.

Results of Cases 3 and 4 are reported in Table 2. Note that, in both cases, the design ensures that $X_{4}$ is marginally independent of but jointly dependent on $Y$. This special design disables the SIS to include $X_{4}$ in the corresponding identified model as confirmed by our numerical results. However, by using ISIS, we are able to select $X_{4}$ for each repetition. Surprisingly, LASSO rarely includes $X_{4}$ even if it is not a marginal screening based method. Case 5 is even more challenging. In addition to the same challenge as case 4, the coefficient $\beta_{5}$ is 3 times smaller than the first four variables. Through the correlation with the first 4 variables, unimportant variables $\{X_{j},j\geq 6\}$ have a larger marginal utility with $Y$ than $X_{5}$. Nevertheless, the ISIS works very well and demonstrates once more that it uses adequately the joint covariate information.

We also compare the computational cost of van-ISIS and LASSO in Table 3 for Cases 1-4. Table 3 shows that it takes LASSO several hours for each repetition, while van-ISIS can finish it in just several minutes. This is a huge improvement. For this reason, for Cases 5 and 6 where $p=1000$, we only report the results for ISIS since it takes LASSO over several days to complete a single repetition. Results of Cases 5 and 6 are reported in Table 4. The table demonstrates similar performance as Cases 2 and 4 even with more covariates.

To conclude the simulation section, we demonstrate the difficulty of our simulated models by showing the distribution, among 100 simulations, of the minimum $|t|$-statistic for the estimates of the true nonzero regression coefficients in the oracle model with only true important predictors included. More explicitly, during each repetition of each simulation setting, we pretend to know the index set $\mathcal{M}^{*}$ the true underlying sparse model, fit the Cox’s proportional hazards model using only predictors with indices in $\mathcal{M}^{*}$ by calling function “coxph” of R package “survival”, and report the smallest absolute value of the $t$-statistic for the regression estimates. For example, for case 1, the model size is only 6 and the minimum $|t|$-statistic is computed based on these 6 estimates for each simulation. This shows the difficulty to recover all significant variables even in the oracle model with the minimum model size. The corresponding boxplot for each case is shown in Figure 2. To demonstrate the effect of including unimportant variables, the minimum $|t|$-statistic for the estimates of the true nonzero regression coefficients is recalculated and shown by the boxplots in Figure 3 for the model with the true important variables and 20 unimportant variables.

As expected, cases 1 and 2 are relatively easy cases, whereas cases 3 and 4 are relatively harder in the oracle model. When we are not in the oracle model with 20 noisy variables are added, the difficulty increases as shown in Figure 3. It has more impact on cases 3, 4 and 6.