Skip to main content
NIHPA Author Manuscripts logoLink to NIHPA Author Manuscripts
. Author manuscript; available in PMC: 2013 Jul 9.
Published in final edited form as: Stat Sin. 2013 Jan 1;23(1):145–167. doi: 10.5705/ss.2011.100

VARIABLE SELECTION FOR CENSORED QUANTILE REGRESION

Huixia Judy Wang 1, Jianhui Zhou 2, Yi Li 3
PMCID: PMC3705953  NIHMSID: NIHMS353982  PMID: 23847407

Abstract

Quantile regression has emerged as a powerful tool in survival analysis as it directly links the quantiles of patients’ survival times to their demographic and genomic profiles, facilitating the identification of important prognostic factors. In view of limited work on variable selection in the context, we develop a new adaptive-lasso-based variable selection procedure for quantile regression with censored outcomes. To account for random censoring for data with multivariate covariates, we employ the ideas of redistribution-of-mass and e ective dimension reduction. Asymptotically our procedure enjoys the model selection consistency, that is, identifying the true model with probability tending to one. Moreover, as opposed to the existing methods, our new proposal requires fewer assumptions, leading to more accurate variable selection. The analysis of a real cancer clinical trial demonstrates that our procedure can identify and distinguish important factors associated with patient sub-populations characterized by short or long survivals, which is of particular interest to oncologists.

Keywords: Conditional Kaplan-Meier, dimension reduction, kernel, quantile regression, survival analysis, variable selection

1. Introduction

Quantile regression, as a valuable alternative to the commonly used Cox proportional hazard model and accelerated failure time (AFT) model (Koenker and Geling, 2001; Portnoy, 2003), directly links the quantiles of subjects’ survival times to their demographic and genomic profiles, and thus can facilitate the identification of important prognostic factors. Direct applications of this model lie in, for example, cancer studies, where physicians are often interested in identifying effective treatments for more severe cases (with shorter survival times). As might be expected, treatments often cause different impacts among patients that fall within the upper or lower quantiles of the survival distribution. It is well known that both Cox and AFT models restrict the covariates to a ect only the location but not the shape of the survival distribution, and thus may overlook interesting forms of heterogeneity. For instance, these two models do not permit the treatment effect to be positive for severe cases while negative for the other cases. In contrast, quantile regression offers a convenient approach to capture the variation caused by heterogeneities by allowing the covariates to exhibit different impacts at different tails of the survival distribution.

Current literature on quantile regression for censored survival outcomes often focus on coefficient estimation (Portnoy, 2003; Peng and Huang, 2008; Huang, 2009; Wang and Wang, 2009), but little effort has been devoted to variable selection. In contrast, various penalization-based variable selection methods have been developed for Cox and AFT models, for instance, Huang et al. (2006), Zhang and Lu (2007), Wang et al. (2008), Engler and Li (2009), among others.

The nature of the existing estimating procedures for censored quantile regression prohibits the direct usage of the popular penalization methods. The point estimation methods of Portnoy (2003) and Peng and Huang (2008) require fitting an entire quantile process, as the estimation at an upper quantile depends on the estimations at all the lower quantiles. Therefore, even if our main interest was to identify variables with strong impacts on the median survival, we would have to estimate and select variables for all the lower quantiles as well. This imposes both computational and theoretical challenges. Alternative estimation methods such as those proposed by Subramanian (2002), Wang and Wang (2009) and Leng and Tong (2011) rely on the kernel-smoothing estimation, and thus are practically feasible only for data with few covariates. Recently, Shows, Lu and Zhang (2010) developed a variable selection approach for censored median regression by using the inverse-probability-weighting scheme of Bang and Tsiatis (2002). However, the method requires the restrictive unconditional independence assumption between survival and censoring times. Moreover, the procedure uses information of only uncensored observations, leading to an efficiency loss in the estimation.

We develop a new and flexible variable selection method based on the adaptive- lasso penalization for censored quantile regression. Our work advances the field in the following three significant ways. First, our method adapts the redistribution of-mass idea of Efron (1967) to account for the censoring in quantile regression. different from existing estimation procedures (Portnoy, 2003; Wang and Wang, 2009), the proposed method estimates the masses for redistribution by using the conditional Kaplan-Meier estimation on a reduced data space. Therefore, the new method is more flexible and is able to accommodate high-dimensional covariates. Secondly, our variable selection procedure enjoys computation readiness and requires much fewer stringent assumptions than those in the literature. As a result, our procedure is model selection consistent under mild conditions and gives better finite sample performance than the existing methods. Lastly, the proposed method is able to capture important heterogeneities in the survival population, and to identify effective biomarkers that have impacts at different tails of the survival distribution. Such results will be of particular interest to physicians who are keen on designing effective treatments for targeted patient subpopulations, often characterized by short survivals.

2. Variable Selection for Censored Quantile Regression

2.1. Model setup

Denote Ti as the uncensored survival outcome, xi as the observable p-dimensional covariates, and Ci as the censoring variable. Consider the following quantile regression model

Ti=xiTβ0(τ)+ei(τ), (2.1)

where 0 < τ < 1 is the given quantile level of interest, β0(τ) is the p-dimensional unknown quantile coefficient vector, and ei(τ) is the random error whose τth conditional quantile given xi equals zero. Without loss of generality, we assume that the first element of xi is 1 corresponding to the intercept. In practice, we only observe (xi, Yi, δi), where Yi = min(Ti, Ci) is the observed response variable and δi = I(Ti ≤ Ci) is the censoring indicator. Our main objective in this paper is to select important predictors that have nonzero effect on the τth conditional quantile of T in the parametric quantile regression model (2.1).

2.2. Variable selection via redistribution-of-mass

We first briefly review the idea of redistribution-of-mass in censored quantile regression. For censored data without covariates, that is, p = 1 in (2.1), Efron (1967) proposed a simple algorithm for deriving the Kaplan-Meier estimator by redistributing the mass of each censored observation uniformly to observations on the right. In regression setup, this means to redistribute the probability masses P(Ti > Ci|Ci, xi) of censored cases to observations on the right.

Let F0(t|x) = P(T < t|x) denote the conditional distribution function of T given x. Define π0i = F0(Ci|xi) as the conditional probability for the ith subject not to be censored. Consider an ideal scenario where π0i are known. Then β0(τ) can be estimated by minimizing the following weighted quantile objective function with respect to β,

L(β,w0)=∑i=1n{w0iρτ(Yi−xiTβ)+(1−w0i)ρτ(Y+∞−xiTβ)}, (2.2)

where ρτ(u) = u{τ – I(u < 0)}, Y+∞ is any value su ciently large to exceed xTi β0(τ) for all i, and

w0i={1δi=10δi=0andπ0i>τ(τ−π0i)∕(1−π0i)δi=0andπ0i≤τ}. (2.3)

It can be shown that a subgradient of the above weighted objective function,

nMn(β,w0)=∑i=1nxi{1−w0iI(Yi<xiTβ)}, (2.4)

is an unbiased estimating function of β0(τ). Therefore, minimizing L(β, w0) with respect to β leads to a consistent estimator of β0(τ). More explanation of the intuition behind the above weighting scheme in (2.2) can be found in Wang and Wang (2009), and Portnoy and Lin (2010).

To select variables, we consider the following penalized objective function

LAL(β,w0)=L(β,w0)+λn∑j=1pνj∣βj∣, (2.5)

where λn is the positive penalization parameter and νj are the adaptive weights. The type of adaptive lasso penalization was first proposed by Zou (2006) for least squares regression and later extended to quantile regression for uncensored data by Wu and Liu (2009) and to Cox's model by Zhang and Lu (2007). The adaptive lasso assigns heavier penalties to the potentially irrelevant variables, so the corresponding effects are shrunk more towards zero. This approach leads to a sparse coefficient estimation, and thus provides a convenient way to conduct model fitting and variable selection simultaneously. The choice of the adaptive weights νj is explained in (2.8) of Section 2.3.

2.3. Estimation of the redistributed mass

In practice, the masses for redistribution, 1–π0i = 1–F0(Ci|xi), are unknown and a variety of attempts have been made to estimate them. For example, Portnoy (2003) proposed to estimate π0i through fitting an entire quantile regression process under the global linearity assumption of the conditional quantile functions. McKeague et al. (2001) suggested to fit a semiparametric regression model such as Cox proportional hazards model to obtain an approximation of π0i. Lindgren (1997), Subramanian (2002) and Wang and Wang (2009) employed a fully nonparametric approach based on the conditional Kaplan-Meier estimator of F0(·|x). Though flexible, this nonparametric approach is only feasible when the covariate dimension is small because of the curse of dimensionality. The theoretical results in Subramanian (2002) and Wang and Wang (2009) were developed only for cases with univariate covariate.

We propose an index-based procedure to obtain a nonparametric estimation of π0i for multivariate covariates. The main idea is to summarize the regression information contained in x by indices through dimension reduction. Specifically, we adopt a global dimension reduction (DR) formulation:

Ti⫫xi∣(xiTγ1,⋯,xiTγq), (2.6)

where ⫫ stands for independence. This formulation stipulates that the dependence of Ti on the p-dimensional xi only comes from q indices, zi,1=xiTγ1,⋯,zi,q=xiTγq, where q is often much smaller than p in practice. For randomly censored data, γj, often referred to as effective dimension reduction (EDR) directions, can be estimated by using the sliced inverse regression (SIR) method of Li et al. (1999) or the hazard-function-based minimum average variance estimation (MAVE) method of Xia et al. (2010). Under some regularity assumptions, both methods lead to estimators γ^j that are root-n consistent to γ0,j, where γ0,j∈Rp,j=1,⋯,q, is a set of EDR directions. Hereafter, we denote the estimated indices as z^i=(z^i,1,⋯,z^i,q) with z^i,j=xiTγ^j, and denote z0i = (z0i,1, · · · , z0i,q) with z0i,j=xiTγ0,j,j=1,⋯,q.

Under model (2.6), we have F0(t|xi) = F0(t|z0i) for any t and i. We then proceed to use Beran's local Kaplan-Meier estimator F^(⋅∣z) (Beran,1981) to estimate F(·|z). Specifically,

F^(t∣z)=1−∏j=1n{1−Bnj(z)∑k=1nI(Yk≥Yj)Bnk(z)}ηj(t), (2.7)

where ηj(t)=IYj≤t,δj=1),Bnk(z)=Kq(z−zkhn)∕∑i=1nKq(z−zihn), hn is the bandwidth, and Kq(z−zihn)=Kq(z1−zi,1hn,⋯,zq−zi,qhn). In this paper, we adopt the commonly used product kernel function Kq(u1,…,uq)=∏i=1qK(ui), where K(·) is a univariate kernel function. Here we opt for Beran's estimator as it is nonparametric and thus flexible, avoiding estimating the entire quantile process assuming global linear models as in Portnoy (2003). Therefore, π0i can be estimated by π^i=F^(Ci∣z^i), the nonparametric estimation of the conditional distribution of T given the indices with a reduced dimension q.

For variable selection, we define the proposed penalized estimator β^(τ) for β0(τ) in model (2.1) as the minimizer of the following penalized objective function

LAL(β,w^)=∑i=1n{w^iρτ(Yi−xiTβ)+(1−w^i)ρτ(Y+∞−xiTβ)}+λn∑j=1pνj∣βj∣, (2.8)

where w^i are the weights for redistribution-of-mass as defined by replacing π0i with π^i in (2.3), and νj are the adaptive weights. We let νj=∣β~j(τ)∣−r, where β~j(τ) is the jth element of the initial consistent estimator of β(τ). We define the initial estimator β~(τ) as the unpenalized estimator, that is, the minimizer of LAL(β,w^) with λn = 0. In our implementation, we choose r = 2.

We further stress that our aim of the paper is to select the important predic tors that have nonzero effect on the τth conditional quantile of T in the parametric quantile regression model (2.1). This variable selection is for the local quantiles of T and thus is different from the global nonparametric dimension reduction in the formulation (2.6). In addition, the proposed penalization procedure does not require the unconditional independence between the survival times Ti and censoring times Ci, which is a significant improvement compared to Shows et al. (2010).

2.4. Computation and tuning

The proposed procedure requires choosing the bandwidth parameter hn. Our empirical experience suggests that the performance of the proposed procedure is not sensitive to the choice of hn (see Section 3). In practice, we can use K-fold cross validation to choose the bandwidth. We first divide the data set randomly into K parts with roughly equal size. For the kth part, k = 1, · · · , K, we fit model (2.1) using the rest K – 1 parts of the data, and then evaluate the quantile loss from predicting the τth conditional quantile of T for the uncensored data that are left out. For computational simplicity, we use the unpenalized estimator β~(τ) with λn = 0 when calculating the quantile loss. We choose the hn that gives the minimum average quantile loss.

The proposed estimation also involves the penalization parameter λn, which determines the sparseness of the resulting estimator. In practice, the penalization parameter is often selected by minimizing some model-selection criterion. One commonly used criterion in variable selection literature is the Bayesian Information Criterion (BIC, Schwarz, 1948), which provides a large-sample approximation to twice the logarithm of the Bayes factor. Specifically, the BIC is defined as

BIC=−2{logL(β^R)−logL(β^F)}+(pR−p)logn,

where L(β^R) and L(β^F) are the maximized likelihoods under a reduced model with pR parameters and under the full model with p parameters, respectively. It is known that Rao's score test statistics (Rao, 1948) are asymptotically equivalent to the likelihood ratio statistics under both null and Pitman's alternative hypotheses (Serfling, 1980, page 156). Koenker and Machado (1999) discussed tests based on these two types of statistics in linear quantile regression setup. Motivated by this, we propose to choose λn that minimizes the following Score-based Bayesian Information Criterion (SBIC):

SBIC(λn)=nMn(β^λn(τ),w^i}Dn−1Mn{β^λn(τ),w^i}+pλnlog(n), (2.9)

where β^λn(τ) is the penalized estimator with the penalization parameter value of λn, pλn is the number of non-zero elements in β^λn(τ), and Dn=n−1∑i=1nxixiT{τw^i2+τ2(1−2w^i)} is the asymptotic covariance matrix of the subgradient based on the plugged-in weights w^i. Similar score-based information criterion was also employed and justified by Leng (2010) for regularized rank regression.

In addition, to obtain the weights for redistribution-of-mass, we need determine the number of indices q. Li et al. (1999) proposed a chi-squared test for determining the number of significant EDR directions obtained by SIR. Xia et al. (2010) developed an alternative selection criterion. This method was shown to be consistent for selecting q but it is based on cross-validation and thus is computationally more intensive than the chi-squared test of Li et al. (1999).

2.5. Asymptotic Properties

To establish the asymptotic results in this paper, we require the following assumptions.

  • A1 The random design vector x is bounded in probability, has a bounded density function, and E(xxT) is a positive definite p × p matrix.

  • A2 Let F0(t|x) and G(t|x) denote the survival functions of Ti and Ci conditional on x, respectively. The first derivatives of F0(t|x) and G(t|x) with respect to t, denoted as f0(t|x) and g(t|x) are uniformly bounded with respect to t and x. The F0(t|x) and G(t|x) have bounded (uniformly in t) second-order partial derivatives with respect to x. In addition, supt|F0(t|x′) – F0(t|x)| = O(||x ′ – x||), where || · || denotes the Euclidean norm.

  • A3 The true coefficient β0(τ) is in the interior of a bounded convex region B. For β in a neighborhood of β0(τ), E{xxT f0(xT β0|x)} is positive definite, and 1 – G(xT β|x) = P (C > xT β|x) > 0 with probability one.

  • A4 There exists an EDR direction γ0,j∈Rq such that for any j = 1, · · · , q, (i) γ^j−γ0,j=Op(n−1∕2); (ii) n−1∕2(γ^k−γk)=n−1∑i=1ndki, where dki are independent p-dimensional with means zero and finite variances.

  • A5 The univariate nonnegative kernel function K(·) has a compact support. It is a νth order kernel function satisfying ∫K(u)du=1,∫K2(u)du<∞,∫ujK(u)du=0 for j < ν and ∫∣u∣vK(u)du<∞, and it is Lipschitz continuous of order ν, where ν ≥ 2 is an integer.

  • A6 The bandwidth hn satisfies hn = O(n–α) with (i) 0 < α < min(1/ν , 1/q); (ii) 1/(2ν) < α < 1/(3q) and ν > 3/2q.

Remark 1 The boundedness condition of x in A1 is posed for technical convenience. It is possible to allow the bound on x to grow slowly in n but this will complicate the technical proof. Assumption A2 is needed to obtain the asymptotic properties of the local Kaplan-Meier estimator. Assumption A3 ensures the identifiability of β0(τ). Assumption A4(i) states the root-n consistency of the estimated EDR direction, which is needed to establish the root-n consistency of β~(τ). This assumption holds for the modified sliced inverse regression estimation in Li et al. (1999) and for the hazard-function-based minimum average variance estimation in Xia et al. (2010). The linear presentation of γ^k is assumed in A4(ii) to help establish the normality of β~(τ), and this condition can be obtained for the dimension reduction methods of Li et al. (1999) and Xia et al. (2010) with more technical endeavors under some higher-level conditions. Assumption A5 requires K(·) to be a νth order kernel, where the requirement of depends on the dimensionality of indices q. For larger q, higher order kernel function is needed in order to control the bias; see Hu and Fan (1992) for discussions on the construction of higher order kernel functions. Assumption A6 specifies the conditions on the bandwidth hn, where the weaker condition (i) is needed for establishing the consistency and the stronger condition (ii) is needed for the normality.

We first establish the consistency and asymptotic normality of the initial unpenalized estimator β~(τ).

Theorem 1 Suppose models (2.1), (2.6) and assumptions A1-A4(i), A5, A6(i) hold, then β~(τ)→β0(τ) in probability as n → ∞. Furthermore, if A4(ii) and A6(ii) both hold, then we have n1∕2{β~(τ)−β0(τ)}→DN(0,Γ1−1VΓ1), where Γ1=E[xxT{1−G(xTβ0(τ)∣x)}f0{xTβ0(τ)∣x}], V = cov(vi) with vi defined in (6.9) of the Appendix.

We next establish the property of consistency in variable selection of the proposed penalized estimator β^(τ). Denote A(τ)={j:βj(τ)≠0} and Ac(τ)={j:βj(τ)=0}.

Theorem 2 Suppose models (2.1), (2.6) and assumptions A1-A6(ii) hold, n–1/2λn → 0 and nr/2 – 1λn → ∞, then P({j:β^j(τ)≠0}=A(τ))→1 as n → ∞.

Remark 2 Theorem 2 suggests that the proposed procedure is able to select the correct model with probability approaching one. To achieve the same efficiency as the oracle estimator obtained under the true model, we can update the estimates for the non-zero coefficients in β^(τ) by minimizing the weighted objective function (2.8) using the selected covariates only with λn = 0. By Theorem 1, the updated estimator is asymptotically normal and it has the same efficiency as the oracle estimator. For finite samples, our numerical studies show that the re-estimation helps reduce the estimation bias of non-zero coefficients caused by shrinkage and thus leads to more e cient estimation.

3. Simulation Study

We study two different examples to investigate the performance of the proposed penalized estimator via redistribution of mass, referred to as PROM. The dimension of covariates is 20 in Example 1 and 100 in Example 2. We focus on two quantile levels τ = 0.25 and 0.5, and sample sizes n = 200 and 500. For each scenario, the simulation is repeated 500 times.

Example 1. The survival times are generated from the following model

Ti=1+1.5xi1+0.7xi2+xi3−0.5xi4+(1+γxi4)∊i,

where i = 1, · · · , n, εi ~ N(0, 1), and γ measures the heteroscedasticity. Besides the 4 relevant covariates, we include another 16 independent noise variables, xi5, · · · , xi20. For j = 1, · · · , 20, xij ~ U(–1, 1). We set γ=-0.742, so that the quantile coefficients are (0.326, 1.5, 0.7, 1.0, 0, · · · , 0) at τ = 0.25 and (1.0, 1.5, 0.7, 1.0, –0.5, 0, · · · , 0) at τ = 0.5. Under this heteroscedastic model, the covariate xi4 has a negative impact on the median but no impact on the first quartile of the conditional distribution of Ti. The observed responses are obtained by Yi = min(Ti, Ci), where Ci ~ U(–2, 18) yielding averaged 15% censoring and Ci ~ U(–2, 8) yielding averaged 30% censoring.

We compare the following estimators. The PROMS and PROMM are two variations of the proposed PROM estimator, where the indices are estimated by using the sliced inverse regression (SIR) estimation of Li et al. (1999) and the minimum average variance estimation (MAVE) of Xia et al. (2010), respectively. The SIR estimation is obtained by using the R function implemented by Sun available at http://www.bios.unc.edu/~wsun/, and the MAVE estimation is obtained by using the matlab program provided by Xia. The oracle estimator is obtained by using the proposed unpenalized method under the true model, that is, with the first three covariates at τ = 0.25 and the first four covariates at τ = 0.50. The oracle estimator serves as a gold standard. Two variations of the oracle estimator, OracleM and OracleS are included corresponding to MAVE and SIR indices estimation, respectively. The PIPW is the penalized estimator developed by Shows et al. (2010) by using the inverse-probability-weighting scheme of Bang and Tsiatis (2002).

For both examples, the number of true indices is q = 2. Therefore, for the PROM and oracle estimators, we use the fourth-order kernel function (Müller, 1984): K(x) = (105/64)(1 – 5x2 + 7x4 – 3x6)I(|x| ≤ 1). The bandwidth hn is selected by 5-fold cross validation as described in Section 2.4. Our numerical studies suggest that the proposed estimator PROM based on the selected q̂ performs very similarly as that based on the the true q; see Table 4 for comparison in Example 2. For computational convenience, in Example 1, we use q = 2 for both MAVE and SIR estimation.

Table 4.

Variable selection results for Example 2.

15% censoring
30% censoring
τ = 0.25
τ = 0.50
τ = 0.25
τ = 0.50
Method TP FP OP TP FP OP TP FP OP TP FP OP
n = 200
PIPW 4.99 39.64 0.00 5.00 25.16 0.01 4.98 68.70 0.00 4.99 62.65 0.00
PROMS (q) 4.91 0.34 0.68 4.97 0.34 0.70 4.78 0.48 0.52 4.85 0.38 0.60
PROMS (q̂) 4.91 0.40 0.64 4.97 0.27 0.76 4.77 0.53 0.50 4.85 0.43 0.58

n = 500
PIPW 5.00 5.04 0.18 5.00 3.19 0.24 5.00 11.37 0.08 5.00 9.48 0.06
PROMS (q) 5.00 0.13 0.90 5.00 0.14 0.89 5.00 0.17 0.88 5.00 0.10 0.93
PROMS (q̂) 5.00 0.16 0.87 5.00 0.13 0.90 5.00 0.19 0.86 5.00 0.09 0.92

TP denotes the average number of relevant variables that are correctly selected. FP denotes the average number of irrelevant variables that are incorrectly selected. OP denotes the oracle proportion, that is, the percentage of times that the true model is correctly selected. The methods PROMS (q) and PROMS (q̂) are based on the true and the estimated number of indices, respectively.

Table 1 summarizes the variable selection results of the penalized estimators PROMM, PROMS and PIPW. Overall, the proposed PROM methods with indices by both MAVE and SIR outperform the PIPW method in variable selection. The PROM and PIPW methods perform comparably for selecting the relevant variables (TP), but the PIPW selects irrelevant variables more often and thus has much lower oracle proportions in all scenarios considered. As n increases to 500, the PROM methods have oracle proportions close to 1 at both quantiles.

Table 1.

Variable selection results for Example 1.

15% censoring
30% censoring
τ = 0.25
τ = 0.50
τ = 0.25
τ = 0.50
Method TP FP OP TP FP OP TP FP OP TP FP OP
n = 200
PIPW 3.96 1.72 0.31 4.92 1.33 0.36 3.94 2.55 0.19 4.89 2.30 0.19
PROMM 3.95 0.19 0.80 4.64 0.12 0.58 3.92 0.16 0.78 4.45 0.12 0.41
PROMS 3.95 0.15 0.82 4.66 0.15 0.58 3.91 0.15 0.78 4.47 0.10 0.41

n = 500
PIPW 4.00 0.71 0.58 5.00 0.52 0.68 4.00 1.27 0.40 5.00 1.07 0.46
PROMM 4.00 0.08 0.92 5.00 0.08 0.91 4.00 0.09 0.92 4.91 0.05 0.87
PROMS 4.00 0.09 0.92 5.00 0.09 0.90 4.00 0.09 0.92 4.90 0.07 0.83

TP denotes the average number of relevant variables that are correctly selected. FP denotes the average number of irrelevant variables that are incorrectly selected. OP denotes the oracle proportion, that is, the percentage of times that the true model is correctly selected.

Table 2 summarizes the mean squared errors of estimates for the non-zero quantile coefficients from all methods. Compared to the other coefficients, β2(0.25) = β2(0.5) = 0.7 and β4(0.5) = –0.5 have smaller magnitudes. When n = 200, the proposed methods PROMM and PROMS tend to over shrink these small coefficients to trade for simpler models, which results in larger mean squared errors than PIPW. However, when n increases to 500, the PROM estimates become more efficient than PIPW , and their mean squared errors become comparable to those of the oracle estimates.

Table 2.

Mean squared errors (×100) of the estimates for nonzero quantile coefficients in Example 1.

τ = 0.25
τ = 0.50
β0 (τ) β1 (τ) β2 (τ) β3 (τ) β0 (τ) β1 (τ) β2 (τ) β3 (τ) β4 (τ)
n = 200, 15% censoring
OracleM 0.67 1.95 2.05 1.79 1.14 1.75 1.69 1.49 2.20
OracleS 0.68 1.92 2.02 1.79 1.12 1.76 1.69 1.54 2.24
PIPW 1.46 2.38 3.33 2.56 1.50 2.01 2.61 2.19 4.79
PROMM 0.77 2.01 4.27 1.85 1.81 1.95 3.25 1.72 9.12
PROMS 0.78 2.03 4.20 1.89 1.76 1.86 2.80 1.65 8.80

n = 200, 30% censoring
OracleM 0.78 2.29 2.36 2.03 1.54 2.06 1.94 1.67 2.83
OracleS 0.81 2.24 2.35 2.06 1.66 1.97 2.00 1.71 3.10
PIPW 1.75 3.08 3.54 3.03 1.80 2.89 3.51 2.89 6.00
PROMM 1.02 2.42 5.95 2.24 2.70 2.32 5.22 2.00 12.86
PROMS 1.08 2.42 6.32 2.12 2.80 2.14 4.91 2.26 13.36

n = 500, 15% censoring
OracleM 0.23 0.75 0.71 0.79 0.44 0.63 0.64 0.71 1.08
OracleS 0.23 0.73 0.72 0.78 0.48 0.62 0.64 0.72 1.12
PIPW 0.48 0.84 1.01 0.96 0.51 0.68 0.87 0.83 1.60
PROMM 0.26 0.75 0.70 0.81 0.49 0.65 0.65 0.71 1.55
PROMS 0.26 0.73 0.71 0.79 0.49 0.64 0.64 0.73 1.63

n = 500, 30% censoring
OracleM 0.27 0.80 0.77 0.88 0.84 0.76 0.67 0.81 1.59
OracleS 0.27 0.79 0.76 0.87 0.88 0.75 0.67 0.79 1.66
PIPW 0.56 0.98 1.10 1.10 0.65 0.91 1.10 1.00 2.03
PROMM 0.30 0.82 0.81 0.89 0.92 0.77 0.69 0.83 3.23
PROMS 0.31 0.78 0.77 0.90 1.03 0.75 0.68 0.84 3.59

To study the sensitivity of the developed PROM method to the bandwidth hn, we choose h = cn–0.15 and apply PROM to the same data sets used in Tables 1 and 2 for c = 0.2, 0.4, 0.6,..., 2.0. The results show that the developed PROM method is very robust to the bandwidth hn in both variable selection and parameter estimation. To save space, we report in Table 3 part of the results for 15% censoring, τ = 0.25, n = 500 and three selected values of c.

Table 3.

Results of PROM methods at different bandwidth values hn = cn–0.15 in Example 1 with 15% censoring, τ = 0.25 and n = 500.

Variable Selection 100×MSE
Method c TP FP OP β 0 β 1 β 2 β 3
PROMM 0.2 4.00 0.10 0.91 0.29 0.72 0.73 0.80
1.0 4.00 0.08 0.93 0.26 0.75 0.71 0.81
2.0 4.00 0.08 0.93 0.25 0.73 0.70 0.80

PROMS 0.2 4.00 0.09 0.92 0.29 0.73 0.72 0.79
1.0 4.00 0.07 0.94 0.27 0.73 0.71 0.79
2.0 4.00 0.08 0.93 0.26 0.73 0.71 0.79

As for the methods to obtain the EDR subspace and the corresponding indices in PROM, Tables 1 and 2 show that the SIR and MAVE estimators perform com parably. Since SIR is computationally more convenient, we use SIR to estimate the dimension q and the dimension reduction directions γi in Example 2.

Example 2. In this study, we increase the number of covariates to 100, and generate data from the following model

Ti=1+xi1+xi2+xi3+xi4+(1−0.5xi4)∊i,

where i = 1, · · · , n, εi ~ N(0, 1), and xij ~ U(–1, 1), j = 1, · · · , 100. The true quantile coefficients equal (0.326, 1, 1, 1, 1.337, 0, · · · , 0) at τ = 0.25 and (1, 1, 1, 1, 1, 0, · · · , 0) at τ = 0.50. The observed responses are obtained by Yi = min(Ti, Ci), where Ci ~ U(–2, 18) yielding averaged 15% censoring, and Ci ~ U(–2, 8) yielding averaged 30% censoring. For the developed PROM method, we use SIR to estimate the indices, as well as the dimension q selected by the chi-square test (Li et al., 1999). We report the results of the proposed method with both the true q = 2 and the estimated q̂ by SIR in Tables 4 and 5.

Table 5.

Mean squared errors (×100) of the estimates for nonzero quantile coefficients in Example 2. The methods PROMS (q) and PROMS (q̂) are based on the true and the estimated number of indices, respectively.

τ = 0.25
τ = 0.50
β0 (τ) β1 (τ) β2 (τ) β3 (τ) β4 (τ) β0 (τ) β1 (τ) β2 (τ) β3 (τ) β4 (τ)
n = 200, 15% censoring
OracleS 1.14 3.02 2.67 2.42 3.42 1.15 2.44 2.35 2.25 2.60
PIPW 2.86 4.25 4.48 4.44 5.69 1.33 3.43 3.94 3.18 3.33
PROMS (q) 1.20 5.10 6.50 5.03 3.72 1.18 3.17 3.40 2.91 3.43
PROMS (q̂) 1.27 5.24 6.03 5.36 4.11 1.17 3.20 3.45 3.10 3.50

n = 200, 30% censoring
OracleS 1.24 3.44 3.21 2.97 3.59 1.53 3.04 2.71 2.67 3.13
PIPW 6.42 11.13 12.24 9.89 11.58 3.03 8.93 8.77 7.73 7.97
PROMS (q) 1.65 8.74 10.59 10.46 7.17 1.61 6.36 7.42 7.05 5.81
PROMS (q̂) 1.61 7.18 12.14 9.60 8.67 1.59 5.90 9.32 6.83 4.91

n = 500, 15% censoring
OracleS 0.45 1.06 1.07 1.11 1.22 0.38 1.02 0.97 0.88 1.04
PIPW 0.58 1.18 1.35 1.36 1.36 0.37 1.30 1.14 1.07 1.28
PROMS (q) 0.45 1.03 1.09 1.11 1.23 0.40 1.01 0.93 0.89 1.02
PROMS (q̂) 0.46 1.02 1.10 1.11 1.25 0.37 1.04 0.99 0.89 1.04

n = 500, 30% censoring
OracleS 0.51 1.21 1.31 1.29 1.38 0.67 1.15 1.04 1.06 1.34
PIPW 0.67 1.47 1.56 1.38 1.69 0.49 1.59 1.46 1.29 1.55
PROMS (q) 0.52 1.22 1.28 1.29 1.36 0.70 1.16 1.07 1.06 1.36
PROMS (q̂) 0.50 1.25 1.28 1.30 1.36 0.51 1.17 1.09 1.08 1.33

Table 4 suggests that the performance of PIPW in variable selection deteriorates with higher dimension of covariates, and PIPW selects many irrelevant variables in all scenarios. In contrast, the PROM methods have much higher accuracy in variable selection, and their oracle proportions approach to 1 as n increases. For those nonzero quantile coefficients, the PROM estimates have mean squared errors smaller than those of PIPW in most cases. For larger sample size n = 500, the PROM estimates are almost as efficient as the oracle estimates. In addition, both Tables 4 and 5 show that the PROM method based on the estimated q̂ behaves very similarly as that based on the true q.

4. Real Data Analysis

We apply the proposed variable selection procedure to a head and neck cancer clinical trial, conducted by Eastern Cooperative Oncology Group and the Southwest Oncology Group (Adelstein et al., 2003). In addition to the evaluation of the overall effectiveness of three treatments, namely, standard radiotherapy (treat ment A), radiotherapy plus simultaneous Cisplatin (treatment B), and split-course radiotherapy plus simultaneous cisplatin and 5-fluorouracil (treatment C), it was of substantial interest to detect the treatment effectiveness over high-risk patient populations (often characterized by the lower quantiles of the survival distribu tion). Such investigations would potentially lead to more effective next generation therapies for targeted subpopulations.

After excluding the ineligible patients and those with missing data, the data set has 171 subjects, of whom 129 died during the follow-up period. We apply the proposed variable selection method PROM to the data set to study the impacts of the predictors on the τth quantile of the overall survival times (in months). We focus on two quantiles τ = 0.25 and 0.5. The standard radiotherapy treatment (treatment A) is treated as the baseline.

Besides the three treatment arms, we also consider the following continuous or ordinal confounders, including age, height, weight, weight loss, tumor di er entiation, size of the primary tumor (cm), and categorical variables including gender (1 = female, 0 = male), race (1 = white, 0 = black), smoking (non- smoker, light cigarette smoker with less than 20 packs a year, moderate cigarette smoker with 20-40 packs a year, heavy cigarette smoker with more than 40 packs a year), alcohol drinking (light drinker: consuming less than 10oz whiskey a week or equivalent, moderate drinker: consuming 10-32 oz whiskey a week or equivalent, heavy drinker: consuming more than 32 whiskey a week or equivalent) and primary tumor site (oralcavity, orapharynx, hypopharynx, larynx with oralcavity). The continuous and ordinal variables are standardized to have mean zero and standard deviation one. For the categorical variables smoking, alcohol drinking and primary tumor site, we treat nonsmoker, light drinker and oralcavity as the baseline. Therefore, the full model contains total p = 19 coefficients including the intercept effect, which presents the τth quantile of survival times for a male, black, nonsmoking and light-alcohol-drinking patient who received standard radiotherapy treatment and had average age, height, weight, weight loss, average tumor size and moderately well differentiated oral cavity tumor.

By using the selection criterion in Xia et al. (2010), the dimension of the central subspace (CS) is selected to be 4 for this data set. The 4 indices are then estimated by the MAVE method of Xia et al. (2010) and the SIR method of Li et al. (1999). Results based on the MAVE indices estimation are very similar to those based on SIR and thus are omitted. The sparse index-based estimation of the coefficients is obtained for τ = 0.25 and 0.5 as in Sections 2.3 and 2.4 with the bandwidth hn selected by 5-fold cross validation and the tuning parameter λn selected by minimizing SBIC(λn). We use the eighth-order kernel function (Hu and Fan, 1992): K(x) = 1/13(1 – x2)(35 – 385x2 + 1001x4 – 715x6)I(|x| ≤ 1).

The penalized coefficient estimates from the PIPW method (Shows et al., 2010) and the proposed PROM method are summarized in Table 6. Previous analysis (without accounting for any confounders) pointed that, in terms of effect, treatment A di ered from C significantly, while only di erring from B marginally (Adelstein et al., 2003). Our quantile regression analysis reveals some interesting heterogeneity in the population. The treatment B tends to be more effective for more severe cases (at the lower quartile), while both treatments B and C show positive effects for the typical cases (at the median). In addition, height shows negative effects at the median, while white patients tend to have longer median survival than black patients. These two effects are not selected for more severe cases. On the other hand, larger size of primary tumor is associated with shorter survivals for more severe cases, but not at the median. The two tumor sites hypopharynx and larynx show positive effects at both quantiles, while oropharynx tends to have no effect. As suggested in the simulation study, PIPW had difficulty identifying the correct model for data sets with larger number of predictors. For this data set, PIPW yields more shrinkage, leading to much more sparse models at both quantiles. More specifically, PIPW suggests that both treatments B and C have no di erence than the baseline treatment at both quantiles, and only tumor size, weight and gender have effects on the lower quartile of the survival distribution.

Table 6.

Spare coefficient estimates for the head and neck cancer data set.

τ = 0.25
τ = 0.5
Variable PIPW PROM PIPW PROM
(Intercept) 7.464 10.056 23.0092 35.296
Treatment B 0.000 2.371 0.000 4.806
Treatment C 0.000 0.000 0.000 5.556
Age 0.000 0.000 0.000 0.000
Tumor differentiation 0.000 0.000 0.000 0.000
Weight loss 0.000 0.000 0.206 0.000
White 0.000 0.000 0.000 2.789
Height 0.000 0.000 -4.222 -2.979
Weight 1.370 0.000 8.898 3.882
Gender -0.088 0.000 -8.108 -5.972
Tumor size -1.286 -0.641 0.000 0.000
Hypopharynx 0.000 1.281 0.000 6.756
Larynx 0.000 3.588 5.032 2.261
Oropharynx 0.000 0.000 0.000 0.000
Light smoker 0.000 0.000 0.000 -14.506
Moderate smoker 0.000 -0.542 0.000 -15.644
Heavy smoker 0.000 -0.268 -4.022 -22.128
Moderate drinker 0.000 -3.139 -4.701 -6.056
Heavy drinker 0.000 -3.531 -5.391 -5.561

To further compare the results from PIPW and PROM, we evaluate the risk prediction accuracy of the models selected by the two methods. For each subject i, we define the risk score as the estimated conditional quantile, m(xi)=xiTβ^(τ), where xi is the covariate vector and β^(τ) is the penalized coefficient estimate from either PIPW or PROM at quantile level τ. In our context, good risk scores are expected to better discriminate among subjects with longer and shorter survivals. In medical studies, the concordance measure Ct0 = P{m(x2) > m(x1)|T2 > T1, T1 < t0} is commonly used to evaluate the overall performance of a risk scoring system, where t0 is a prespecified follow-up time point. To account for the censoring, we employ the C-statistic of Pencina and D'Agostino (2004):

C^t0=∑i≠jδiI(Yi<Yj,Yi<t0)I{m(xi)<m(xj)}∑i≠jδiI(Yi<Yj,Yi<t0).

The Ct0-statistics with different values of t0 are summarized in Table 7. Results show that the model selected by PROM has higher risk prediction accuracy than that selected by PIPW at both quantile levels for all values of t0 examined.

Table 7.

Ct0 statistics of the models selected by PIPW and PROM for the head and neck cancer data set.

τ = 0.25
τ = 0.5
t 0 PIPW PROM PIPW PROM
20 0.580 0.613 0.612 0.647
40 0.566 0.602 0.600 0.638
60 0.565 0.599 0.599 0.637
80 0.565 0.598 0.598 0.637

5. Discussion

We have developed a new variable selection approach for censored quantile regression. Our models depict more completely the survival distribution of interest, identifying important factors leading to poor prognosis in survival. This is of particular interest to physicians, who are keen on designing effective treatments for targeted patient sub-population, often characterized by short survival. In contrast, such a small group may be overlooked by using the more popularly used Cox and AFT models. As opposed to the existing methods for censored quantile regression, our developed methods require fewer stringent assumptions, and enjoys computational readiness.

We proposed to estimate the censoring probabilities nonparametrically based on effective dimension reduction indices. The proposed method can be extended to situations where the number of predictors grow in the sample size. However, the theoretical justification requires studying properties of the dimension reduction methods for censored regression data with growing dimensions, and this is beyond the scope of the current paper. To avoid the curse of dimensionality, an alternative way is to estimate the censoring probabilities by fitting a semiparametric regression model (e.g. Cox proportional hazards model). McKeague et al. (2001) showed that this method works reasonably well even when the Cox model is slightly misspecified. However, this approach requires the semiparametric and the quantile regression models to be compatible. We leave the formal investigation of this semiparametric approach to a future study.

Acknowledgement

Wang's research is supported by National Science Foundation (DMS-1007420). Zhou's research is supported by National Science Foundation (DMS-0906665). Li's research is supported by National Institutes of Health (R01CA95747). The authors would like to thank Professor Robert Gray and the Eastern Cooperative Oncology Group for providing the head and neck cancer data set, and Professor Yingcun Xia for providing the matlab program of the MAVE method. We would like to thank two anonymous reviewers, the associate editor and the editor for constructive comments and helpful suggestions.

APPENDIX

To simplify the presentation, we will omit τ in various expressions such as β(τ), ei(τ) and A(τ), and we focus on the cases with q = 1. Proof for q > 1 follows the same line by using the asymptotic properties of the local Kaplan-Meier estimate F^(⋅∣z) for general cases with q ≥ 1, but the notations are more complicated. Define z^i=xiTγ^, zi=xiTγ and z0i=xiTγ0. To reflect the dependence of the weights wi for redistribution of masses on F and γ, we write wi as wi(F, γ). In addition, we define Mn(β,F,γ)=n−1∑i=1nmi(β,F,γ) as the subgradient of the weighted quantile objective function n–1L(β, w) defined in (2.2), where

mi(β,F,γ)=xi{τ−wi(F,γ)I(Yi≤xiTβ)}=xi[{τ−I(Ci>xiTβ,Ti≤xiTβ)−I(Ci≤xiTβ,Ti≤Ci)}−I(Ci≤xiTβ,Ti≥Ci)I{F(Ci∣zi)<τ}{τ−F(Ci∣zi)1−F(Ci∣zi)}].

Let M(β,F,γ)=Ex{m(β,F,γ)}=E[x{τ−H(xTβ∣x)−R(β,F,γ∣x)}], where H(t∣x)={1−G(t∣x)}F0(t∣x)+∫−∞tF0(u∣x)g(u∣x)du, The proof of Theorem 1 depends on the following Lemma.

Lemma 1 Suppose assumptions A1-A5 hold, then we have for any q ≥ 1

(i)‖F^−F0‖H≐suptsupx∣F^(t∣z^)−F0(t∣x)∣=suptsupx∣F^(t∣z^)−F0(t∣z0)∣=Op({logn∕(nhnp)}1∕2+hnν);(ii)F^(t∣z)−F0(t∣z)=∑j=1nBnj(z)ξ(Yj,δj,t,z)+Op({logn∕(nhnp)}3∕4+hnν)a.s.,

where ξ(Y,δ,t,z)={1−F0(t∣z)}[−∫0min(Y,t)f0(s∣z){1−F0(s∣z)}−2{1−G(s∣z)}−1ds+I(Y≤t,δ=1){1−F0(Y∣z)}−1{1−G(Y∣z)}−1]

Proof. Note that

F^(t∣z^)−F0(t∣z0)={F^(t∣xTγ^)−F0(t∣xTγ^)}+{F0(t∣xTγ^)−F0(t∣xTγ0)}. (6.1)

By extending Theorem 2.1 of Gonzalez-Manteiga and Cadarso-Suarez (1994) to general cases with q ≥ 1, we have

suptsupz∣F^(t∣z)−F0(t∣z)∣=Op({logn∕(nhnp)}1∕2+hnν). (6.2)

Lemma 1(i) thus follows by combining (6.1), (6.2), assumptions A1, A2 and A4. Lemma 1(ii) gives the linear representation of F^(⋅) for general cases with q ≥ 1. Its proof is similar to that of Theorem 2.3 in Gonzalez-Manteiga and Cadarso- Suarez (1994). The main difference is that the bias influence and the variance influence are hν and (nhnq)−1 for a νth order kernel function in the q-dimensional context, in contrast to h2 and (nhn)−1 for a second order Kernel function in the one-dimensional context.

Proof of Theorem 1

The consistency of β~ can be easily shown by using Lemma 1(i) and the similar arguments as in the proof Theorem 1 in Wang and Wang (2009). Therefore, we omit the details.

To establish the asymptotic normality, we first prove the √n-consistency of β~ to β0. Define

Γ1=∂M(β,F0,γ0)∂β∣β=β0=−E[xxT{1−G(xTβ0∣x)f0(xTβ0∣x)}].

Note that Γ1 is continuous at β = β0, and has a full rank under assumption A3. Therefore, there exists a constant K such that ‖β~−β0‖≤K‖M(β~,F0,γ0)‖ with probability tending to one. It then suffices to show that ‖M(β~,F0,γ0)‖=Op(n−1∕2).

Using similar arguments as in the proof of Lemma 2 in Wang and Wang (2009), we have

‖Mn(β,F,γ)−Mn(β0,F0,γ0)−M(β,F,γ)‖=op(n−1∕2), (6.3)

uniformly over (β,F,γ) such that ‖β−β0‖≤an,‖F−F0‖H≤an and ‖γ−γ0‖≤an, where an → 0 as n → ∞. Therefore, by the consistency of β~, Lemma 1 and assumption A4,

‖M(β~,F^,γ^)−Mn(β~,F^,γ^)+Mn(β0,F0,γ0)‖=op(n−1∕2). (6.4)

Under assumptions A1, A2 and A4, ‖Mn(β~,F^,γ^)−Mn(β~,F^,γ0)‖=op(n−1∕2). In addition, combining the subgradient condition (Koenker, 2005) and assumption A1 gives ‖Mn(β~,F^,γ^)‖=Op(n−1). Therefore,

‖Mn(β~,F^,γ0)‖≤‖Mn(β~,F^,γ^)−Mn(β~,F^,γ0)‖+‖Mn(β~,F^,γ^)‖=Op(n−1∕2). (6.5)

Let ε > 0 and Fε(t|z) = F0(t|z)+ε{F(t|z)–F0(t|z)}. Following some routine algebra, we can derive the functional derivative of M(β0,F,γ0) at F0 in the direction [F – F0] as

Γ2(β0,F0,γ0)[F−F0]=lim∊→01∊[M{β0,F0+∊(F−F0)},γ0)−M(β0,F0,γ0)]=(1−τ)E[x∫−∞xTβ0F(t∣z0)−F0(t∣z0)1−F0(t∣z0)g(t∣x)dt].

Therefore,

Γ2(β0,F0,γ0)[F^−F0]−(1−τ)E[x∫−∞xTβ0F^(t∣z0)−F0(t∣z0)1−F0(t∣z0)g(t∣x)dt].

By plugging in the linear representation of F^(t∣z)−F0(t∣z) in Lemma 1(ii) and by applying the Taylor expansion, we obtain

Γ2(β0,F0,γ0)[F^−F0]=n−1∑i=1n(1−τ)ϕi+op(n−1∕2), (6.6)

where ϕi=xi∫−∞xiTβ0g(t∣z0i)ξ(Yi,δi,z0i){1−F0(t∣xi)}−1dt, and ξ(Y,δ,t,z) is defined in Lemma 1(ii). Here ϕi are independent random variables with mean zero, and Γ2(β0,F0,γ0)[F^−F0]=Op(n−1∕2).

Let δn = o(1) be a positive sequence such that P(‖β~−β0∣≥δn,‖F^−F0‖H≥δn)→0. Following the routine Taylor expansion, we can show that under assumptions A1-A2, ‖Γ2(β~,F0,γ0)[F^−F0]−Γ2(β0,F0,γ0)[F^−F0]‖=‖β~−β0‖op(1), and

‖M(β~,F^,γ0)−M(β~,F0,γ0)−Γ2(β~,F0,γ0)[F^−F0]‖≤K‖F^−F0‖H2 (6.7)

for a constant K ≥ 0. By Lemma 1(i), for any α > 0 such that 1∕(4ν)<α<1∕(2q),‖F^−F0‖H2=op(n−1∕2). This together with (6.6) gives

‖M(β~,F0,γ0)‖≤‖M(β~,F0,γ0)−M(β~,F^0,γ0)‖+‖M(β~,F^0,γ0)≤‖M(β~,F^0,γ0)−M(β~,F0,γ0)−Γ2(β~,F0,γ0)[F^−F0]‖+‖M(β~,F^,γ0)+‖Γ2(β~,F0,γ0)[F^−F0]−Γ2(β0,F0,γ0)[F^−F0]‖+‖Γ2(β0,F0,γ0)[F^−F0]‖≤‖F^−F0‖H2+‖β~−β0‖op(1)+Op(n−1∕2)≤‖M(β~,F0,γ0)‖op(1)+Op(n−1∕2). (6.8)

Therefore, combining (6.4)-(6.8) gives β~−β0≤K‖M(β~,F0,γ0)‖=Op(n−1∕2).

Recall from (6.5) that Mn(β~,F^,γ0)=Op(n−1∕2), The rest of the proof of normality is similar to those of Theorem 3.3 in Pakes and Pollard (1989) and Theorem 2 of Chen et al. (2003). Therefore, we just sketch the main steps here. Let Γ3=∂M(β0,F0,γ)∕∂γ∣γ=γ0. By Assumption 4(i) and Taylor expansion, we get M(β~,F^,γ^)−M(β~,F^,β0)=Γ3(γ^−γ0)+op(n−1∕2). define Ln(β)=Mn(β0,F0,γ0)+Γ1(β−β0)+Γ2(β0,F0,γ0)[F^−F0]+Γ3(γ^−γ0). By the root-n consistency result shown above, Lemma 1, (6.3) and (6.7), we have

‖Mn(β~,F^,γ^)−Ln(β~)‖=‖Mn(β0,F0,γ0)+M(β~,F^,γ^)−Ln(β~)+Mn(β~,F^,γ^)−M(β~,F^,γ^)−Mn(β0,F0,γ0)‖≤‖M(β~,F0,γ0)−Γ1(β~−β0)‖+‖M(β~,F^,γ^)−M(β~,F0,γ0)−Γ2(β0,F0,γ0)[F^−F0]‖+‖M(β~,F^,γ^)−M(β~,F^,γ0)−Γ3(γ^−γ0)‖+‖Mn(β~,F^,γ^)−M(β~,F^,γ^)−Mn(β0,F0,γ0)‖=op(n−1∕2).

Similarly, ‖Mn(β~,F^,γ^)−Ln(β‒‖=op(n−1∕2), where β‒ is the minimizer of Ln(β) satisfying n1∕2(β‒−β0)=−Γ1−1{M0(β0,F0,γ0)+Γ2(β0,F0,γ0)[F^−F0]+Γ3(γ^−γ0)}. Since Γ1 is of full rank under assumption A3, with a bit more work, we get n1∕2(β~−β‒)=op(1). Note that Mn(β0,F0,γ0),Γ2(β,F0,γ0)[F^−F0] and γ^−γ0 are the averages of independent random vectors of means zero. Therefore, applying the central limit theorem gives

n1∕2{Mn(β0,F0,γ0)+Γ2(β0,F0,γ0)[F^−F0]+Γ3(γ^−γ0)}→DN(0,V), (6.9)

where V = Cov(vi) with vi=mi(β0,F0,γ0)+ϕi+di. The asymptotic normality of β~ is thus proven.

Proof of Theorem 2

Let A^n={j:β^j≠0}. We first show that for any j∉A,P(j∈A^n)→0 as n → ∞. Suppose there exists a k∈Ac such that ∣β^k∣≠0. Let β* be a vector constructed by replacing β^k with 0 in β^. For notational simplicity, we write w^i=wi(F^,γ^). Note that ∣ρτ(a)−ρτ(b)∣≤∣a−b∣max{τ,1−τ}<∣a−b∣. Therefore, for large enough n,

LAL(β^,w^)−LAL(β∗,w^)=∑i=1nw^i{ρτ(yi−xiTβ^)−ρτ(yi−xiTβ∗)}+∑i=1nw^i{ρτ(y+∞−xiTβ^)−ρτ(y+∞−xiTβ∗)}+λnvk∣β^k∣≥−∑i=1n∣ρτ(y+∞−xiTβ^)−ρτ(y+∞−xiTβ∗)∣−∑i=1nτ∣xiTβ^−xiTβ∗∣+λnvk∣β^k∣≥−2∑i=1n‖xi‖⋅∣β^k∣+λn∣β~k∣−r∣β^k∣>0, (6.10)

where the last inequality holds as ∑i=1n‖xi‖=Op(n) by Assumption A1, and n−1λn∣β~k∣−r≥nr∕2−1λn→∞. The (6.10) thus contradicts the fact that LAL(β^,w^)≤LAL(β∗,w^).

We next show that for any j∈A, P(j∉A^n)→0. We write bA=(bj,j∈A) for any vector b ∈ Rp, and BAA as the sub-matrix of a p × p matrix B with both row and column indices in A. Recall from the proof of Theorem 1 that

Mn(βA,F,γ)=Mn(β0A,F0,γ0)+Γ1AA(βA−β0A)+Γ2AA(βA0,F0,γ0)[F−F0]+op(n−1∕2) (6.11)

uniformly over βA, F and γ such that ‖βA−β0A‖=O(n−1∕2), ‖F−F0‖H=o(n−1∕4) and γ – γ0 = O(n–1/2). Denote βA−β0A=n−1∕2u and let K be some positive constant. By (6.11), for ||u|| = K, we have

nuTMn(βA,F^,γ^)=nuT{Mn(β0A,F0,γ0)+Γ2AA}+n1∕2uTΓ1AAu+op(n1∕2), (6.12)

where Γ2AA=Γ2AA(βA0,F0,γ0)[F^−F0]. Therefore, with probability tending to one,

−nuTMn(βA,F^,γ^)≥−nuT{Mn(β0A,F0,γ)+Γ2}−n1∕2uTΓ1u+o(n1∕2)≥k0n1∕2 (6.13)

for some positive k0. However, the subgradient condition (see the proof of Theorem 1 in Wang and Wang, 2009) requires that

‖nuTMn(β^A,F^,γ^)‖+λn∑j∈Avj∣τ−I(β^j<0)∣≤Op(maxi‖xi‖). (6.14)

When λn = o(n1/2) and assumption A1 holds, (6.13) and (6.14) suggest that the subgradient condition cannot hold if ‖β^A−β0A‖=Kn−1∕2. Using the monotonicity argument in Jurečková (1977), we can show that the gradient condition also can not hold if ‖β^A−β0A‖>Kn−1∕2. Therefore ‖β^A−β0A‖≤Kn−1∕2 with probability tending to one. The proof of Theorem 2 is thus complete.

References

  1. Adelstein DJ, Li Y, Adams GL, Wagner HJ, Kish JA, Ensley JF, Schuller DE, Forastiere AA. An intergroup phase iii comparison of standard radiation therapy and two schedules of concurrent chemoradiotherapy in patients with unresectable squamous cell head and neck cancer. Journal of Clinical Oncology. 2003;21:92–98. doi: 10.1200/JCO.2003.01.008. [DOI] [PubMed] [Google Scholar]
  2. Bang H, Tsiatis A. Median regression with censored cost data. Biometrics. 2002;58:643–649. doi: 10.1111/j.0006-341x.2002.00643.x. [DOI] [PubMed] [Google Scholar]
  3. Beran R. Nonparametric regression with randomly censored survival data. Tech. rep. University of California; Berkeley: 1981. [Google Scholar]
  4. Efron B. The two-sample problem with censored data. In: 88 Le Cam L, Neyman J, editors. Porcessdings Fifth Berkeley Symposium in Mathematical Statistics. IV. Prentice-Hall; New York: 1967. pp. 831–853. [Google Scholar]
  5. Engler D, Li Y. Survival analysis with high-dimensional covariates: an application in microarray studies. Statistical Applications in Genetics and Molecular Biology. 2009;8(Issue 1) doi: 10.2202/1544-6115.1423. Article 14. [DOI] [PMC free article] [PubMed] [Google Scholar]
  6. Hu TC, Fan J. Bias correction and higher order kernel functions. Statistics and Probability Letters. 1992;13:235–243. [Google Scholar]
  7. Huang J, Ma S, Xie H. Regularized estimation in the accelerated failure time model with high-dimensional covariates. Biometrics. 2006;62:813–820. doi: 10.1111/j.1541-0420.2006.00562.x. [DOI] [PubMed] [Google Scholar]
  8. Huang Y. Quantile calculus and censored regression. Annals of Statistics. 2009;38:1607–1637. doi: 10.1214/09-aos771. [DOI] [PMC free article] [PubMed] [Google Scholar]
  9. Koenker R, Geling O. Reappraising medfly longevity: a quantile regression survival analysis. Journal of the American Statistical Association. 2001;96:458–468. [Google Scholar]
  10. Koenker R, Machado JAF. Goodness of fit and related inference processes for quantile regression. Journal of the American Statistical Association. 1999;94:1296–1310. [Google Scholar]
  11. Leng C. Variable selection and coefficient estimation via regularized rank regression. Statistica Sinica. 2010;20:167–181. [Google Scholar]
  12. Leng C, Tong X. A quantile regression estimator for censored data. 2011. Manuscript.
  13. Li KC, Wang JL, Chen CH. Dimension reduction for censored regression data. Annals of Statistics. 1999;27:1–23. [Google Scholar]
  14. Lindgren A. Quantile regression with censored data using generalized l1 minimization. Computational Statistics and Data Analysis. 1995;23:509–524. [Google Scholar]
  15. McKeague IW, Subramanian S, Sun Y. Median regression and the missing information principle. Journal of Nonparametric Statistics. 2001;13:709–727. [Google Scholar]
  16. Pencina M, D'Agostino R. Overall c as a measure of discrimination in survival analysis: model specific population value and confidence interval estimation. Statistics in Medicine. 2004;23:2109–2123. doi: 10.1002/sim.1802. [DOI] [PubMed] [Google Scholar]
  17. Peng L, Huang Y. Survival analysis with quantile regression models. Journal of the American Statistical Association. 2008;103:637–649. [Google Scholar]
  18. Portnoy S. Censored regression quantiles. Journal of the American Statistical Association. 2003;98:1001–1012. [Google Scholar]
  19. Portnoy S, Lin G. Asymptotics for censored regression quantiles. Journal of Nonparametric Statistics. 2010;22:115–130. [Google Scholar]
  20. Rao CR. Large sample tests of statistical hypotheses concerning several parameters with applications to problems of estimation. Mathematical Proceedings of the Cambridge Philosophical Society. 1948;44:50–57. [Google Scholar]
  21. Schwarz GE. Estimating the dimension of a model. Annals of Statistics. 1948;6:461–464. [Google Scholar]
  22. Serfling R. Approximation Theorems of Mathematical Statistics. John Wiley & Sons; 1980. [Google Scholar]
  23. Shows JH, Lu W, Zhang HH. Sparse estimation and inference for censored median regression. Journal of Statistical Planning and Inference. 2010;140:1903–1917. doi: 10.1016/j.jspi.2010.01.043. [DOI] [PMC free article] [PubMed] [Google Scholar]
  24. Subramanian S. Median regression using nonparametric kernel estimation. Nonparametric Statistics. 2002;14:583–605. [Google Scholar]
  25. Wang H, Wang L. Locally weighted censored quantile regression. Journal of American Statistical Association. 2009;104:1117–1128. [Google Scholar]
  26. Wang S, Nan B, Z. J., Beer DG. Doubly penalized buckleyjames method for survival data with high-dimensional covariates. Biometrics. 2008;64:132–140. doi: 10.1111/j.1541-0420.2007.00877.x. [DOI] [PubMed] [Google Scholar]
  27. Wu Y, Liu Y. Variable selection in quantile regression. Statistica Sinica. 2009;19:801–817. [Google Scholar]
  28. Xia Y, Zhang D, Xu J. Dimension reduction and semiparametric estimation of survival models. Journal of the American Statistical Association. 2010;105:278–290. [Google Scholar]
  29. Zhang HH, Lu W. Adaptive lasso for cox's proportional hazards model. Biometrika. 2007;94:691–703. [Google Scholar]
  30. Zou H. The adaptive lasso and its oracle properties. Journal of the American Statistical Association. 2006;101:1418–1429. [Google Scholar]

RESOURCES