Variable selection and estimation for accelerated failure time model via seamless-$L_0$ penalty

1.   Introduction

Analyzing high-dimensional survival data has become an important topic in statistics, among which finding covariates with good predictive power of survival is a fundamental step. In the variable selection, penalized least squares procedures is an attractive approach. Penalized least squares procedures are used for variable selection and estimation which help predict estimators.

As a useful alternative to the Cox model ^[2], the AFT model^[10] based on linear regression models has an intuitive form compared to Cox model. The AFT model with an unspecified error distribution has been studied commonly for right-censored data. Two approaches in this aspect have gained attractive attention. One uses the Kaplan-Meier estimator to obtain the ordinary least squares estimator. The other is the rank-based estimator, which is motivated by the score function of the partial likelihood. See for examples in ^[1,13,17].

Identifying significant factors with predictive power, many techniques for linear regression models have been extended to the Cox regression and the AFT model. Penalized methods have drawn extensive attentions, which are for imposing some penalties to the regression coefficients. By balancing the goodness of fit and model complexity, penalization approaches lead the complex models to a profile. There exists plenty of methods used in gene expression analysis with survival data; see for examples in ^[16,19]. Moreover, various penalization methods of consistent selection have also been proposed. Examples include the adaptive Lasso ^[19], the smoothly clipped absolute deviations (SCAD) ^[6], the minimax concave penalty (MCP)^[20] and the bridge penalty. It has been shown that the bridge penalty had the oracle estimation in the linear regression models having divergent number of covariates. For the AFT models, there also exists much literature (e.g., ^[7,9,22]). To name but a few, Huang et al. ^[7] considered the regularization approaches for estimation in the AFT model with high-dimension covariates based on Stute's weighted least squares method. Huang and Ma ^[8] considered variable selection for AFT model with bridge method. Wang and Song^[18] applied adaptive lasso to the AFT models. In recent years, there are still a lot of studies on the AFT models. For example, Chai et al.^[3] considered a set of low-dimensional covariates of main interest and a set of high-dimensional covariates that may also affect survival under the accelerated failure time model. Choi and Choi ^[4] proposed the logistic-kernel smoothing procedure for the semi-parametric AFT model with high-dimensional right-censored data. Li et al.^[12] proposed a unified Expectation-Maximization approach combined with the $L_1$-norm penalty to perform variable selection and parameter estimation simultaneously in the accelerated failure time model with right-censored survival data of moderate sizes.

This article is motivated by the need for considering a seamless-$L_0$ (SELO) penalty ^[5] under the AFT model, which is a smooth function similar to the $L_0$ penalty. Under appropriate conditions, we show that the SELO selects a model whose dimension is comparable to the underlying model and prove that the proposed estimators is asymptotically normal. Monte Carlo simulations to evaluate the finite sample performance of the proposed procedure are computed. The proposed method is also demonstrated through an empirical analysis.

The rest of this paper is organized as follows. The AFT model is based on SELO penalization and computational algorithm are introduced in Section 2. In Section 3, we further propose an accurate variable selection for high dimensional sparse AFT model based on seamless $L_0$. The root $n$ consistency and the asymptotic normality of the resulting estimate are established. We simulate Monte Carlo simulation study to examine the finite sample performance of the proposed estimate in Section 4. A real data example is used to illustrate the proposed methodology in Section 5.

2.   SELO estimation in the AFT model

Let $T_{i}$ be the logarithm of the failure time and $X_{i}$ be the $p$-dimensional covariate vector. The AFT model assumes

where $\alpha$ is the intercept, $\beta\in\mathcal{R}^{p}$ is an unknown vector of interest, and ${\varepsilon}_{i}$ is the random error. When $T_{i}$ is subject to right censoring, we can only observe $(Y_{i}, \delta_{i}, X_{i})$, where $Y_{i} = \min(T_{i}, C_{i})$, $X_{i}$ be the $p$-dimensional covariate vector for the $i$th row of the $n\times p$ matrix $X$ which is the covariate matrix, $C_{i}$ is the logarithm of the censoring time, and $\delta_{i} = I\{T_{i}\le C_{i}\}$ is the censoring indicator. We assume that $(Y_{i}, \delta_{i}, X_{i})$, $i = 1, \ldots, n$, come from the same distribution.

Let $\hat{F}_{n}$ be the Kaplan-Meier estimator of the distribution function $F$ of $T$. $\hat{F}_{n}$ can be written as

where the $w_{i}$'s are the jumps in the Kaplan-Meier estimator expressed as $w_{1} = \frac{\delta_{(1)}}{n}$, $w_{i} = \frac{\delta_{(i)}}{n-i+1}\prod_{j = 1}^{i-1}\Big(\frac{n-j}{n-j+1}\Big)^{\delta_{(j)}}$, $j = 2, \ldots, n$. $w_{i}$'s are also called the Kaplan-Meier weights; see for examples in Stute and Wang^[14]. Here $Y_{(1)}\le\cdots\le Y_{(n)}$ are the order statistics of $Y_{i}$'s and $\delta_{(1)}, \ldots, \delta_{(n)}$ are the associated censoring indicators. Similarly, let $X_{(1)}, \ldots, X_{(n)}$ be the associated covariates of the ordered $Y_{i}$'s. The weighted least square(WLS) loss function is

Let $\bar{X}_{w} = \sum_{i = 1}^{n}w_{i}X_{(i)}/\sum_{i = 1}^{n}w_{i}$, and $\bar{Y}_{w} = \sum_{i = 1}^{n}w_{i}Y_{(i)}/\sum_{i = 1}^{n}w_{i}$, denoted by $X_{(i)}^{*} = (nw_{i})^{1/2}(X_{(i)}-\bar{X}_{w})$ and $Y_{(i)}^{*} = (nw_{i})^{1/2}(Y_{(i)}-\bar{Y}_{w})$. The weighted least square(WLS) objective function (2.2) can be written as

Penalized regression problem has been studied extensively. LASSO is one of the most popular and widely studied $L_{1}$ penalty. But it has been proved that its estimator may be inconsistent for model selection. The smoothly clipped absolute deviations (SCAD) and the minimax concave penalty (MCP) are another two popular penalties. SCAD is a continuous penalty and its estimator has oracle property. MCP also performs well in variable selection, whose estimator is consistent. However, $L_{0}$ penalty directly penalizes the non-zero parameters, whose drawback is the difficulty of computing because of its discontinuity. Seamless-$L_{0}$ (SELO) was proposed in Dicker ^[5], which was explicitly designed to minic $L_{0}$ penalty. It has been found that SELO possessed good theoretical properties.

We now describe the variable selection for AFT model via SELO. Coordinate descent is introduced to solve this problem. We propose tuning parameter $\lambda$ using cross-validation. The SELO penalized objective function is,

where SELO($\beta_{j}$) is defined as,

and $\lambda$ is tuning parameter. When $\lambda$ is large, SELO may select small estimators. In the paragraph, $\lambda$ is determined by Cross Validation. It is easy to see that when $\tau$ is enough small, $p_{SELO}(\beta_{j})\approx\lambda I\{\beta_{j}\neq0\}$, which is similar to $L_{0}$ penalty.

To minimize (2.3), we utilize coordinate descent algorithm. Coordinate descent algorithm^[21] has been widely used in penalized regression problem, which optimizes an objective function by calculating a single parameter at a time until convergence is reached. Dicker ^[5] described this algorithm for obtaining SELO estimators. The algorithm is formulated in terms of the tuning parameter $\lambda$. For a fixed value of $\lambda$, it can be implemented in the following steps.

3.   Theoretical properties

In this section, we prove the consistency and asymptotic normality of WLS estimator via SELO under some conditions. Following the notation of Stute ^[14,15], let $H$ denote the distribution function of $Y$. Under the assumption of independence between $T$ and $C$, $1-H(y) = (1-F(y))(1-G(y))$, where $F$ and $G$ are the distribution functions of $T$ and C. Let $\tau_{Y}$, $\tau_{T}$ and $\tau_{C}$ be the endpoints of the support of $Y$, $T$ and $C$. We put

Now, introduce the following sub-distribution functions:

Under random censoring the limit variance becomes much more complicated. Let

and

We assume that

(A1) (a) $E[(Y-X^{T}\beta^{*})^{2}XX^{T}\delta] < \infty$, (b) $\int|(w-x^{T}\beta^{*})x_{j}|D^{1/2}(w)\tilde F^{0}(dx, dw) < \infty$, for $j = 1, ..., p$ and $D(y) = \int_{0}^{y-}[(1-H(w))(1-G(w))]^{-1}G(dw)$.

(A2) $\lambda = O(1)$, $\tau = O(p^{1/2}/n^{3/2})$, $\lambda\sqrt{n/p}\rightarrow \infty$ and $p\sigma^{2}/n\rightarrow 0$ when $n\rightarrow \infty$.

(A3) $r \le\lambda(E(XX^{T}))\le R$, where $r$ and $R$ are positive constant.

(A4) $\lim_{n\rightarrow \infty} n^{-1}\max_{1\le i\le n}\sum_{j = 1}^{p}w_{i}^{2}x_{ij}^{2} = 0.$

(A5) $E\left(\Big|\frac{\epsilon}{\sigma}\Big|^{2+\frac{2\delta}{1+\delta}}\right) < M$ for some $M < 0$ and $\delta> 0$.

Condition (A1) is usually used for the proof of consistency in Stute^[14]. Condition (A2) restricts the size of $\lambda$ and $\tau$. Condition (A3) gives the bound of the eigenvalues of $E(XX^{T})$, which is used in Theorem 3.1. Conditions (A4) and (A5) are used for the proof the asymptotic normality of SELO estimators and are related to the Lindeberg condition of Lindeberg-Feller CLT.

Theorem 3.1. Suppose that conditions (A1)–(A5), then

(i) $\lim_{n\rightarrow \infty}P(\{j:\hat {\bf{ \pmb{\mathsf{ β}}}}_{j}\not = 0\} = A) = 1$, $A = \{j; \beta_{j}\neq 0\}$.

(ii) $\sqrt{n}(n^{-1}X_{A}^{T}WX_{A}/\sigma^{2})^{1/2}(\hat\beta_{A}-\beta_{A}^{*})\stackrel{D}\rightarrow (\sigma^{2}(X_{A}^{T}WX_{A}))^{-\frac{1}{2}}G_{A}$ where $G\sim N(0, \Sigma)$ and $\Sigma = Var\{\delta \gamma_{0}(Y)(Y-X\beta^{*})X+(1-\delta)\gamma_{1}(Y; \beta^{*})-\gamma_{2}(Y; \beta^{*})\}/n$ and $G_{A}$ is the part of G corresponding to $\beta_{A}^{*}$.

The proof is put in the supplementary.

4.   Numerical results

Let $T$ be generated from $T = X\tt\beta+\epsilon$, where $\epsilon\sim N(0, \sigma)$. The covariates $X = (X_{1}, ..., X_{p})$ are standard normal. Here we set $\sigma = 0.1$. The censoring variables are generated as uniformly distributed $U(0, C_{0})$ and independent of the events, where $C_{0}$ is chosen to obtain the censoring rate $25\%$ and $40\%$. Set two sample sizes $n = 200$ and $n = 400$. The tuning parameter $\lambda$ is chosen by cross validation. For each value of $n$, we simulated 1000 independent datasets $\{(y_{1}, x_{1}^{T}), ..., (y_{n}, x_{n}^{T})\}$. For each dataset, we calculated estimates of $\beta$. For each estimator $\hat\beta$, we recorded: the model size, $\hat{A} = \{j; \hat\beta_{j}\neq 0\}$; an indicator of whether or not the true model was selected, $I\{\hat A = A\}$; the false positive rate, $|\hat{A}-A|/|\hat{A}|$; the false negative rate, $|A-\hat{A}|/(p-|\hat A|)$; and the model error, $(\hat\beta-\beta^{*})^{T}(\hat\beta-\beta^{*})$. The column labeled "size", "rate", "F+", "F-" and "MSE" represent the above indicators. Results for SELO, LASSO, SCAD and MCP are summarized in the tables. Furthermore, we use the $V$-fold cross-validation to determine the tuning parameter. The CV score is $\sum_{v = 1}^{V}[\ell_n(\hat\beta^{(-v)})-\ell_{n}^{(-v)}(\hat\beta^{(-v)})].$ In this article, we set $V = 5$.

4.1.   Simulation I

The example was conducted with $p = 8$ and set $\beta = (3, 1.5, 0, 0, 1, 0, 0, 0)\in R^{8}$, where Table 1 summarizes the variable selection results based on SELO, LASSO, SCAD and MCP when censoring rates are $25\%$ and $40\%$. Overall, SELO performs better than other three methods, which selects the correct model more frequently. For instance, when the censoring rate is $25\%$, SELO selects the true model most accurately. The true model size is 3 and the average size from SELO is 3.43. LASSO performs worse than other three methods both in model size and correct rate. LASSO, SCAD and MCP select model with average size 4.12, 4.05 and 4.04 and select the correct model in $42\%$, $46.7\%$ and $47\%$. Similar results perform when the censoring rate is $40\%$. But we can easily see that the situation when the censoring rate is $25\%$ is better than that when the censoring rate is $40\%$. When $n$ increases, we can see that the results are better. For instance, SELO selects 3.07 variables when $n = 400$, which is more accurate compared to 3.43 when $n = 200$. Other indicators can also prove it.

4.2.   Simulation II

The example was conducted with $p = 50$ and set $\beta = (3, 1.5, 0, 0, 2, 0, 3, 0, 0, 2, 0, \ldots, 0)\in R^{50}$, the rest of $\beta$ is zero. Other settings were similar to the case in simulation I and the results are listed in Table 2. It is easy to see that SELO remains better performance compared to other three methods. The model size from SELO is 5.68, which is the closest to the true model when the censoring rate is $25\%$ compared to 9.64, 9.27 and 8.86 for LASSO, SCAD and MCP. And it also performs better in correct rate and other indicators. For instance, SELO selects $52\%$ correct variables which is better compared to $3\%$, $5\%$ and $7\%$ for LASSO, SCAD and MCP. The indicators "F+" and "F-" of SELO are the smallest among four methods. Similar results perform when $p = 50$ compared to $p = 8$. It can be concluded that the results are worse when the censoring rate is $40\%$. However, when $n$ increases from 200 to 400, SELO selects more correct models.

4.3.   Simulation III

The example is set under $25\%$ and $40\%$ censoring, and we also estimate the mean estimated variance across 1000 simulated datasets when $n$ is 200 and 400. From the Tables 3 and 4, we see that SELO with tuning over $\tau\in\{0.001, 0.01, 0.1, 0.5\}$ seems to give better variance when compared to SELO with $\tau = 0.01$ fixed. We can see that the 1000 estimators remain stable whenever censoring rate are $25\%$ and $40\%$ and the simulation results perform better when $n$ increases to 400.

4.4.   PBC Data

PBC data was collected in the Mayo Clinic trial of primary biliary cirrhosis of liver conducted between 1974 and 1984. A total of 424 PBC patients, referred to Mayo Clinic during that ten-year interval, met eligibility criteria for the randomized placebo controlled trial of the drug D-penicillamine. The first 312 cases in the data were participated in the randomized trial and contained complete data. The additional 112 cases did not participate in the clinical trial. After deleting the missing data, the remaining 276 datasets are used for the analysis. We consider 17 covariates: age, albumin, alk.phos, ascites, ast, bili, chol, copper, platelet, edema, hepato, protime, sex, spiders, stage, trt, trig.

Among 276 samples without losing data, we calculated Table 5. The optimal values of lambda with SELO is small, which is 0.008. And the optimal values of lambda that are chosen by CV for LASSO, SCAD and MCP are 0.011. LASSO selects 6 variables and SELO selects 3 variables. LASSO selects the genes including sex, heptato, bili, albumin, protime, stage. SCAD selects the same variables like MCP. However, SELO selects sex, albumin and protime. The variables selected by SELO are also contained by other three methods. Meanwhile, we also calculate the AIC (Akaike Information Criteria) $AIC = n\log(\hat\sigma^{2})+2(d+1)$ ($d$ is the non-zero parameter and $\hat\sigma$ is the error between the estimator and true value) for four methods which show that SELO results better with smaller AIC. In Table 6, we also calculated $p$ values for the coefficients of variables selected by SELO. We calculated $p$ values for several steps. Firstly, we calculated t values called $t_{k}$ where $t_{k} = \hat\beta/s$ and $s$ is the sum of squares error of the estimators. Secondly, we found the value of $t_{\alpha/2}(n-K)$ where $\alpha = 0.05$ and $n-K$ is the degree. Finally, we calculated the values of $P(T> t_{\alpha/2}(n-K)) = P(T < t_{\alpha/2}(n-K))$.

Table 6 indicates the significance test results. From the results, we can see that the $p$ value of sex, albumin and protime are all less than 0.05. The variables selected by SELO are all significant. Overall, SELO selected a simple model compared to another three methods.

5.   Conclusions

Statistical analysis of failure time with high dimension covariates is an important topic. In this article, we investigate a new method (SELO) for the AFT model with high dimension covariates, for simultaneous variable selection and estimation. A real dataset (PBC) is analyzed and SELO selects some important covariates. Our numerical results indicate that SELO performs better than another three methods. In this article, we address the situation where $p\ll n$ and prove the oracle property under the condition $p/n\rightarrow 0$. We allow both $n$ and $p$ to diverge but $p$ goes to infinity more slowly than $n$. The situation where $p$ is much larger than $n$ will be extended in the future research.

Acknowledgments

The authors would like to thanks the editors and three reviewers for their valuable and helpful comments. The authors also thanks Guangren Yang and Yiming Liu for their advice for revising the paper.

Conflict of interest

The authors declare no conflict of interest.

Supplementary

In order to complete the proofs of Theorem 3.1 (i), we first show two lemmas.

Lemma 1. Recall that

Then for every $r\in(0, 1)$, there exists a constant $C_{0}> 0$ such that

where $C\ge C_{0}$.

Proof. Let $\alpha_{n} = \sqrt{p\sigma^{2}/n}$ and fix $r\in(0, 1)$. To prove the Lemma 1, it suffices to show that if $C> 0$ is large enough, then

Furthermore, define $Q_{n}(u) = Q(\beta^{*}+C\alpha_{n}u)-Q(\beta^{*})$.Then,

By the results of Stute (1993, 1996), we have

where $W\sim N(0, \Sigma)$, with $\Sigma$ defined in the theorem. The last term in $Q_{n}(u)$ where $K(u) = \{j; P_{SELO}(\beta_{j}^{*}+C\alpha_{n}u_{j})-P_{SELO}(\beta_{j}^{*}) < 0\}$. The fact that $P_{SELO}$ is concave on $[0, \infty)$ imply that, for each $\beta$, $P_{SELO}(\beta_{j}^{*}+C\alpha_{n}u_{j})-P_{SELO}(\beta_{j}^{*})\ge -C\alpha_{n}|u_{j}|P_{SELO}'(\beta_{j}+C\alpha_{n}u_{j})$

Under conditions (A1) and (A3), for $I_{1}$,

For $I_{2}$,

For $I_{3}$, by condition (A2), we have $I_{3} = o_{p}(C\alpha_{n})$. We conclude that if $C> 0$ is large enough, then $\inf_{\|u\| = 1}Q_{n}(u)> 0$ holds for all $n$ sufficiently large, with probability at least $1-r$. This finishes the proof of the Lemma 1.

Lemma 2. We assume that $C> 0$, $Q(\beta)$ is similar to that in Lemma 1. Under the conditions (A1)–(A3),

where $A^{c} = \{1, ..., p\}/A$ is the complement of $A$ in $\{1, ..., p\}$.

Proof. Suppose that $\beta\in R^{p}$ and that $\|\beta-\beta^{*}\| < C\alpha_{n}u$. Define $\tilde\beta\in R^{p}$ by $\tilde {\bf{ \pmb{\mathsf{ β}}}}_{A^{c}} = 0$ and $\tilde\beta_{A} = \beta_{A}$. Similar to the proof of Lemma 1, if $D_{n}(\beta, \tilde\beta) = Q_{n}(\beta)-Q_{n}(\tilde\beta)$, then

For $I_{1}$ and $I_{2}$, under the conditions (A1) and (A3),

For $I_{3}$, $P_{SELO}$ is concave $\|\beta\| < C\sqrt{p\sigma^{2}/n}\rightarrow |\beta_{j}| < C\sqrt{p\sigma^{2}/n}$, and

Under the condition (A2), we combine the results of $I_{1}$, $I_{2}$ and $I_{3}$. So, $D_{n}(\beta, \tilde\beta)> 0$.

Combining the proof of Lemmas 1 and 2, we can have the conclusion of Theorem 3.1 (i).

Proof of Theorem 3.1 (ii). We consider the proof related to the Lindeberg condition of Lindeberg-Feller CLT. Under the conditions (A1)–(A5), let $\hat\beta_{A}$ be a estimator where $A = \{j; \hat\beta_{j}\neq 0\}$.

We can easily have the following form,

and

To prove,

where $w_{i, n} = \left(\sigma^{2}X_{A}^{T}WX_{A}\right)^{-\frac{1}{2}}w_{i}x'_{(i), A}\epsilon_{i}.$ Let $\eta_{i, n} = w_{i}x_{(i), A}^{T}\left(X_{A}^{T}WX_{A}\right)^{-\frac{1}{2}}(X_{A}^{T}WX_{A})^{-\frac{1}{2}}w_{i}x_{(i), A}.$

By the condition of Lindeberg-Feller CLT, we have

By Holder inequality, set $\frac{1}{p} = \frac{2}{2+\delta}$, $\frac{1}{q} = \frac{\delta}{\delta+2}$,

By Markov inequality,

and

We showed that $\sum_{i = 1}^{n}\eta_{i, n} = \sum_{i = 1}^{n}w_{i}$ and $\eta_{i, n}\le \|(n^{-1}X_{A}^{T}WX_{A})^{-\frac{1}{2}}\|^{2}\underset{1\le i\le n}{\max}\sum_{j = 1}^{q}\frac{1}{n}w_{i}^{2}x_{ij}^{2}$

By conditions (A4) and (A5), $\sum_{i = 1}^{n}E[\|w_{i, n}\|^{2};\|w_{i, n}\|^{2}>\delta_{0}]\rightarrow 0$.

By conditions (A1)–(A3), $\|\left(\sigma^{2}\sum_{i = 1}^{n}X_{A}^{T}WX_{A}\right)^{-1/2}p'_{A}(\hat\beta)\| = o_{p}(1)$,

and

where, $\Sigma = Var\{\delta \gamma_{0}(Y)(Y-X\beta^{*})X+(1-\delta)\gamma_{1}(Y; \beta^{*})-\gamma_{2}(Y; \beta^{*})\}$.

So, $(\sigma^{2}X_{A}^{T}WX_{A})^{-\frac{1}{2}}X_{A}^{T}W\epsilon\stackrel{D}\rightarrow (\sigma^{2}X_{A}^{T}WX_{A})^{-\frac{1}{2}}G_{A}$, where $G\sim N(0, \Sigma)$ follows.