Optimal Model Averaging Estimation in

arXiv:1612.08365v3 [stat.ME] 21 Jun 2017

High-dimensional Censored Linear Models
Xiaodong Yan and Xingqiu Zhao∗
Yunnan University and The Hong Kong Polytechnic University

Abstract: This article considers ultrahigh dimensional prediction problems with censored response variables. We propose a two-step model averaging procedure for improving prediction
accuracy of the true conditional mean of a censored response variable. The first step is to
construct a class of candidate models, each with low-dimensional covariates. For this, a feature screening procedure is developed to separate the active and inactive predictors through
a fused mean-variance index and to group covariates with similar size of index together
to form regression models with censored response variables. The new model-free screening
method can easily deal with many types of predictors and response variables, such as discrete, categorical and continuous variables, still works well when predictors have heavy-tailed
distributions or strongly depend on each other, and enjoys rank consistency properties under
mild regularity conditions. The second step is to find the optimal model weights for averaging by adapting a delete-one Mallows criterion, where the standard constraint that weights
sum to one is removed. The theoretical results show that the delete-one Mallows criterion
achieves the lowest possible prediction loss asymptotically. Numerical studies demonstrate
the superior performance of the proposed variable screening and model averaging procedures
over existing methods.
Key words and phrases: Asymptotic optimality, high-dimensional censored regression, feature screening, model averaging.
∗

Correspondence to: Xingqiu Zhao, Department of Applied Mathematics, The Hong Kong Polytechnic

University, Hong Kong, China (E-mail: xingqiu.zhao@polyu.edu.hk)

1. INTRODUCTION
Improving prediction accuracy for the true conditional mean of censored response variable with ultrahigh-dimensional predictors is a challenging problem in many scientific fields
such as genomics, medicine, economics, finance, and public health. Examples include gene
expression data such as the mantle cell lymphoma data that motivated this research. It has
long been recognized that model selection and model averaging are two popular methods for
enhancing prediction accuracy in regression analysis. However, model selection procedures
such as Akaike information criterion (AIC) and Bayesian information criterion (BIC) offer
one best model for prediction, which may ignore the additional uncertainty or even lead
to biased prediction (Hjort and Claeskens 2003). The model averaging approach compromises across a set of competing models by taking into account the model uncertainty and
optimizing model weights to minimize prediction errors.
For the analysis of high-dimensional data, many penalized methods, such as the least
absolute shrinkage and selection operator (Lasso; Tibshirani 1996), the smoothly clipped
absolute deviation (SCAD; Fan and Li 2001), the adaptive Lasso (Zou 2006), the Dantzig
selector (Candes and Tao 2007), and the minimax concave penalty (MCP; Zhang 2010)
have been developed to simultaneously select the important covariates and estimate parameters in various statistical models when the number of covariates diverges. High-dimensional
sparse modeling with censored survival data is of great practical importance. It is commonly
assumed that only a small number of covariates actually contributes to survival models considered, which leads to the well-known sparse survival models for helping interpretation and
improving prediction accuracy (Bradic, Fan, and Wang 2011). Several regularization methods originally developed for linear regression have been adapted to survival models. For
example, Tibshirani (1997) and Fan and Li (2002) extended the Lasso and nonconcave penalized likelihood, respectively, to the Cox model, while Zhang and Lu (2007) and Zou (2008)
developed the adaptive Lasso and the efficient and adaptive shrinkage methods for variable
selection in the Cox model. Antoniadis, Fryzlewicz, and Letué (2010) studied the Dantzig
selector for the Cox model in a high-dimensional setting. Variable selection techniques have
also been extended to other survival models including the additive hazards model (Leng and
Ma 2007; Martinussen and Scheike 2009; Lin and Lv 2013) and the accelerate failure time
1

model (Huang et al. 2006; Wang et al. 2008).
On the other hand, a number of model averaging procedures are developed for uncensored
response variables under the standard setting in which the number of predictors is fixed and
much smaller than the sample size. These methods include the forecasting model averaging
(Newbold and Granger 1974), the least squares model averaging (Hansen 2007; Wan, Zhang,
and Zou 2010), the predictive likelihood-based model averaging (Ando and Tsay 2010), the
frequentist model averaging (Liang et al. 2011), the jackknife model averaging (Hansen and
Racine 2012), the heteroskedasticity-robust Cp model averaging (Liu and Okui 2013), the
optimal model averaging for linear mixed-effects models and generalized linear mixed-effects
models (Zhang, Zou, and Liang 2014; Zhang et al. 2015). In contrast, very little research has
been done on the development of model averaging in high-dimensional settings, particularly
when the number of predictors is much larger than the sample size. For example, Ando
and Li (2014) developed a novel model averaging procedure for high-dimensional regression
models using a delete-one cross-validation criterion. This effective approach is designed for
uncensored response variables, but not for censored response variables, and the formulation is
not applicable to high-dimensional censored regression models. To the best of our knowledge,
there does not seem to exist any established model averaging method for improving prediction
accuracy of censored response variables in high-dimensional sparse regression models.
In this article, we present a two-step model averaging approach for high-dimensional
censored regression. The first step is to prepare a class of candidate models for averaging.
Instead of using the existing variable screening procedures, we propose a new variable screening procedure through the fused mean-variance index to separate the active regressors and
the inactive ones and group regressors with similar size of index together to form regression models with censored response variables. The second step is to find the optimal model
weights for averaging a class of pre-constructed censored regression models. Instead of using
a delete-one cross-validation criterion, we propose to use a delete-one Mallows criterion to
achieve this goal.
The main contributions of this article are threefold. First, we develop a fully nonparametric variable screening method based on a fused mean-variance index. This method can
easily deal with many types of predictors and response variables, such as discrete, categorical

2

and continuous variables, still works well when predictors have heavy-tailed distributions or
strongly depend on each other, and enjoys rank consistency properties under mild regularity
conditions. The superior performance of the new method over the existing variable screening
methods is demonstrated by simulations. Second, the adapted delete-one Mallows criterion
to censored situations can asymptotically minimize the squared error between the true mean
and the predicted value, where the standard constraint that the weights sum to one is removed. Third, simulation results show that the combination of the new variable screening
method for preparing a class of candidate models and the delete-one Mallows criterion for averaging yields a superior method over the existing methods, including AIC model averaging,
BIC model averaging, Lasso, Group Lasso, SCAD, and MCP.
The reminder of the article is organized as follows. We introduce models and present a
two-step model averaging procedure in Section 2, and establish the theoretical properties of
the proposed model averaging method in Section 3. The new variable screening procedure
through a fused mean-variance index and its corresponding rank consistency property are
summarized in Section 4. Simulation and real data analysis results are reported in Sections 5
and 6, respectively, to assess the performance of the proposed procedure compared to several
existing methods. Some concluding remarks are made in Section 7. All technical details are
provided in the Appendix and Supplementary Materials.

2. MODELS AND ESTIMATION PROCEDURE
Consider a survival study with n subjects and p covariates denoted by X = (X1 , · · · , Xp )>,
where p > n. Let T and C denote the survival and censoring times, respectively. The
observations consist of {Ti∗ , δi , Xi , i = 1, . . . , n}, independent copies of {T ∗ , δ, X} with T ∗ =
T ∧ C and δ = I(T ≤ C). Let V be a known monotone transformation of T and U be
the corresponding transformation of C. The goal is to predict the true mean of censored
response variable V based on the observed data {Yi , δi , Xi , i = 1, . . . , n}, independent copies
of {Y, δ, X} with Y = V ∧ U . Suppose that V takes the following semiparametric linear
regression model:
Vi = X>
i β + i , i = 1, . . . , n,

(1)

where β = (β1 , . . . , βp )> is a vector of unknown regression coefficients and i ’s are independent

3

and identically distributed random errors with mean 0 and variance σ2 . Here, we have
dropped the intercept term by assuming that the means of the response variable and the
input variables are already subtracted out. Assume that i ’s are independent of (Xi , Ui )’s.
To improve prediction accuracy in high-dimensional censored regression, we propose a
two-step model averaging procedure as follows.
Step 1: Construct the Candidate Models.
We develop a new feature screening method which will be presented in Section 4 and use
this method to partition the p predictors into K + 1 groups, where the first group has the
highest values of the marginal utility and the (K + 1)th group has values closest to zero.
The (K + 1)th group is removed. Let Ak be the index set of predictors to be included in the
kth group. Then we can construct K candidate models M1 , · · · , MK by
Mk : Vi =

X

Xij βj + i , i = 1, . . . , n.

(2)

j∈Ak

For estimation of βj , we use the weighted least squares method described below.
Step 2: Determine the Optimal Model Weights for Averaging.
To achieve the optimality of model weights, we propose to use a delete-one Mallows
criterion without the classical restriction that weights sum to one. The details are given
below.

2.1 Weighted Least Squares Estimation for Model Mk
Now we rewrite Mk in the matrix form
V = X k β k + ,
where β k is a pk -dimensional vector of unknown parameters, and  is a n-dimensional noise
vector.
For the estimation of β k , we use the weighted least squares (WLS) method proposed by
Zhou (1992). The weights, considering censoring information, are given by
πi = π(Yi ) =

δi
U
Ḡ (Y

i −)

,

(3)

where GU (t) = P (U ≤ t), and ḠU (t) = 1 − GU (t). Let GU (t−) denote the left limit of the
4

function GU (·) at t. For the largest observation of Yi , we adjust its weight by
π(n) =

1
ḠU (Y

(n) −)

.

(4)

In practice, we approximate GU (t) by the Kaplan-Meier estimator GUn (t), and replace
πi with πni , accordingly. Based on the weights in (3) and (4), the WLS estimator of β k is
defined by
b = arg min(Y − X k β )>π n (Y − X k β ) = arg min ||Y ∗ − X ∗ β ||2 ,
β
k
k
k
n
nk k

(5)

where π n is a n-dimensional diagonal matrix with the ith diagonal element being πni , Y ∗n is
√
a n-dimensional vector with the ith component being πni Yi , and X ∗nk is a matrix with the
√
ith row vector being πni X ik , for i = 1, . . . n. Therefore, the WLS estimator of coefficients
in kth model is given by
b = (X ∗> X ∗ )−1 X ∗> Y ∗ = (X ∗> X ∗ )−1 X>Ỹ n ,
β
k
nk
nk
nk n
nk
nk
k

(6)

where Ỹ n is a n-dimensional vector whose components is πni Yi . So the true conditional
b for k = 1, . . . , K.
b nk = X k β
mean of V can be estimated by u
k

2.2 Optimal Model Weight Selection for Averaging
After applying the WLS method to each candidate model, we achieve a list of predictors
b nK }. Now we determine the weight of each candidate model. The kth synthetic
{b
un1 , · · · , u
∗ −1 >
hat matrix X k (X ∗>
nk X nk ) X k π n of each candidate model is denoted by Hnk .

Let ω = (ω1 , · · · , ωK ) be the weight vector of the K models and
V = {ω ∈ [0, 1]K : 0 ≤ ωk ≤ 1}.
Here, the standard restriction

K
P

ωk = 1 is removed. The advantages without this restriction

k=1

are as given in Ando and Li (2014).
We denote the true conditional mean of V on X by u. Then its model average predictor
is given by
b n (ω) =
u
=

K
P
k=1
K
P

b nk =
ωk u

K
P

∗ −1 >
ωk X k (X ∗>
nk X nk ) X k π n Y

k=1

ωk Hnk Y = Hn (ω)Y ,

k=1

5

K
P

where Hn (ω) =

ωk Hnk is the corresponding weighted hat matrix.

k=1

Ando and Li (2014) used the delete-one cross-validation to estimate the model weights.
To consider the influence of variance σ2 on weights, we propose to use the delete-one Mallows
criterion (Wan, Zhang, and Zou 2010) to determine the weights. We first introduce more no(−m)

tation. Let ũnk

be the predicted value of the mth observation from the kth model Mk , which
(−1)

(−n)

is derived from the observations except for (Ym , δm , Xm ). Define ũnk = (ũnk , . . . , ũnk )>.
As shown in Li (1986), we have ũnk = H̃nk Y , where H̃nk is the smoothing matrix given by
H̃nk = Dnk (Hnk − I) + I and Dnk is the n × n diagonal matrix with jth diagonal element
equal to (1−hnkj )−1 , where hnkj is the jth diagonal element of Hnk . The delete-one predictor
is given by
ũn (ω) =

K
X
k=1

K
P

where H̃n (ω) =

ωk ũnk =

K
X

ωk H̃nk Y = H̃n (ω)Y ,

k=1

ωk H̃nk . Thus, the delete-one Mallows criterion can be written as

k=1

M(ω) = (Y − ũn (ω))>π n (Y − ũn (ω))
= (Y − H̃n (ω)Y )>π n (Y − H̃n (ω)Y ).
Then, minimizing the M(ω) over the space V yields the selected weights
b = arg min M(ω).
ω
ω∈V

(7)

3. THEORETICAL RESULTS OF MODEL AVERAGING
K
P
k=1

b n (ω) =
Let X = (X1 , · · · , Xn ) and u = E(V |X). Define a model averaging estimator u
b n (ω))>π n (u − u
b n (ω)). The correωk Hnk Y , and consider loss function Ln (ω) = (u − u

sponding risk function is
R(ω) = E[Ln (ω)|X].

(8)

Let ξn = inf R(ω) Let E(A) be the eigenvalue of matrix A, and ϑ(A) be the maximum
ω∈V

diagonal element of matrix A. In the following, B and Bj ’s are finite constants. To establish
the property of proposed model averaging approach with the adapted delete-one Mallows
criterion, we need the following regularity conditions.
6

(C1) Let bV = sup{t, P (V ≤ t) < 1}, bU = sup{t, P (U ≤ t) < 1}, and τ0 denote the
maximum follow-up time. The first derivative of ḠU (t) is uniformly bounded, bV ≥ bU ,
and for any t ∈ [0, τ0 ], P (t ≤ V ≤ U ) ≥ λ0 for some positive constant λ0 .
k)
(C2) (i) sup ϑ(H
≤
pk

(iv)

k
||u||2
n

B1
;
n

k
(ii) sup np3/4
≤ B2 ; (iii) For some fixed integer κ,

k

K 4κ+2 ||u||2κ
2κ
ξn

= op (1);

≤ B3 .

(C3) (i) E(X>k X k ) is nonsingular, and supk E(X>k X k ) = O(n).
Condition (C2) is the same as conditions required in Theorem 1 of Ando and Li (2014). In
condition (C3)(i), the nonsingular assumption is necessary for hat matrix. The bound of
largest eigenvalue can ensure that the distance between

X>
kXk
nḠU
n (t)

and

X>
kXk
nḠU (t)

can be controlled

only through the distance between 1/ḠU (t) and 1/ḠUn (t). If the order of ξn is n1−γ with
−1

γ ≥ 0, condition (C2)(iii) is reduced to K 1+0.5κ

= op (n(1−2γ)/4 ). Since κ is fixed and the

term 0.5κ−1 is ignorable, K is allowed to grow to infinity if γ < 1/2.
Theorem 1. Under conditions (C1)-(C3), we have
Ln (b
ω)
→1
inf Ln (ω)

(9)

ω∈V

in probability.

4. FUSED MEAN-VARIANCE FILTER
4.1 Screening Method
Instead of using the existing feature screening, we propose a new nonparametric screening
procedure for ultrahigh-dimensional survival data, motivated by the independence feature
screening for ultrahigh dimensional discriminant analysis in Cui, Li, and Zhong (2015). Let
Fj (x|T ) = P (Xj ≤ x|T ) denoting the conditional distribution function of Xj given T , and
Fj (x) = P (Xj ≤ x) being the unconditional distribution Fj (x) of Xj . Note that Xj is
independent of T if and only if Fj (x|T ) = Fj (x). To measure the dependence between Xj
and T , following Cui, Li, and Zhong (2015), we consider a mean-variance based index for
the jth covariate as
ZZ
M Vj = E{Var(F (Xj |T ))} =

{Fj (x|T = t) − Fj (x)}2 dFj (x)dFT (t).
7

(10)

Our marginal utility M Vj also has a remarkable property that M Vj = 0 if and only if Xj
and T are statistically independent, which motivates us to utilize it for feature screening to
characterize both linear and nonlinear relationships in ultrahigh dimensional data analysis.
In order to use M Vj , we can slice the response to obtain an approximation of M Vj . Define
a partition:
n
o
S = [ag , ag+1 ) : ag < ag+1 , g = 1, · · · , S .

(11)

Here we set a1 = 0, aS+1 = bF , where bF = sup{t, F (t) < 1} with F (t) = P (T ≤ t). Here,
[ag , ag+1 ) is called a slice. We also define a random variable G ∈ {1, · · · , S} such that G = g
if and only if T is in the gth slice. The multiclass problem considered by Cui, Li, and Zhong
(2015) can be regarded as a particular case, that is, T = G taking values from 1, · · · , S.
Define
M VjS

=

S
X

Z
pSg

{FjS (x|G = g) − Fj (x)}2 dFj (x),

(12)

g=1

where pSg = P (G = g), FjS (x|G = g) =

P (Xj ≤x,G=g)
.
pSg

The form of (12) implies that M VjS can be

denoted as the weighted average of Cramér-von Mises between the conditional distribution
of X given the slice G = g and the unconditional distribution function of X.
Now we provide the sample version of M VjS . Write
n

pbSg

n

1X
1X
=
I{ag ≤ Ti∗ < ag+1 }πni , Fbj (x) =
I{Xi ≤ x},
n i=1
n i=1

(13)

n

1X
I{Xi ≤ x, ag ≤ Ti∗ < ag+1 }πni /b
FbjS (x|G = g) =
pSg ,
n i=1

(14)

where πni ’s are as defined in Section 2. For censored survival response, it is obvious to use
its sample counterpart to estimate M VjS as follows:
n

S

S
1 X X S bS
d
pbg {Fj (xi |G = g) − Fbj (xi )}2 .
M
Vj =
n i=1 g=1

(15)

Furthermore, motivated by Cook and Zhang (2014) and Mai and Zou (2015), we use
the idea of fusion to improve the efficiency of the mean-variance based measure. Taking L
different slice schemes, where each slice scheme Sl contains Sl intervals, l = 1, · · · , L, we

8

propose the fused mean-variance filter given by
F\
MV j =

L
X

Sl
d
M
Vj

(16)

M VjSl .

(17)

l=1

as an estimate of
F M Vj =

L
X
l=1

Mai and Zou (2015) suggested that Sl ≤ dlog ne for all 1 ≤ l ≤ L so that there is a decent
sample size within each slice for all slicing schemes, where dxe denotes the integer part of x.
In practice, by a technical condition, we take Sl = 3, . . . , dn1/3 e with each partition Sl .
To evaluate F\
M V j , we still need to determine Sl . If the distribution of T is known, then
we can consider an oracle uniform slicing to form partitions Sl through Sl intervals with
ag = F −1 ((g − 1)/Sl ), g = 1, . . . , Sl . In practice, F is unknown and can be estimated by the
Kaplan-Meier estimator Fbn (t). So, we can estimate ag by b
ag = Fb−1 ((g − 1)/Sl ). Write
n

b
Sl = {[âg , âg+1 ) : âg < âg+1 , g = 1, . . . , Sl , â0 = 0, âSl +1 = bF }
as an intuitive uniform slicing. For the oracle uniform slicing, set
(◦)
Qj

=

L
X

Sl
j

MV

b(◦) =
and Q
j

L
X

Sl
d
M
Vj

l=1

l=1

and for the intuitive uniform slicing, set
Qj =

L
X

MV

b
Sl
j

bj =
and Q

l=1

L
X

Sl
d
M
Vj .
b

l=1

When p  n, we define the active predictor subset by
D = {j : F (t|Xj ) functionally depends on Xj for some t}.

(18)

bj as a marginal utility to measure the importance of Xj for censored
Then, we choose Q
response variable T , and estimate D by
b = {j : Q
bj is among the dn th largest},
D
9

(19)

where dn is the predefined positive integer. We refer this procedure to the fused MV-based
sure independence screening, FMV for short.

4.2 Theoretical Property
In this section, we establish the rank consistency properties of our screening approach.
Obviously, we aim at partitioning each of the “observed predictors” into signal and nonsignal components. The signal ones contain the set of predictors to be included in the model
for model averaging. Throughout this section, B denotes a generic positive constant.
To establish the ranking consistency of the FMV, we need the following regularity conditions. Define Smin = min{Sl }, and Smax = max{Sl }.
l

l

(C4) There exists a set E such that D ⊂ E and ∆E = minl {min M VjSl − max M VjSl } > 0,
j∈E

−3/4

n

3/4

log(n)

j ∈E
/

c

Smax /∆E = op (1), and Smax = O(n ).

(C5) For any b1 , b2 such that P {t ∈ [b1 , b2 )} ≤

1+∆E
,
Smin

we have sup |Fj (x|t1 )−Fj (x|t2 )| ≤
x∈RX

∆E
8

for all j and t1 , t2 ∈ [b1 , b2 ).
(C6) P (T ≤ C|X, T ) = P (T ≤ C|T ).
Here, conditions (C4) and (C5) are similar to conditions (C1) and (C2) given in Mai and
Zou (2015) with a detailed discussion.
Condition (C6) makes the censored model flexible enough to allow for a dependence
between X and C. In other words, given the time of death, the covariates do not provide
any further information whether censoring will occur or not. So condition 7 also shows the
R R
1
independence of δ and X on T and implies that F(t, x) = y≤t w≤x Ḡ(y−)
F ∗ (dy, dw), where
F(t, x) = P (T ≤ t, X ≤ x), F ∗ (t, x) = P (Y ≤ t, X ≤ x, δ = 1) (Stute 1996).
Theorem 2. (Ranking Consistency Property) Suppose that Conditions (C1), (C4)-(C6)
bj − max Q
bj } > 0
hold. If n2c−1 ∆−2 = o(1) and nc−1 log(Lp)∆−2 = o(1), then lim inf {min Q
E

E

n→∞

j∈D

j6∈D

a.s..
bj values of active predictors are always ranked beyond
The theorem implies that the Q
those of inactive ones with high probability. Thus, we can separate the active and inactive
predictors through taking an ideal threshold value.
10

5. SIMULATION STUDIES
To assess the performance of the proposed two-step model averaging methods, we conducted simulations with two objectives. The first objective is to compare the proposed new
screening procedure (FMV) with some existing ones for constructing candidate models for
averaging; these include marginal correlation screening (SIS; Fan and Lv 2008), nonparametric independence screening (NIS; Fan, Feng, and Song 2011), sure independent ranking
and screening (SIRS; Zhu et al. 2011), quantile-adaptive screening (QAS; He, Wang, and
Hong 2013) and the fused Kolmogorov filter screening (FKS; Mai and Zou 2015). The
second objective is to compare the proposed model averaging method with the classical
weight selection-based model averaging methods and the penalized variable selection methods; these include model averaging with the Akaike information criterion under the restric10
P
ωk = 1 (MAIC), model averaging with the Bayesian information criterion under the
tion
k=1

restriction

10
P

ωk = 1 (MBIC), FMV-based model averaging with the delete-one Mallows

k=1

criterion under the restriction

10
P

ωk = 1 (MCV1), FMV-based model averaging with the

k=1

delete-one Mallows criterion without the restriction

10
P

ωk = 1 (MCV2), the Lasso method

k=1

(Tibshirani 1996), the group lasso procedure (Glasso; Yuan and Lin 2006), the penalized
regression by SCAD approach (Fan and Li 2001), and the penalized regression by the MCP
approach (Zhang 2010).
In this simulation study, we suppose that censored response variable T took the following
accelerated failure time (AFT) model
log(Ti ) = X>
i θ + i ,
where Xi = (xi1 , . . . , xip ) was assumed to be generated from a multivariate normal distribution with zero mean and covariance matrix Φ = (djl ) with djl = a|j−l| , and i was assumed
to follow the standard normal distribution N (0, σ2 ). Furthermore, we generated censoring
time Ci from min(l, Unif(0, l + 2)) where l controls the censoring rate. The true coefficients
θj were generated from the normal with mean of 0 and standard deviation (SD) of 0.5. Set
sample size n = 100 and number of regressors p = 2000. Set the number of true regressors
11

s = 50 and let the true regressors xk be spaced evenly, k = 40(j − 1) + 1, j = 1, . . . , 50.
We took a = 0.5, σ = 0.5, and l = 8 for about 35% censoring proportion averagely. We
took the largest number of slices for our FMV approach to be dn1/3 e based on Condition 5.
Therefore, three slicing schemes are involved with Sl ∈ {3, 4, 5}.
To prepare candidate models for averaging, we ordered the regressors for grouping using
the proposed FMV and the above mentioned screening methods. For the QAS, we considered
two quantiles 0.2 and 0.4, denoted by QAS1 and QAS2. Since these screening methods
except FMV are designed for uncensored data, we performed the SIRS and FKS for response
variable T ∗ = min(T, C) and the SIS, QAS1, QAS2 and NIS for response variable log(T ∗ ),
respectively. We separated the first 100 active predictors by each screening method, and set
K = 10 to yield a class of 10 candidate models, each with 10 regressors. Then, we conducted
the proposed model averaging procedure presented in Section 2 for the obtained candidate
models from each screening approach.
We used weighted MSE

1
(u
n

b n (b
b n (b
−u
ω ))>π n (u − u
ω )) as the performance measure for

each method. Figure 1 shows the boxplots of weighted MSEs after 200 replications. As
shown in this figure, our screening method FMV, combing with the proposed model weight
selection procedure, seems to have a superior performance in the sense that it achieves the
smallest MSE median.
F igure 1 about here
Furthermore, we used R package program ncvreg to implement the MCP and SCAD
algorithms and perform 10-fold cross-validation for these penalized regression models over a
grid of values for the regularization parameter to select an optimal size of penalty. We implemented Lasso through using lars package in R and the BIC for selecting appropriate tuning
parameter. Group Lasso was implemented by grplasso package in R. We also partitioned
the regressors into K + 1 groups. The first K groups are the same as those obtained by
MCV2, MAIC, and MBIC. The last group consists of all the remaining regressors. We also
used MSE as the performance measure for each approach. Figure 2 displays the boxplots
of MSE values for each method based on the same 200 simulated data sets. The simulation results demonstrate that our proposed model averaging procedure has the superior
performance over the classical model averaging and the penalized methods considered.

12

F igure 2 about here

6. APPLICATION
In this section, we applied the proposed method to the mantle cell lymphoma (MCL)
microarray data, which is available from the web site http://llmpp.nih.gov/MCL. These
data have been studied by Rosenwald et al. (2003). The primary goal of this study was to
identify genes that have significant influence on patients’ survival risk. Among 101 untreated
patients with no history of previous lymphoma, 92 were classified as having MCL based on
the morphologic and immunophenotypic criteria. The data include the expression values of
6312 genes for each patient. During the follow-up, 64 patients died of MCL and the other
28 patients were censored. Here, the sample size n = 92, the number of predictors p = 6312
and the censoring rate is 36%.
We applied the proposed FMV independence screening procedure to the censored MCL
n
data and took the first 4d log(n)
e = 80 ranked genes as active predictors. Suppose that the

survival time T follows the AFT model considered in our simulations. Then, we set K = 10
to yield a class of 10 candidate models, each with 8 genes.
To evaluate the prediction performance of various methods, we randomly selected 50
observations for model fitting, and used the rest of the data as the testing set. Let A = {i :
observation i belongs to the testing set}. We used the weighted averaged squared prediction
errors (WASPE)
1

X δi
| log(Ti ) − µ
bi (b
ω )|2
δ
π
i∈A i i∈A ni

ASPE = P

as the prediction performance measure for each method. We calculated the ASPE from the
estimation results by applying the proposed optimal weight selection method to the obtained
candidate models based on the randomly selected subset. We also calculated the ASPEs
by using MCV1, MCV2, MAIC, MBIC, Glasso, SCAD, MCP and LASSO as described in
Section 5. Figure 3 shows the boxplots of ASPE values obtained from these methods based
on 200 replications. The prediction results demonstrate the superior performance of the
proposed two-step model averaging method over the classical weight selection methods and
the penalized variable selection methods for the MCL data.
F igure 3 about here
13

7. CONCLUDING REMARKS
For the analysis of high-dimensional data, many penalized methods such LASSO, SCAD
and MCP, have been developed. These penalized procedures, however, may not perform well
for a very large number of covariates because ultrahigh dimensionality brings simultaneous
challenges of computational expediency, statistical accuracy and algorithmic stability (Fan,
Samworth, and Wu 2009). In this article, to address the new challenges from ultrahighdimensional data in the presence of censoring, we have proposed a two-step model averaging
procedure for improving prediction accuracy. The first step is to construct candidate models
for averaging, while the second step is to find the optimal model weights for predicting. For
the first step, to accommodate censoring in the ultrahigh-dimensional survival data, we have
designed a new nonparametric screening procedure based on a fused mean-variance filter
and established its rank consistency property under very weak regularity conditions. The
proposed screening approach is invariant under the monotone transformation of the response
and can still be powerful when covariates are strongly dependent on each other. For the
second step, a weighted least squares method has been used to estimate regression parameters
for each candidate model subject to censoring, and a delete-one Mallows criterion has been
adapted to determine the optimal model weights, where the standard constraint that weights
sum to one has been successfully removed. The superior performance of the proposed model
averaging approach over the classical model selection methods and the penalized methods
has been demonstrated through simulation studies and real data analysis.
Note that the weighted least squares approach is used for estimation of regression coefficients with each candidate model. To improve the efficiency, we can replace the weighted
least squares estimator with the sieve maximum likelihood estimator (SMLE) presented in
Ding and Nan (2011). The complicated structure of the SMLE would bring new challenges
to the establishment of model weight selection optimality, which would deserve further research. Another direction is to extend the proposed model averaging method to other highdimensional survival models such as the Cox model and the additive hazards model.

14

References
Ando, T. and Li, K. C. (2014). A model-averaging approach for high-dimensional regression. J. Am. Statist. Assoc. 109, 254–265.
Ando, T. and Tsay, R. (2010). Predictive likelihood for the Bayesian model selection
and averaging. Int. J. Forecast. 26, 744–763.
Antoniadis, A., Fryzlewicz, P. and Letué, F. (2010). The Dantzig selector in Cox’s
proportional hazards model. Scand. J. Statist. 37, 531–552.
Bradic, J., Fan, J. and Wang, W. (2011). Penalized composite quasi-likelihood for
ultrahigh-dimensional variable selection. J. R. Statist. Soc. B 73, 325–349.
Candes, E. and Tao, T. (2007). The Dantzig selector: statistical estimation when p is
much larger than n. Ann. Statist. 35, 2313–2351.
Cook, R. D. and Zhang, X. (2014). Fused estimators of the central subspace in sufficient
dimension reduction. J. Am. Statist. Assoc. 109, 815–827.
Cui, H., Li, R. and Zhong, W. (2015). Model-free feature screening for ultrahigh dimensional discriminant analysis. J. Am. Statist. Assoc. 110, 630–641.
Ding, Y. and Nan, B. (2011). A sieve M-theorem for bundled parameters in semiparametric models, with application to the efficient estimation in a linear model for censored
data. Ann. Statist. 39, 3032–3061.
Fan, J., Feng, Y. and Song, R. (2011). Nonparametric independence screening in sparse
ultrahigh-dimensional additive models. J. Am. Statist. Assoc. 106, 544–557.
Fan, J. and Li, R. (2001). Variable selection via nonconcave penalized likelihood and its
oracle properties. J. Am. Statist. Assoc. 96, 1348–1360.
Fan, J. and Li, R. (2002). Variable selection for Cox’s proportional hazards model and
frailty model. Ann. Statist. 30, 74–99.

15

Fan, J. and Lv, J. (2008). Sure independence screening for ultrahigh dimensional feature
space. J. R. Statist. Soc. B 70, 849–911.
Fan, J., Samworth, R. and Wu, Y. (2009). Ultrahigh dimensional feature selection:
beyond the linear model. J. Mach. Learn. Res. 10, 2013–2038.
Hansen, B. E. (2007). Least squares model averaging. Econometrica 75, 1175–1189.
Hansen, B. E. and Racine, J.(2012). Jackknife model averaging. J. Econometrics 167,
38–46.
He, X., Wang, L. and Hong, H. G. (2013). Quantile-adaptive model-free variable screening for high-dimensional heterogeneous data. Ann. Statist. 41, 342–369.
Hjort, N. L. and Claeskens, G.(2003). Frequentist model average estimators. J. Am.
Statist. Assoc. 98, 879–899.
Huang, J., Ma, S. and Xie, H. (2006). Regularized estimation in the accelerated failure
time model with high-dimensional covariates. Biometrics 62, 813–820.
Leng, C. and Ma, S. (2007). Path consistent model selection in additive risk model via
Lasso. Stat. Med. 26, 3753–3770.
Li, K. C. (1986). Asymptotic optimality of CL and generalized cross-validation in ridge
regression with application to spline smoothing. Ann. Statist. 14, 1011–1112.
Liang, H., Zou, G., Wan, A. T. K. and Zhang, X. (2011). Optimal weight choice for
frequentist model average estimators. J. Am. Statist. Assoc. 106, 1053–1066.
Lin, W. and Lv, J. (2013). High-dimensional sparse additive hazards regression. J. Am.
Statist. Assoc. 108, 247–264.
Liu, Q. and Okui, R. (2013). Heteroskedasticity-robust Cp model averaging. Economet.
J. 16, 463–472.
Mai, Q. and Zou, H. (2015). The Fused Kolmogorov Filter: a nonparametric model-free
screening method. Ann. Statist. 43, 1471–1497.
16

Martinussen, T. and Scheike, T. H. (2009). Covariate selection for the semiparametric
additive risk model. Scand. J. Statist. 4, 602–619.
Newbold, P. and Granger, C. W. J. (1974). Experience with forecasting univariate
time series and the combination of forecasts. J. R. Statist. Soc. A 137, 131–149.
Rosenwald, A., Wright, G., Wiestner, A.,Chan, W. C., Connors, J. M., Campo,
E. and Chiorazzi, M. (2003). The proliferation gene expression signature is a quantitative integrator of oncogenic events that predicts survival in mantle cell lymphoma. Cancer
Cell 3, 185–197.
Stute W. (1996). Distributional convergence under random censorship when covariables
are present. Scand. J. Statist. 23, 461-471.
Tibshirani, R. (1996). Regression shrinkage and selection via the Lasso. J. R. Statist. Soc.
B 58, 267–288.
Tibshirani, R. (1997). The lasso method for variable selection in the Cox model. Stat.
Med. 16, 385–395.
Wan, A. T. K., Zhang, X. and Zou, G. (2010). Least squares model averaging by
Mallows criterion. Biometrics 156, 277–283.
Wang, S., Nan, B., Zhu, J. and Beer, D. G. (2008). Doubly penalized Buckley-James
method for survival data with high-dimensional covariates. Biometrics 64, 132–140.
Yuan, M. and Lin, Y. (2006). Model selection and estimation in regression with grouped
variables. J. R. Statist. Soc. B 68, 49–67.
Zhang, C. H. (2010). Nearly unbiased variable selection under minimax concave penalty.
Ann. Statist. 38, 894–942.
Zhang, H.-H. and Lu, W. (2007). Adaptive lasso for Cox’s proportional hazards model.
Biometrika 94, 691–703.
Zhang, X., Zou, G. and Liang, H. (2014). Model averaging and weight choice in linear
mixed-effects models. Biometrika 101, 205–218.
17

Zhang, X., Yu, D., Zou, G. and Liang, H. (2015). Optimal model averaging estimation
for generalized linear models and generalized linear mixed-effects models. J. Am. Statist.
Assoc. DOI: 10.1080/01621459.2015.1115762.
Zhou, M. (1992). M-estimation in censored linear models. Biometrika 79, 837–841.
Zhu, L. P., Li, L., Li, R. and Zhu, L. X. (2011). Model-free feature screening for
ultrahigh-dimensional data. J. Am. Statist. Assoc. 106, 1464–1475.
Zou, H. (2006). The adaptive lasso and its oracle properties. J. Am. Statist. Assoc. 101,
1418–1429.
Zou, H. (2008). A note on path-based variable selection in the penalized proportional
hazards model. Biometrika 95, 241–247.

18

SUPPLEMENTARY MATERIALS
The supplementary materials include the proof of Theorema 1 and 2.

19

30
25

●

20

●

●
●

15
5

10

MSE

●

●
●
●

●

●

0

●
●

●

FMV

SIS

SIRS

FKS

NIS

QAS1

QAS2

Figure 1: Performance of the proposed model averaging method via different independence
screening methods.

20

15
5

10

MSE

20

25

30

●

●
●

0

●
●

●

●

MCV2

Glasso

●

MAIC

MBIC

MCV1

SCAD

MCP

LASSO

Figure 2: Comparison of the proposed model averaging method with the classical methods
and the penalized methods based on simulated data sets.

21

200

●
●

●
●

●
●

●
●

●
●
●
●

●
●

●
●
●

●
●

●
●
●

●
●
●

●
●
●
●
●
●

●
●
●

●

●

150

●

●
●
●

●
●
●

●
●

●

●
●

●
●
●

100
0

50

ASPE

●

MAIC

MBIC

MCV1

MCV2

Glasso

SCAD

MCP

LASSO

Figure 3: Boxplots for the prediction performance of the proposed model averaging method,
the classical methods and the penalized methods based on the MCL data.

22

