Variable screening with multiple studies
Tianzhou Ma1 , Zhao Ren2 and George C. Tseng1
1

Department of Biostatistics, University of Pittsburgh
2
Department of Statistics, University of Pittsburgh

arXiv:1710.03892v1 [stat.ME] 11 Oct 2017

Abstract
Advancement in technology has generated abundant high-dimensional data that allows integration of multiple relevant studies. Due to their huge computational advantage, variable screening methods based on marginal correlation have become promising
alternatives to the popular regularization methods for variable selection. However, all
these screening methods are limited to single study so far. In this paper, we consider
a general framework for variable screening with multiple related studies, and further
propose a novel two-step screening procedure using a self-normalized estimator for highdimensional regression analysis in this framework. Compared to the one-step procedure
and rank-based sure independence screening (SIS) procedure, our procedure greatly reduces false negative errors while keeping a low false positive rate. Theoretically, we
show that our procedure possesses the sure screening property with weaker assumptions on signal strengths and allows the number of features to grow at an exponential
rate of the sample size. In addition, we relax the commonly used normality assumption
and allow sub-Gaussian distributions. Simulations and a real transcriptomic application illustrate the advantage of our method as compared to the rank-based SIS method.
Key words and phrases: Multiple studies, Partial faithfulness, Self-normalized estimator, Sure screening property, Variable selection

1

Introduction

In many areas of scientific disciplines nowadays such as omics studies (including genomics,
transcriptomics, etc.), biomedical imaging and signal processing, high dimensional data with
much greater number of features than the sample size (i.e. p >> n) have become rule rather
than exception. For example, biologists may be interested in predicting certain clinical
outcome (e.g. survival) using the gene expression data where we have far more genes than
the number of samples. With the advancement of technologies and affordable prices in recent
biomedical research, more and more experiments have been performed on a related hypothesis
or to explore the same scientific question. Since the data from one study often have small
sample size with limited statistical power, effective information integration of multiple studies
can improve statistical power, estimation accuracy and reproducibility. Direct merging of the
data (a.k.a. “mega-analysis”) is usually less favored due to the inherent discrepancy among
the studies (Tseng et al., 2012). New statistical methodologies and theories are required to
solve issues in high-dimensional problem when integrating multiple related studies.
Various regularization methods have been developed in the past two decades and frequently used for feature selection in high-dimensional regression problems. Popular methods
include, but are not limited to, Lasso (Tibshirani, 1996), SCAD (Fan and Li, 2001), elastic
1

net (Zou and Hastie, 2005) and adaptive Lasso (Zou, 2006). When group structure exists
among the variables (for example, a set of gene features belonging to a pre-specified pathway), group version of regularization methods can be applied (Yuan and Lin, 2006; Meier
et al., 2008; Nardi et al., 2008). One can refer to Fan and Lv (2010) and Huang et al. (2012)
for a detailed overview of variable selection and group selection in high-dimensional models.
When the number of features grows significantly larger than the sample size, most regularization methods perform poorly due to the simultaneous challenges of computation expediency,
statistical accuracy and algorithmic stability (Fan et al., 2009). Variable screening methods
become a natural way to consider by first reducing to a lower or moderate dimensional problem and then performing variable regularization. Fan and Lv (2008) first proposed a sure
independent screening (SIS) method to select features based on their marginal correlations
with the response in the context of linear regression models and showed such fast selection
procedure enjoyed a “sure screening property”. Since the development of SIS, many screening methods have been proposed for generalized linear models (Fan et al., 2009, 2010; Chang
et al., 2013), nonparametric additive models or semiparametric models (Fan et al., 2011;
Chang et al., 2016), quantile linear regression (Ma et al., 2017), Gaussian graphical models
(Luo et al., 2014; Liang et al., 2015) or exploit more robust measures for sure screening (Zhu
et al., 2011; Li et al., 2012, 2017). However, all these screening methods are limited to single
study so far.
In this paper, we first propose a general framework for simultaneous variable screening
with multiple related studies. Compared to single study scenario, inclusion of multiple studies
gives us more evidence to reduce dimension and thus increases the accuracy and efficiency
of removing unimportant features during screening. To our knowledge, our paper is the first
to utilize multiple studies to help variable screening in high-dimensional linear regression
model. Such a framework provides a novel perspective to the screening problem and opens
a door to the development of methods using multiple studies to perform screening under
different types of models or with different marginal utilities. In this framework, it is natural
to apply a selected screening procedure to each individual study, respectively. However,
important features with weak signals in some studies may be incorrectly screened out if only
such a one-step screening is performed. To avoid such false negative errors and fully take
advantage of multiple studies, we further propose a two-step screening procedure, where one
additional step of combining studies with potential zero correlation is added to the one-step
procedure for a second check. This procedure has the potential to save those important
features with weak signals in individual studies but strong aggregate effect across studies
during the screening stage. Compared to the naive multiple study extension of SIS method,
our procedure greatly reduces the false negative errors while keeping a low false positive
rate. These merits are confirmed by our theoretical analysis. Specifically, we show that
our procedure possesses the sure screening property with weaker assumptions on signals and
allows the number of features to grow at an exponential rate of the sample size. Furthermore,
we only require the data to have sub-Gaussian distribution via using novel self-normalized
statistics. Thus our procedure can be applied to more general distribution family other than
Gaussian distribution, which is considered in Fan and Lv (2008) and Bühlmann et al. (2010)
for a related screening procedure under single study scenarios. After screening, we further
apply two general and applicable variable selection algorithms: the multiple study extension
of PC-simple algorithm proposed by Bühlmann et al. (2010) as well as a two-stage feature
selection method to choose the final model in a lower dimension.
2

The rest of the paper is organized as follows. In Section 2, we present a framework for
variable screening with multiple related studies as well as notations. Then we propose our
two-step screening procedure in Section 3. Section 4 provides the theoretical properties of our
procedure, and demonstrates the benefits of multiple related studies as well as the advantages
of our procedure. General algorithms for variable selection that can follow from our screening
procedure are discussed in Section 5. Section 6 and 7 include the simulation studies and
a real data application on three breast cancer transcriptomic studies, which illustrate the
advantage of our method in reducing false negative errors and retaining important features
as compared to the rank-based SIS method. We conclude and discuss possible extensions of
our procedure in Section 8. Section 9 provides technical proofs to the major theorems.

2

Model and Notation

Suppose we have data from K related studies, each has n observations. Consider a random
design linear model in each study k ∈ [K] ([K] = 1, . . . , K):
Y

(k)

=

p
X

(k)

(k)

βj Xj + (k) ,

(2.1)

j=1
(k)

(k)

(k)

where each Y (k) ∈ R, each X (k) = (X1 , . . . , Xp )T ∈ Rp with E(X (k) ) = µX and
(k)
cov(X (k) ) = ΣX , each (k) ∈ R with E((k) ) = 0 and var((k) ) = σ 2 such that (k) is
(k)
(k)
(k)
(k)
uncorrelated with X1 , . . . , Xp , and β (k) = (β1 , . . . , βp )T ∈ Rp . We assume implicitly
(k)
with E(Y (k)2 ) < ∞ and E{(Xj )2 } < ∞ for j ∈ [p] ([p] = 1, . . . , p).
When p is very large, we usually assume that only a small set of covariates are true
predictors that contribute to the response. In other words, we assume most of βj =
(1)
(K)
(βj , . . . , βj )T , where j ∈ [p], are equal to a zero vector. In addition, in this paper, we
(k)
assume βj ’s are either zero or non-zero in all K studies. This framework is partially motivated by a high-dimensional linear random effect model considered in literature (e.g.,Jiang
T
, 0T )T , where β(1) is the vector of the
et al. (2016)). More specifically, we can have β = (β(1)
first s0 non-zero components of β (1 ≤ s0 ≤ p). Consider a random effect model where only
(k)
the true predictors of each study are treated as the random effect, that is, β (k) = (β(1) , 0)T
(k)

and β(1) is distributed as N (β(1) , τ 2 Is0 ), where τ 2 is independent of  and X. Consequently,
(k)

βj ’s are either zero or non-zero in all K studies with probability one. Such assumption
fits the reality well, for example, in a typical GWAS study, a very small pool of SNPs are
reported to be associated with a complex trait or disease among millions (Jiang et al., 2016).
With n i.i.d. observations from model (2.1), our purpose is to identify the non-zero β(1) ,
thus we define the following index sets for active and inactive predictors:
(k)

6= 0 for all k};

(k)

= 0 for all k},

A = {j ∈ [p]; βj 6= 0} = {j ∈ [p]; βj

AC = {j ∈ [p]; βj = 0} = {j ∈ [p]; βj

(2.2)

where A is our target. Clearly, under our setting, A and AC are complementary to each
other so that the identification of AC is equivalent to the identification of A. Let |A| = s0 ,
where | · | denotes the cardinality.
3

3

Screening procedure with multiple studies

3.1

Sure independence screening

For a single study (K = 1), Fan and Lv (2008) first proposed the variable screening method
called sure independence screening (SIS) which ranked the importance of variables according
to their marginal correlation with the response and showed its great power in preliminary
screening and dimension reduction for high-dimensional regression problems. Bühlmann
et al. (2010) later introduced the partial faithfulness condition that a zero partial correlation
for some separating set S implied a zero regression coefficient and showed that it held almost
surely for joint normal distribution. In the extreme case when S = ∅, it is equivalent to the
SIS method.
The purpose of sure screening is to identify a set of moderate size d (with d << p) that
will still contain the true set A. Equivalently, we can try to identify AC or subsets of AC
which contain unimportant features that need to be screened out. There are two potential
errors that may occur in any sure screening methods (Fan and Lv, 2010):
1. False Negative (FN): Important predictors that are marginally uncorrelated but
jointly correlated with the response fail to be selected.
2. False Positive (FP): Unimportant predictors that are highly correlated with the
important predictors can have higher priority to be selected than other relatively weaker
important predictors.
The current framework for variable screening with multiple studies is able to relieve us
from the FP errors significantly. Indeed, we have multiple studies in our model setting thus
we have more evidence to exclude noises and reduce FP errors than single study. In addition,
sure screening is used to reduce dimension at a first stage, so we can always include a second
stage variable selection methods such as Lasso or Dantzig selection to further refine the set
and reduce FP errors.
The FN errors occur when signals are falsely excluded after screening. Suppose ρj is
the marginal correlation of the jth feature with the response, with which we try to find the
set {j : ρj = 0} to screen out. Under the assumption of partial faithfulness (for explicit
definition, see Section 4.3), these variables have zero coefficients for sure so the FN errors
are guaranteed to be excluded. However, this might not be true for the empirical version
of marginal correlation. For a single study (K = 1), to rule out the FN errors in empirical
case, it is well-known that the signal-to-noise ratio has to be large (at least of an order of
(log p/n)1/2 after Bonferroni adjustment). In the current setting with multiple studies, the
requirement on strong signals remains the same if we naively perform one-step screening in
each individual study. As we will see next, we propose a novel two-step screening procedure
which allows weak signals in individual studies as long as the aggregate effect is strong
enough. Therefore our procedure is able to reduces FN errors in the framework with multiple
studies.
Before closing this section, it is worthwhile to mention that, to perform a screening test,
one usually applies Fisher’s z-transformation on the sample correlation (Bühlmann et al.,
2010). However, this will require the bivariate normality assumption. Alternatively, in this
paper, we propose to use the self-normalized estimator of correlation that works generally well
even for non-Gaussian data (Shao, 1999). Similar ideas have been applied in the estimation
of large covariance matrix (Cai and Liu, 2016).
4

3.2

Two-step screening procedure with multiple studies
(k)

In the presence of multiple studies, we have more evidence to reduce dimension and ρj = 0
for any k will imply a zero coefficient for that feature. On one hand, it is possible for features
(k)
with zero βj to have multiple non-zero ρj ’s. On the other hand, a non-zero βj will have
(k)
non-zero ρj ’s in all studies. Thus, we aim to identify the following two complementary sets
while performing screening with multiple studies:
(k)

A[0] = {j ∈ [p];

min |ρj | = 0},

A[1] = {j ∈ [p];

min |ρj | =
6 0}.

k

(k)

(3.1)

k

We know for sure that A[0] ⊆ AC and A ⊆ A[1] with the partial faithfulness assumption.
For j ∈ A[0] , the chance of detecting a zero marginal correlation in at least one study has
been greatly increased with increasing K, thus unimportant features will more likely be
screened out as compared to single study scenario.
(k)
One way to estimate A[1] is to test H0 : ρj = 0 of each k for each feature j. When any
of the K tests is not rejected for a feature, we will exclude this feature from Â[1] (we call it
the “One-Step Sure Independence Screening” procedure, or “OneStep-SIS” for short). This
can be viewed as an extension of the screening test to multiple study scenario. However, in
(k)
reality, it is possible for important features to have weak signals thus small |ρj |’s in at least
one study. These features might be incorrectly classified into Â[0] since weak signals can be
indistinguishable from null signals in individual testing. It will lead to the serious problem
of false exclusion of important features (FN) from the final set during screening.
This can be significantly improved by adding a second step to combine those studies
(k)
with potential zero correlation (i.e., fail to reject the null H0 : ρj = 0) identified in the
first step and perform another aggregate test. For the features with weak signals in multiple
studies, as long as their aggregate test statistics is large enough, they will be retained. Such
procedure will be more conservative in screening features as to the first step alone, but will
guarantee to reduce false negative errors.
(k)
(k)
For simplicity, we assume n i.i.d. observations (Xi , Yi ), i ∈ [n], are obtained from
all K studies. It is straightforward to extend the current procedure and analysis to the scenarios with different sample sizes across multiple studies, and thus omitted. Our proposed
“Two-Step Aggregation Sure Independence Screening” procedure (“TSA-SIS” for short) is
formally described below:
Step 1. Screening in each study
In the first step, we perform screening test in each study k ∈ [K] and obtain the estimate
of study set with potential zero correlations ˆlj for each j ∈ [p] as:
√ (k)
nσ̂
(k)
(k)
−1
ˆlj = {k; |T̂ | ≤ Φ (1 − α1 /2)} and T̂ = q j ,
(3.2)
j
j
(k)
θ̂j
P
P
(k)
(k)
(k)
(k)
(k)
(k)
where σ̂j = n1 ni=1 (Xij −X̄j )(Yi −Ȳ (k) ) is the sample covariance and θ̂j = n1 ni=1 [(Xij −
(k)
(k)
(k)
(k)
(k)
X̄j )(Yi − Ȳ (k) ) − σ̂j ]2 . T̂j is the self-normalized estimator of covariance between Xj
5

and Y (k) . Φ is the CDF of standard normal distribution and α1 the pre-specified significance
level.
(k)
In each study, we test if |T̂j | > Φ−1 (1 − α1 /2), if not, we will include study k into
ˆlj . This step does not screen out any variables, but instead separates potential zero and
non-zero study-specific correlations for preparation of the next step. Define the cardinality
of ˆlj as κ̂j = |ˆlj |. If κ̂j = 0 (i.e., no potential zero correlation), we will for sure retain feature
j and not consider it in step 2; Otherwise, we move on to step 2.
(k)

Remark 1. By the scaling property of T̂j , it is sufficient to impose assumptions on the stan(k)

(k)

(k)
dardized variables: W (k) = Y√ −E(Y(k) ) , Zj =
var(Y

)

(k)

(k)

Xj −E(Xj )
q
.
(k)
var(Xj )

(k)

Thus T̂j

(k)

can also be treated
(k)

as the self-normalized estimator of correlation. We thus can define θj = var(Zj W (k) ) and
(k)
(k)
(k)
σj = cov(Zj , W (k) ) = ρj .
Remark 2. In our analysis, the index set in (3.2) is shown to coincide with lj (j ∈ A[0] ) and
lj (j ∈ A[1] ) which will be introduced in more details in Section 4.
Step 2. Aggregate screening
In the second step, we wish to test whether the aggregate effect of potential zero correlations in ˆlj identified in step 1 is strong enough to be retained. Define the statistics
P (k) 2
L̂j =
(T̂j ) and this statistics will approximately follow a χ2κ̂j distribution with degree
k∈l̂j

of freedom κ̂j under null. Thus we can estimate Â[0] by:
Â[0] = {j ∈ [p]; L̂j ≤ ϕ−1
κ̂j (1 − α2 ) and κ̂j 6= 0},

(3.3)

or equivalently estimate Â[1] by:
Â[1] = {j ∈ [p]; L̂j > ϕ−1
κ̂j (1 − α2 ) or κ̂j = 0},

(3.4)

where ϕκ̂j is the CDF of chi-square distribution with degree of freedom equal to κ̂j and α2
the pre-specified significance level.
(k)
The second step takes the sum of squares of T̂j from studies with potential zero corP (k) 2
relation as the test statistics. For each feature j, we test if
(T̂j ) > ϕ−1
κ̂j (1 − α2 ). If
k∈l̂j

rejected, we conclude that the aggregate effect is strong and the feature needs to be retained,
otherwise, we will screen it out. This step performs a second check in addition to the individual testing in step 1 and potentially saves those important features with weak signals in
individual studies but strong aggregate effect.
In Table 1, we use a toy example to demonstrate our idea and compare the two approaches
(“OneStep-SIS” vs. “TSA-SIS”). In this example, suppose we have five studies (K = 5)
and three features (two signals and one noise). “S1” is a strong signal with β = 0.8 in all
studies, “S2” is a weak signal with β = 0.4 in all studies and “N1” is a noise with β = 0.
In hypothesis testing, both small β and zero β can give small marginal correlation and are
sometimes indistinguishable. Suppose T = 3.09 is used as the threshold (corresponding to
6

Table 1: Toy example to demonstrate the strength of two-step screening procedure.

k=1
k=2
k=3
k=4
k=5
ˆlj
κ̂j
TSA-SIS
L̂j
Â[0]
Â[1]
Â[0]
OneStep-SIS
Â[1]

S1 (signal)
S2 (signal)
N1 (noise)
(1)
(1)
(1)
|T̂1 | = 3.71
|T̂2 | = 3.70
|T̂3 | = 0.42
(2)
(2)
(2)
|T̂1 | = 3.16
|T̂2 | = 2.71
|T̂3 | = 0.54
(3)
(3)
(3)
|T̂1 | = 3.46
|T̂2 | = 2.65
|T̂3 | = 0.56
(4)
(4)
(4)
|T̂1 | = 3.63
|T̂2 | = 2.68
|T̂3 | = 0.12
(5)
(5)
(5)
|T̂1 | = 3.24
|T̂2 | = 1.94
|T̂3 | = 0.69
∅
{2, 3, 4, 5}
{1, 2, 3, 4, 5}
0
4
5
25.31 > ϕ4 (0.95) 1.27 < ϕ5 (0.95)
N
N
Y
Y
Y
N
N
Y
Y
Y
N (FN)
N

α1 = 0.001). For the strong signal “S1”, all studies have large marginal correlations, so both
“OneStep-SIS” and “TSA-SIS” procedures include it correctly. For the weak signal “S2”,
since in many studies it has small correlations, it is incorrectly screened out by “OneStepSIS” procedure (False Negative). However, the “TSA-SIS” procedure saves it in the second
step (with α2 = 0.05). For the noise “N1”, both methods tend to remove it after screening.

4
4.1

Theoretical properties
Assumptions and conditions

We impose the following conditions to establish the model selection consistency of our procedure:
(C1) (Sub-Gaussian Condition) There exist some constants M1 > 0 and η > 0 such that for
all |t| ≤ η, j ∈ [p], k ∈ [K]:
(k)2

E{exp(tZj

)} ≤ M1 ,

E{exp(tW (k)2 )} ≤ M1 .
(k)

In addition, there exist some τ0 > 0 such that min θj ≥ τ0 .
j,k

(C2) The number of studies K = O(pb ) for some constant b ≥ 0. The dimension satisfies:
log3 (p) = o(n) and κj log2 p = o(n), where κj is defined next.
(k)

(k)

(C3) For j ∈ A[0] , lj (j ∈ A[0] ) = {k; ρj = 0} and κj = |lj |. If k ∈
/ lj , then |ρj | ≥
q
q
(k)
log p
1.01θj , where C3 = 3(L + 1 + b).
C3
n
7

q
q
(k)
(k)
(C4) For j ∈ A[1] , lj (j ∈ A[1] ) = {k; |ρj | < C1 logn p 0.99θj } and κj = |lj |, where
q
q
(k)
(k)
C1 = L + 1 + b. If k ∈
/ lj , then |ρj | ≥ C3 logn p 1.01θj . In addition, we require
√
P (k) 2 C2 (log2 p+ κj log p)
|ρj | ≥
, where C2 is some large positive constant.
n
k∈lj

(k)

The first condition (C1) assumes that each standardized variable Zj or W (k) , j ∈ [p],
k ∈ [K], marginally follow a sub-Gaussian distribution in each study. This condition relaxes
the normality assumption in (Fan and Lv, 2008; Bühlmann et al., 2010). The second part
of (C1) assumes there always exist some positive τ0 not greater than the minimum variance
(k)
(k)
of Zj W (k) . In particular, if (Xj , Y (k) ) jointly follows a multivariate normal distribution,
(k)
(k)2
then θj = 1 + ρj ≥ 1, so we can always pick τ0 = 1.
The second condition (C2) allows the dimension p to grow at an exponential rate of
sample size n, which is a fairly standard assumption in high-dimensional analysis. Many
sure screening methods like “SIS”, “DC-SIS” and “TPC” have used this assumption (Fan
and Lv, 2008; Li et al., 2012, 2017). Though the PC-simple algorithm (Bühlmann et al.,
2010) assumes a polynomial growth of pn as a function of n, we notice that it can be readily
relaxed to an exponential of n level. Further, we require the product κj log2 p to be small,
which is used to control the errors in the second step of our screening procedure. It is always
true if K log2 p = o(n).
Conditions (C3) assumes a lower bound on non-zero correlation (i.e. k ∈
/ lj ) for features
(k)
[0]
from A . In other words, if the marginal correlation |ρj | is not zero, then it must have a
large enough marginal correlation to be detected. While this has been a key assumption for
a single study in many sure screening methods (Fan and Lv, 2008; Bühlmann et al., 2010;
Li et al., 2012, 2017), we only impose this assumption for j ∈ A[0] rather than all j ∈ [p].
This condition is used to control for type II error in step 1 for features from A[0] .
Condition (C4) gives assumptions on features from A[1] . We assume the correlations to
be small for those k ∈ lj and large for those k ∈
/ lj so that studies with strong or weak
signals can be well separated in the first step. This helps control the type II error in step
1 for features from A[1] . For those studies in lj , we further require their sum of squares of
correlations to be greater than a threshold, so that type II error can be controlled in step 2.
This condition is different from other methods with single study scenario, where they usually
assume a lower bound on each marginal correlation for features from A[1] just like (C3). We
relax this condition and only put restriction on their L2 norm. This allows features from A[1]
to have weak signals in each study but combined strong signal. To appreciate this relaxation,
we compare the minimal requirements with and without step 2. For each j ∈ A[1] , in order to
(k)
detect this feature, we need |ρj | ≥ C(log p/n)1/2 with some large constant C for all k ∈ lj ,
P (k) 2
and thus at least
|ρj | ≥ C 2 κj log p/n. In comparison, the assumption in (C4) is much
k∈lj

weaker in reasonable settings κj >> log p.

4.2

Consistency of the two-step screening procedure

We state the first theorem involving the consistency of screening in our step 1:

8

Theorem 1. Consider a sequence of linear models as in (2.1) which satisfy assumptions and
conditions (C1)-(C4), define the event A = {ˆlj = lj for √
all j ∈ [p]}, there exists a sequence
α1 = α1 (n, p) → 0 as (n, p) → ∞ where α1 = 2{1 − Φ(γ log p)} with γ = 2(L + 1 + b) such
that:
P (A) = 1 − O(p−L ) → 1 as (n, p) → ∞.
(4.1)
The proof of Theorem 1 can be found in Section 9. This theorem states that the screening
in our first step correctly identifies the set lj for features in both A[0] and A[1] (in which strong
and weak signals are well separated) and the chance of incorrect assignment is low. Given
the results in Theorem 1, we can now show the main theorem for the consistency of the
two-step screening procedure:
Theorem 2. Consider a sequence of linear models as in (2.1) which satisfy assumptions
and conditions (C1)-(C4), we know there exists a sequence
α1 = α1 (n, p) → 0 and α2 =
√
α2 (n, p) → 0 as (n, p) → ∞ where α1 = 2{1 −
pΦ(γ log p)} with γ = 2(L + 1 + b) and
α2 = 1 − ϕκj (γκj ) with γκj = κj + C4 (log2 p + κj log p) and some constant C4 > 0 such
that:
P {Â[1] (α1 , α2 ) = A[1] } = 1 − O(p−L ) → 1 as (n, p) → ∞.
(4.2)
The proof of Theorem 2 can be found in Section 9. The result shows that the two-step
screening procedure enjoys the model selection consistency and identifies the model specified
in (3.1) with high
√ probability. The choice of significance level that yields consistency is
α1 = 2{1 − Φ(γ log p)} and α2 = 1 − ϕκj (γκj ) .

4.3

Partial faithfulness and Sure screening property

Bühlmann et al. (2010) first came up with the partial faithfulness assumption which theoretically justified the use of marginal correlation or partial correlation in screening as follows:
ρj|S = 0 for some S ⊆ {j}C implies βj = 0,

(4.3)

where S is the set of variables conditioned on. For independence screening, S = ∅.
Under the two conditions: the positive definiteness of ΣX and non-zero regression coefficients being realization from some common absolutely continuous distribution, they showed
that partial faithfulness held almost surely (Theorem 1 in Bühlmann et al. (2010)). Since
the random effect model described in Section 2 also satisfies the two conditions, the partial
faithfulness holds almost surely in each study.
Thus, we can readily extend their Theorem 1 to a scenario with multiple studies:
Corollary 1. Consider a sequence of linear models as in (2.1) satisfying the partial faithfulness condition in each study and true active and inactive set defined in (2.2), then the
following holds for every j ∈ [p]:
(k)

ρj|S = 0 for some k for some S ⊆ {j}C implies βj = 0.

9

(4.4)

(k)

The proof is straightforward and thus omitted: if ρj|S = 0 for some study k, then with
(k)

partial faithfulness, we will have βj = 0 for that particular k. Since we only consider
(k)
features with zero or non-zero βj ’s in all studies in (2.2), we will have βj = 0. In the case
(k)
of independence screening (i.e. S = ∅), ρj = 0 for some k will imply a zero βj .
With the model selection consistency in Theorem 2 and the extended partial faithfulness
condition in Corollary 1, the sure screening property of our two-step screening procedure
immediately follows:
Corollary 2. Consider a sequence of linear models as in (2.1) which satisfy assumptions
and conditions (C1)-(C4) as well as the extended partial faithfulness condition in Corollary
1, there exists a sequence
√ α1 = α1 (n, p) → 0 and α2 = α2 (n, p) → 0 as (n, p) → ∞
where α1 = 2{1 − Φ(γ log p)} with γ = 2(L + 1 + b) and α2 = 1 − ϕκj (γκj ) with γκj =
p
κj + C4 (log2 p + κj log p) such that:
P {A ⊆ Â[1] (α1 , α2 )} = 1 − O(p−L ) → 1 as (n, p) → ∞.

(4.5)

The proof of this Corollary simply combines the results of Theorem 2 and the extended
partial faithfulness and is skipped here.

5

Algorithms for variable selection with multiple studies

Usually, performing sure screening once may not remove enough unimportant features. In our
case since there are multiple studies, we expect our two-step screening procedure to remove
many more unimportant features than in single study. If the dimension is still high after
applying our screening procedure, we can readily extend the two-step screening procedure
to an iterative variable selection algorithm by testing the partial correlation with gradually
increasing size of the conditional set S. Since such method is a multiple study extension of
the PC simple algorithm in Bühlmann et al. (2010), we call it “Multi-PC” algorithm (Section
5.1).
On the other hand, if the dimension has already been greatly reduced with the two-step
screening, we can simply add a second stage group-based feature selection techniques to
select the final set of variables (Section 5.2).

5.1

Multi-PC algorithm

We start from S = ∅, i.e., our two-step screening procedure and build a first set of candidate
active variables:
Â[1,1] = Â[1] = {j ∈ [p]; L̂j > ϕ−1
κ̂j (1 − α2 ) or κ̂j = 0}.

(5.1)

We call this set stage1 active set, where the first index in [, ] corresponds to the stage
of our algorithm and the second index corresponds to whether the set is for active variables
([, 1]) or inactive variables ([, 0]). If the dimensionality has already been decreased by a large

10

amount, we can directly apply group-based feature selection methods such as group lasso to
the remaining variables (to be introduced in Section 5.2).
However, if the dimension is still very high, we can further reduce dimension by increasing
the size of S and considering partial correlations given variables in Â[1,1] . We follow the
similar two-step procedure but now using partial correlation of order one instead of marginal
correlation and yield a smaller stage2 active set:
Â[2,1] = {j ∈ Â[1,1] ; L̂j|q > ϕκ̂−1j|q (1 − α2 ) or κ̂j|q = 0, for all q ∈ Â[1,1] \{j}},

(5.2)

where each self-normalized estimator of partial correlation can be computed by taking the
residuals from regressing over the variables in the conditional set.
We can continue screening high-order partial correlations, resulting in a nested sequence
of m active sets:
Â[m,1] ⊆ . . . ⊆ Â[2,1] ⊆ Â[1,1] .
(5.3)
Note that the active and inactive sets at each stage are non-overlapping and the union
of active and inactive sets at a stage m will be the active set in a previous stage m − 1, i.e.,
Â[m,1] ∪ Â[m,0] = Â[m−1,1] . This is very similar to the original PC-simple algorithm, but now
at each order-level, we perform the two-step procedure. The algorithm can stop at any stage
m when the dimension of Â[m,1] already drops to low to moderate level and other common
group-based feature selection techniques can be used to select the final set. Alternatively,
we can continue the algorithm until the candidate active set does not change anymore. The
algorithm can be summarized as follows:
Algorithm 1. Multi-PC algorithm for variable selection.
Step 1. Set m = 1, perform the two-step screening procedure to construct stage1 active set:
Â[1,1] = {j ∈ [p]; L̂j > ϕ−1
κ̂j (1 − α2 ) or κ̂j = 0}.
Step 2. Set m = m + 1. Construct the stagem active set:
Â[m,1] = {j ∈ Â[m−1,1] ; L̂j|S > ϕ−1
κ̂j|S (1 − α2 ) or κ̂j|S = 0,
for all S ⊆ Â[m−1,1] \{j} with |S| = m − 1}.
Step 3. Repeat Step 2 until m = m̂reach , where m̂reach = min{m : |Â[m,1] | ≤ m}.

11

5.2

Two-stage feature selection

As an alternative to “Multi-PC” algorithm for variable selection, we also introduce here a
two-stage feature selection algorithm by combining our two-step screening procedure and
other regular feature selection methods together. In single study, for example, Fan & Lv
(2008) performed sure independence screening in the first stage followed by model selection
techniques including Adaptive Lasso, Dantzig Selector and SCAD, etc., and named those
procedures as “SIS-AdaLasso”,“SIS-DS”, “SIS-SCAD” , accordingly.
In our case, since the feature selection is group-based, we adopt a model selection technique using group Lasso penalty in the second stage:
min
β

K
X

(k)

(k)

||y (k) − XÂ[1] βÂ[1] ||22 + λ

k=1

X

||βj ||2

,

(5.4)

j∈Â[1]

where Â[1] is the active set identified from our two-step screening procedure and the tuning
parameter λ can be chosen by cross-validation or BIC in practice just like for a regular group
Lasso problem. We call such two-stage feature selection algorithm as “TSA-SIS-groupLasso”.
In addition, at any stages of the “Multi-PC” algorithm when the dimension has already
been dropped to a moderate level, the group Lasso-based feature selection techniques can
always take over to select the final set of variables.

6

Numerical evidence

In this section, we demonstrate the advantage of TSA-SIS procedure in comparing to the
multiple study extension of SIS (named “Min-SIS”), which ranks the features by the minimum absolute correlation among all studies. We simulated data according to the linear
(k)
model in (2.1) including p covariates with zero mean and covariance matrix Σi,j = r|i−j|
(k)
(k)
where Σi,j denotes the (i, j)th entry of ΣX .
In the first part of simulation, we fixed the sample size n = 100, p = 1000, the number of
studies K = 5 and performed B = 1000 replications in each setting. We assumed that the
true active set consisted of only ten variables and all the other variables had zero coefficients
(i.e., s0 = 10). The indices of non-zero coefficients were evenly spaced between 1 and p.
The variance of the random error term in linear model was fixed to be 0.52 . We randomly
drew r from {0, 0.2, 0.4, 0.6} and allowed different r’s in different studies. We considered the
following four settings:
1. Homogeneous weak signals across all studies: nonzero βj generated from Unif(0.1, 0.3)
(1)
(2)
(K)
and βj = βj = . . . = βj = βj .
2. Homogeneous strong signals across all studies: nonzero βj generated from Unif(0.7, 1)
(1)
(2)
(K)
and βj = βj = . . . = βj = βj .
3. Heterogeneous weak signals across all studies: nonzero βj generated from Unif(0.1, 0.3)
(k)
and βj ∼ N (βj , 0.52 ).

12

Table 2: Sensitivity analysis on the choice of α1 and α2 in simulation (Sensitivity/Specificity)
Sensitivity/Specificity
α1 =0.01
0.001
0.0001

α2 = 0.15
0.05
0.01
0.001
0.793/0.901 0.525/0.984 0.210/0.999 0.142/1.000
0.947/0.826 0.864/0.943 0.691/0.990 0.373/0.999
0.966/0.816 0.922/0.932 0.840/0.985 0.681/0.998

Note: All value are based on average results from B = 1000 replications.

4. Heterogeneous strong signals across all studies: nonzero βj generated from Unif(0.7, 1)
(k)
and βj ∼ N (βj , 0.52 ).
We evaluated the performance of Min-SIS using receiver operating characteristic (ROC)
curves which measured the accuracy of variable selection independently from the issue of
choosing good tuning parameters (for Min-SIS, the tuning parameter is the top number of
features d). The OneStep-SIS procedure we mentioned above was actually one special case
of the Min-SIS procedure (by thresholding at α1 ). In presenting our TSA-SIS procedure,
we fixed α1 = 0.0001 and α2 = 0.05 so the result was just one point on the sensitivity vs.
1-specificity plot. We also performed some sensitivity analysis on the two cutoffs based on
the first simulation (see Table 2) and found the two values to be optimal since they had
both high sensitivity and high specificity. Thus we suggested fixing these two values in all
the simulations.
Figure 1 showed the results of simulation 1-4. When the signals were homogeneously weak
in all studies as in (1), TSA-SIS clearly outperformed the Min-SIS procedure (above its ROC
curve). It reached about 90% sensitivity with controlled false positive errors (specificity
∼ 95%). In order to reduce false negatives, Min-SIS had to sacrifice the specificity and
increased the false positives, which in the end lost the benefits of performing screening
(i.e. end up keeping too many features). When the signals became strong as in (2), both
procedures performed equally well. This fit our motivation and theory and showed the
strength of our two-step procedure in saving weak signals without much increase in false
positive rates. When the signals became heterogeneous as in (3) and (4), both procedures
performed worse than before. But the Min-SIS procedure never outperformed the TSASIS procedure since it only examined the minimum correlation among all studies while the
two-step procedure additionally considered the aggregate statistics.

7

Real data application

We next demonstrated our method in three microarray datasets of triple-negative breast
cancer (TNBC, sometimes a.k.a. basal-like), an aggressive subtype of breast cancer usually
with poor prognosis. Previous studies have shown that the tumor suppressor protein “p53”
played an important role in breast cancer prognosis and its expression was associated with
both disease-free survival and overall survival in TNBC (Yadav et al., 2015). Our purpose
was to identify the genes most relevant and predictive to the response - the expression level
of TP53 gene, which encodes p53 protein. The three datasets are publicly available on
authors’ website or at GEO repository including METABRIC (a large cohort consisting
of roughly 2000 primary breast tumours), GSE25066 and GSE76250 (Curtis et al., 2012;
13

Figure 1: Simulation results 1-4: the ROC curve is for Min-SIS, the black point is for our
TSA-SIS using α1 = 0.0001 and α2 = 0.05.

1.0

(2) Homogeneous strong

1.0

(1) Homogeneous weak
●

0.8
0.6

Sensitivity

0.0

0.2

0.4

0.6
0.4
0.0

0.2

Sensitivity

0.8

●

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

(3) Heterogeneous weak

(4) Heterogeneous strong

1.0

1.0

1−Specificity

1.0

1−Specificity

0.8
0.6

Sensitivity

0.0

0.2

0.4

0.6
0.4
0.0

0.2

Sensitivity

0.8

●
●

0.0

0.2

0.4

0.6

0.8

1.0

0.0

1−Specificity

0.2

0.4

0.6

0.8

1.0

1−Specificity

Itoh et al., 2014; Liu et al., 2016). We subset the data to focus on the TNBC cases only
and ended up with 275, 178 and 165 TNBC samples in each dataset, respectively. After
routine preprocessing and filtering by including genes sufficiently expressed and with enough
variation, a total of 3377 genes remained in common for the analysis.
We applied our Multi-PC algorithm and compared to the OneStep-SIS procedure as
well as the Min-SIS method by using d = n/ log(n) = 49 (as suggested by their paper).
We used α1 = 0.0001 and α2 = 0.05 (as determined by sensitivity analysis in simulation)
and the “Multi-PC” algorithm only ran up to the first order (i.e. m = 2) and stopped
with six features. This again showed the power of screening with multiple studies. After
feature selection, we fit the linear model in each study to obtain the coefficient estimates and
adjusted R2 . Table 3 showed the coefficient estimates and standard errors of the final set of
six genes selected by our procedure. We added three columns to indicate whether they were
also retained by the Min-SIS (and their relative rank) or OneStep-SIS procedures. As we
can see from the table, all the six genes selected by our procedure were missed by the other
14

Table 3: The six genes selected by our TSA-SIS procedure.
Gene

METABRIC
Est (SE)
Intercept
7.600 (1.502)
EXOC1
0.251 (0.081)∗∗
ITGB1BP1 -0.134 (0.045)∗∗
RBM23
0.168 (0.078)∗
SETD3
-0.166 (0.081)∗
SQSTM1 -0.114 (0.050)∗
TRIOBP -0.126 (0.062)∗
Adjusted-R2
0.151

GSE25066
Est (SE)
0.213 (0.553)
0.278 (0.157).
0.003 (0.111)
0.144 (0.167)
0.366 (0.184)∗
0.029 (0.099)
0.084 (0.118)
0.522

GSE76250 Min-SIS Rank in OneStep-SIS
Est (SE)
d=49 Min-SIS
|S|=25
-1.783 (0.971)
.
0.293 (0.167)
N
164
N
-0.178 (0.194)
N
123
N
0.367 (0.168)∗
N
152
N
-0.080 (0.175)
N
101
N
0.245 (0.183)
N
98
N
0.628 (0.261)∗
N
91
N
0.359

Note: “.” indicates significant level of 0.1, “∗” for level of 0.05, “∗∗” for level of 0.01.

methods. Those genes typically had weak signals in one or more studies thus were very likely
to be incorrectly excluded if only one step screening is performed. Since the METABRIC
study had a larger sample size, all the coefficients appeared to be more significant than the
other two studies.
The gene EXOC1 and p53 are both components of the Ras signaling pathway which
is responsible for cell growth and division and can ultimately lead to cancer (Rajalingam
et al., 2007). RBM23 encodes for an RNA-binding protein implicated in the regulation of
estrogen-mediated transcription and has been found to be associated with p53 indirectly via
a heat shock factor (Asano et al., 2016). ITGB1BP1 encodes for an integrin protein which is
essential for cell adhesion and other downstream signaling pathways that are also modulated
by p53 (Brakebusch et al., 2002).

8

Discussion

In this paper, we proposed a two-step screening procedure for high-dimensional regression
analysis with multiple related studies. In a fairly general framework with weaker assumptions
on the signal strength, we showed that our procedure possessed the sure screening property
for exponentially growing dimensionality without requiring the normality assumption. We
have shown through simulations that our procedure consistently outperformed the rankbased SIS procedure independent of their tuning parameter d. As far as we know, our paper
is the first proposed procedure to perform variable screening in high-dimensional regression
when there are multiple related studies. In addition, we also introduced two applicable
variable selection algorithms following the two-step screening procedure.
Variable selection in regression with multiple studies have been studied in a subfield of
machine learning called multi-task learning (MTL) before and the general procedure is to
apply regularization methods by putting group Lasso penalty, fused Lasso penalty or trace
norm penalty, etc. (Argyriou et al., 2007; Zhou et al., 2012; Ji and Ye, 2009). However, at
ultra-high dimension, such regularization methods usually fail due to challenges in computation expediency, statistical accuracy and algorithmic stability. Instead, sure screening can be

15

used as a fast algorithm for preliminary feature selection, and as long as it exhibits comparable statistical performance both theoretically and empirically, its computational advantages
make it a good choice in application (Genovese et al., 2012). Our method has provided an
alternative to target the high-dimensional multi-task learning problems.
The current two-step screening procedure is based on the linear models but relaxes the
Gaussian assumption to sub-Gaussian distribution. One can apply a modified Fisher’s ztransformation estimator rather than our self-normalized estimator to readily accommodate
general elliptical distribution families (Li et al., 2017). In biomedical applications, noncontinuous outcomes such as categorical, count or survival outcomes are more commonly
seen. Fan et al. (2010) extended SIS and proposed a more general independent learning
approach for generalized linear models by ranking the maximum marginal likelihood estimates. Fan et al. (2011) further extended the correlation learning to marginal nonparametric
learning for screening in ultra-high dimensional additive models. Other researchers exploited
more robust measure for the correlation screening (Zhu et al., 2011; Li et al., 2012; Balasubramanian et al., 2013). All these measures can be our potential extension by modifying the
marginal utility used in the screening procedure. Besides, the idea of performing screening
with multiple studies is quite general and is applicable to relevant statistical models other
than the regression model, for example, Gaussian graphical model with multiple studies. We
leave these interesting problems in future study.

9

Proofs

We start by introducing three technical lemmas that are essential for the proofs of the main
(k)
results. By the scaling property of T̂j and Remark 1, without loss of generality, we can
(k)
(k)
assume E(Xj ) = E(Y (k) ) = 0 and var(Xj ) = var(Y (k) ) = 1 for all k ∈ [K], j ∈ [p].
(k)
(k)
Therefore in the proof we do not distinguish between σj and ρj . The first lemma is on
(k)
the concentration inequalities of the self-normalized covariance and θ̂j .
Lemma 1. Under the assumptions (C1) and (C2), for any δ ≥ 2 and M > 0, we have:
q
(k)
(k)
σ̂ −σj
(i) P (max | j (k) 1/2
| ≥ δ logn p ) = O((log p)−1/2 p−δ+1+b ),
j,k

(θ̂j )
(k)

(k)

(ii) P (max |θ̂j − θj | ≥ Cθ

q

j,k

log p
)
n

= O(p−M ),

where Cθ is a positive constant depending on M1 , η and M only.
The second and third lemmas, which will be used in the proof of Theorem 2, describe
n
P
r (k)
(k)
(k)
(k)
(k)
√ (k)
√1
[(Xij −X̄j )(Yi −Ȳ (k) )−ρj ]
n
θ̂j
nρ
(k)
(k)
i=1
q
q j
the concentration behaviors of Ĥj :=
=
T̂
−
(k)
j
(k)
(k)
θj

n
P

(k)

and Ȟj

:=

(k) (k)
(k)
√1
(Xij Yi −ρj )
n
i=1
q
(k)
θj

.

16

θj

θj

Lemma 2. There exists some constant c > 0 such that,
P (|

X

(k)2

[Ȟj

− 1]| > t) ≤ 2 exp(−c min[

k∈lj

t2 1/2
, t ]),
κj

where c depends on M1 and η only.
Lemma 3. There exists some constant CH > 0 such that,
s
log2 p
(k)
(k)
P (max |Ȟj − Ĥj | > CH
) = O(p−M ),
j,k
n
s
log3 p
(k)2
(k)2
P (max |Ȟj − Ĥj | > CH
) = O(p−M ),
j,k
n
where CH depends on M1 , η, M and τ0 only.
The proofs of the three lemmas are provided in the Appendix.
Proof of Theorem 1. We first define the following error events:
[0]

I,A
Ej,k

[0]

II,A
Ej,k

[1]

I,A
Ej,k

[1]

II,A
Ej,k

(k)

= {|T̂j | > Φ−1 (1 − α1 /2) and j ∈ A[0] , k ∈ lj },
(k)

/ lj },
= {|T̂j | ≤ Φ−1 (1 − α1 /2) and j ∈ A[0] , k ∈
(k)

= {|T̂j | > Φ−1 (1 − α1 /2) and j ∈ A[1] , k ∈ lj },
(k)

/ lj }.
= {|T̂j | ≤ Φ−1 (1 − α1 /2) and j ∈ A[1] , k ∈

To show Theorem 1 that P (A) = 1 − O(p−L ), it suffices to show that,
[ I,A[0]
II,A[0]
P { (Ej,k
∪ Ej,k
)} = O(p−L ),

(9.1)

j,k

and
P{

[ I,A[1]
II,A[1]
(Ej,k ∪ Ej,k
)} = O(p−L ).

(9.2)

j,k

One√can apply Lemma 1 to bound each component in (9.1) and (9.2) with α1 = 2{1 −
Φ(γ log p)} and γ = 2(L + 1 + b). Specifically, we obtain that,
[ I,A[0]
P ( Ej,k
) = P ( max
j,k

j∈A[0] ,k∈lj

(k)

|T̂j | ≥ γ

p
1
log p) = O( √
p−γ+1+b ) = o(p−L ),
log p

17

(9.3)

(k)

where the second equality is due to Lemma 1 (i) with δ = γ, noting that σj
√ (k) q (k)
(k)
T̂j = nσ̂j / θ̂j . In addition, we have that,
[ I,A[1]
P ( Ej,k
) =P { max

j∈A[1] ,k∈lj

j,k

(k)

|T̂j | ≥ γ
(k)

≤P ( max

|

p

j∈A[1] ,k∈lj

log p}

(k)

σ̂j − ρj

r
| ≥ (γ − C1 )

(k)
(θ̂j )1/2

= 0 and

log p
) + O(p−L )
n

(9.4)

1
p−(γ−C1 )+1+b ) + O(p−L )
log p
=O(p−L ),
=O( √

where the inequality on the second line is due to assumption (C4) on lj for j ∈ A[1] , Lemma
(k)
(k)
(k)
1 (ii) with M = L, and assumption (C1) minj,k θj ≥ τ0 , i.e., θ̂j ≥ θj − Cθ (log p/n)1/2 ≥
(k)
0.99θj . The equality on the third line follows from Lemma 1 (i) where δ = γ −C1 = L+1+b.
In the end, we obtain that,
[ II,A[0]
p
(k)
II,A[1]
P { (Ej,k
∪ Ej,k
)} =P (max |T̂j | < γ log p)
j,k∈l
/j

j,k

(k)

≤P (max |
j,k∈l
/j

(k)

σ̂j − ρj

(k)
(θ̂j )1/2

r
| ≥ (C3 − γ)

log p
) + O(p−L )
n

(9.5)

1
p−(C3 −γ)+1+b ) + O(p−L )
log p
−L
=O(p ),
=O( √

where the inequality on the second line is due to assumptions (C3) and (C4) on lj , Lemma
(k)
(k)
1 (ii) with M = L and assumption (C1) on sub-Gaussian distributions, i.e., θ̂j ≤ θj +
(k)
(k)
Cθ (log p/n)1/2 ≤ 1.01θj . In particular, we have implicitly used the fact that maxj,l θj is
upper bounded by a constant depending on M1 and η only. The equality on the third line
follows from Lemma 1 (i) where δ = C3 − γ = L + 1 + b.
Finally, we complete the proof by combining (9.3)-(9.5) to show (9.1)-(9.2).
Proof of Theorem 2. We first define the following error events:
[0] ,2

EjA

[1] ,2

EjA

= {|L̂j | > ϕ−1 (1 − α2 ) or κ̂j = 0} for j ∈ A[0] ,

= {|L̂j | < ϕ−1 (1 − α2 ) and κ̂j 6= 0} for j ∈ A[1] .

To prove Theorem 2, we only need to show that,
[ A[0] ,2
[ A[1] ,2
Ej
Ej
P(
) = O(p−L ) and P (
) = O(p−L ),
j∈A[0]

with α2,κj := 1 − ϕκj [κj + C4 (log2 p +

j∈A[1]

p
κj log p)] := 1 − ϕκj (γκj ).
18

(9.6)

Recall the event A defined in Theorem 1. Thus we have that,
[
[0]
[1]
P {(∪j∈A[0] EjA ,2 ) (∪j∈A[1] EjA ,2 )}
X (k)2
X
≤P (AC ) + p max P (
T̂j > γκj ) + p max P (
T̂ (k)2 < γκj ).
j∈A[0]

j∈A[1] ,κj 6=0

k∈lj

k∈lj

Therefore, given the results in Theorem 1, it suffices to show,
X
P(
T̂ (k)2 > γκj ) = O(p−L−1 ) for any j ∈ A[0] ,

(9.7)

k∈lj

and
P(

X

T̂ (k)2 < γκj ) = O(p−L−1 ) for any j ∈ A[1] and κj > 0.

(9.8)

k∈lj
[0]

We first prove equation (9.7). Since j ∈ A , we have
P
(k)2
to bound the probability of k∈lj T̂j > γκj below.
P(

X

(k)2

T̂j

> γκj )

(k)2
Ĥj

Cθ
> (1 −
τ0

(k)
Ĥj

=

(k)
T̂j

r

(k)

θ̂j

(k)

θj

. We are ready

k∈lj

≤P (

X
k∈lj

r

log p
)γκj ) + O(p−L−1 )
n
r

s

log3 p
) + O(p−L−1 )
n
j
s
X (k)2
p
κ2j log p
C
θ
2
=P ( (Ȟj − 1) > κj + C4 (log p + κj log p) −
τ0
n
k∈lj
s
s
s
Cθ C4
log5 p
κj log2 p
log3 p
−
(
+
) − κj − κj CH
) + O(p−L−1 )
τ0
n
n
n
X (k)2
p
≤P ( (Ȟj − 1) > C20 (log2 p + κj log p)) + O(p−L−1 )
X (k)2
Cθ
≤P ( (Ȟj − 1) > (1 −
τ0
k∈l

log p
)γκj − κj − κj CH
n

k∈lj

=O(p−L−1 ).
(k)

The inequality on the second line is due to assumption (C1) that min θj
j,k

≥ τ0 > 0 and

Lemma 1 (ii) with M = L + 1. The inequality on the third line follows from Lemma 3 with
M = L + 1. The inequality on the fifth line is by the choice of γκj with a sufficiently large
C4 > 0 and the assumption (C2) that log3 p = o(n) and κj log2 p = o(n). The last equality
follows from Lemma 2.

19

Lastly, we prove (9.8) as follows,
X (k)2
P(
T̂j < γκj )
k∈lj

√ (k)
(k)
X (k)
nρj 2 θj
) (k) < γκj )
=P ( (Ĥj + q
(k)
θ̂j
k∈lj
θj
r
√ (k)
X (k)
nρj 2
Cθ log p
)γκj ) + O(p−L−1 )
≤P ( (Ĥj + q
) ≤ (1 +
(k)
τ
n
0
k∈lj
θj
s
r
X (k)2
X (k)2
log3 p
Cθ log p
≤P ( (Ȟj − 1) ≤ κj CH
− κj + (1 +
)γκj − Cm n
ρj
n
τ
n
0
k∈lj
k∈lj
s
√ (k)
√
(k)
X (k) nρj
log2 p X n|ρj |
q
Ȟj q
+ 2CH
) + O(p−L−1 ).
−2
(k)
(k)
n
k∈lj
k∈lj
θj
θj
(k)

The inequality on the third line is due to assumption (C1) that min θj
j,k

(9.9)

≥ τ0 > 0 and

Lemma 1 (ii) with M = L + 1. The inequality on the fourth line follows from Lemma 3
(k)
(both equations) and min(θj )−1 := Cm > 0, guaranteed by the sub-Gaussian assumption
j,k

in assumption (C1).
We can upper bound the term 2CH
s
2CH

q

log2 p
n

√ (k)
n|ρ |
q j
k∈lj
(k)
θj

P

in (9.9) as follow,

s
√ (k)
√
s X
sX
X
n|ρj |
log p
log2 p n √
(k)2
(k)2
q
ρj = o( n
≤ 2CH
κj
ρj ).
√
(k)
n k∈l
n
τ0
k∈lj
k∈lj
θj
j
2

(9.10)

The first inequality is by the Cauchy-Schwarz inequality and assumption (C1), and the
second equality by the assumption (C2) that κj log2 p = o(n).
√ (k)
P
(k)
(k) nρj
We next upper bound the term −2 k∈lj Ȟj q (k)
with high probability. Note that θj
θj

(k)
θj

(k)

−1
Cm

is bounded below and above, i.e., τ0 ≤
≤
by assumption (C1). In addition, Ȟj has
zero mean and is sub-exponential with bounded constants by assumption (C1). By Bernstein
inequality (Proposition 5.16 in Vershynin (2010)), we have with some constant c0 > 0,
P (|2

X
k∈lj

We pick t = CB

(k)
|Ȟj

√ (k)
n|ρ |
t2
t
q j | > t) ≤ 2 exp(−c0 min[ P (k)2 ],
).
√
(k)
(k)
n
ρj
max n|ρj |
θj
k∈lj
k∈lj

q P
(k)2
n k∈lj ρj log2 p with a large constant CB in the inequality above and

20

apply (9.10) to reduce (9.9) as follows,
X (k)2
P(
T̂j < γκj )
k∈lj

s X
X (k)2
X (k)2
(k)2
≤P ( (Ȟj − 1) ≤ −Cm n
ρj log2 p
ρj + 2CB n
k∈lj

k∈lj

k∈lj

p
+ 2C4 κj log p + 2C4 log2 p) + O(p−L−1 )
q
X (k)2
p
p
2
≤P ( (Ȟj − 1) ≤ −Cm C2 (log p + κj log p) + 2CB C2 log2 p(log2 p + κj log p)
k∈lj

p
+ 2C4 κj log p + 2C4 log2 p) + O(p−L−1 )
X (k)2
p
≤P ( (Ȟj − 1) ≤ −C20 (log2 p + κj log p)) + O(p−L−1 )
k∈lj

=O(p−L−1 ).
The inequality on the first line is obtained by the choice of γκj with the chosen C4 > 0
and the assumption (C2) that κj log2 p = o(n). The inequalities
√ on the second line and third
P
C2 (log2 p+ κj log p)
(k) 2
for a sufficiently large
line are by the assumption (C4) that k∈lj |ρj | ≥
n
C2 > 0. The last equality is by Lemma 2.
This completes the proof of (9.7) and (9.8), which further yields to
[
[0]
[1]
P {(∪j∈A[0] EjA ,2 ) (∪j∈A[1] EjA ,2 )} = O(p−L ),
with the results from Theorem 1. Therefore we complete the proof of Theorem 2.

Appendix
S1. Proof of Lemma 1
Proof. Part (i) immediately follows from Lemma 2 (i) equation (25) in Cai and Liu (2011).
To prove part (ii), we need to bound the three terms on the right side of the following
inequality,
(k)

(k)

(k)

(k)

(k)

(k)

(k)

(k)

max |θ̂j − θj | ≤ max |θ̂j − θ̃j | + max |θ̃j − θ̌j | + max |θ̌j − θj |,
j,k

(k)

where θ̃j
(k)

:=

j,k

1
n

n
P

(k)

(k)

(Xij Yi

j,k

(k)

(k)

− ρ̃j )2 with ρ̃j =

i=1
(k)

1
n

(A1)

j,k

n
P

(k)

(k)

(k)

Xij Yi , and θ̌j

i=1

(k)

:=

1
n

n
P

(k)

(k)

(Xij Yi

−

i=1

ρj )2 . Note that E(θ̌j ) = θj .
By the marginal sub-Gaussian distribution assumption in assumption (C1), we have that
(k) (k)
(k)
(k)
(Xij Yi − ρj )2 has mean θj and finite Orlicz ψ1/2 -norm (see, e.g., Adamczak et al.
(2011)). Thus we can apply equation (3.6) of Adamczak et al. (2011), i.e.,
√ (k)
t2
(k)
P (max n|θ̌j − θj | > t) ≤ 2 exp(−c min[ , t1/2 ]),
j,k
n
21

√
with t = (Cθ /3) n log p for a large enough constant Cθ > 0 depending on M1 , η and M only
to obtain that,
r
log p
(k)
(k)
P (max |θ̌j − θj | > (Cθ /3)
) = O(p−M ).
(A2)
j,k
n
2

We have used the assumption log p = o(n1/3 ) in assumption (C2) to make sure tn ≤ t1/2 .
By applying equation (1) in supplement of Cai and Liu (2011), we obtain that,
r
log p
(k)
(k)
P (max |θ̃j − θ̂j | > (Cθ /3)
) = O(p−M ).
(A3)
j,k
n
In addition, by a similar truncation argument as that in the proof of Lemma 2 in Cai and
Liu (2011) and equation (7) therein, we obtain that by picking a large enough Cθ > 0,
r
log p
(k)
(k)
P (max |θ̃j − θ̌j | > (Cθ /3)
) = O(p−M ).
(A4)
j,k
n
We complete the proof by combining (A1)-(A4) with a union bound argument.

S2. Proof of Lemma 2
(k)

(k)

Proof. It is easy to check that E(Ȟj ) = 0 and var(Ȟj ) = 1. The marginal sub-Gaussian
(k)
distribution assumption in assumption (C1) implies that Ȟj has finite Orlicz ψ1 -norm (i.e.,
(k)
sub-exponential distribution with finite constants). Therefore, (Ȟj )2 −1 is centered random
(k)
variable with finite Orlicz ψ1/2 -norm. Note that Ȟj are independent for k ∈ [K]. The result
follows from equation (3.6) of Adamczak et al. (2011).

S3. Proof of Lemma 3
√
(k) √
√ (k)
nX̄j
nȲ (k)
(k)
(k)
q
Proof. Note that Ȟj −Ĥj =
. By assumption (C1), we have that E( nX̄j ) =
(k)
nθj
√ (k)
√ (k)
√
√ (k)
√
E( nȲ ) = 0, var( nX̄j ) = var( nȲ (k) ) = 1, and both nX̄j and nȲ (k) are subGaussian with bounded constants. Therefore, the first equation follows from Bernstein
inequality (e.g., Definition 5.13 in Vershynin (2010)) applied to centered sub-exponential
√ (k) √
(k)
variable nX̄j · nȲ (k) , noting θj ≥ τ0 by assumption (C1). The second equation follows
from the first one, log3 p = o(n), and a Bernstein inequality (e.g., Corollary 5.17 in Vershynin
(k)
(2010)) applied to the sum of centered sub-exponential variables Ȟj .

References
Adamczak, R., Litvak, A. E., Pajor, A., and Tomczak-Jaegermann, N. (2011). Restricted
isometry property of matrices with independent columns and neighborly polytopes by
random sampling. Constructive Approximation, 34(1):61–88.

22

Argyriou, A., Evgeniou, T., and Pontil, M. (2007). Multi-task feature learning. In Advances
in neural information processing systems, pages 41–48.
Asano, Y., Kawase, T., Okabe, A., Tsutsumi, S., Ichikawa, H., Tatebe, S., Kitabayashi,
I., Tashiro, F., Namiki, H., Kondo, T., et al. (2016). Ier5 generates a novel hypophosphorylated active form of hsf1 and contributes to tumorigenesis. Scientific reports,
6.
Balasubramanian, K., Sriperumbudur, B., and Lebanon, G. (2013). Ultrahigh dimensional
feature screening via rkhs embeddings. In Artificial Intelligence and Statistics, pages 126–
134.
Brakebusch, C., Bouvard, D., Stanchi, F., Sakai, T., and Fassler, R. (2002). Integrins in
invasive growth. The Journal of clinical investigation, 109(8):999.
Bühlmann, P., Kalisch, M., and Maathuis, M. H. (2010). Variable selection in highdimensional linear models: partially faithful distributions and the pc-simple algorithm.
Biometrika, 97(2):261–278.
Cai, T. and Liu, W. (2011). Adaptive thresholding for sparse covariance matrix estimation.
Journal of the American Statistical Association, 106(494):672–684.
Cai, T. T. and Liu, W. (2016). Large-scale multiple testing of correlations. Journal of the
American Statistical Association, 111(513):229–240.
Chang, J., Tang, C. Y., and Wu, Y. (2013). Marginal empirical likelihood and sure independence feature screening. Annals of statistics, 41(4).
Chang, J., Tang, C. Y., and Wu, Y. (2016). Local independence feature screening for
nonparametric and semiparametric models by marginal empirical likelihood. Annals of
statistics, 44(2):515.
Curtis, C., Shah, S. P., Chin, S.-F., Turashvili, G., Rueda, O. M., Dunning, M. J., Speed, D.,
Lynch, A. G., Samarajiwa, S., Yuan, Y., et al. (2012). The genomic and transcriptomic
architecture of 2,000 breast tumours reveals novel subgroups. Nature, 486(7403):346–352.
Fan, J., Feng, Y., and Song, R. (2011). Nonparametric independence screening in sparse
ultra-high-dimensional additive models. Journal of the American Statistical Association,
106(494):544–557.
Fan, J. and Li, R. (2001). Variable selection via nonconcave penalized likelihood and its
oracle properties. Journal of the American statistical Association, 96(456):1348–1360.
Fan, J. and Lv, J. (2008). Sure independence screening for ultrahigh dimensional feature space. Journal of the Royal Statistical Society: Series B (Statistical Methodology),
70(5):849–911.
Fan, J. and Lv, J. (2010). A selective overview of variable selection in high dimensional
feature space. Statistica Sinica, 20(1):101.

23

Fan, J., Samworth, R., and Wu, Y. (2009). Ultrahigh dimensional feature selection: beyond
the linear model. Journal of Machine Learning Research, 10(Sep):2013–2038.
Fan, J., Song, R., et al. (2010). Sure independence screening in generalized linear models
with np-dimensionality. The Annals of Statistics, 38(6):3567–3604.
Genovese, C. R., Jin, J., Wasserman, L., and Yao, Z. (2012). A comparison of the lasso and
marginal regression. Journal of Machine Learning Research, 13(Jun):2107–2143.
Huang, J., Breheny, P., and Ma, S. (2012). A selective review of group selection in highdimensional models. Statistical science: a review journal of the Institute of Mathematical
Statistics, 27(4).
Itoh, M., Iwamoto, T., Matsuoka, J., Nogami, T., Motoki, T., Shien, T., Taira, N., Niikura, N., Hayashi, N., Ohtani, S., et al. (2014). Estrogen receptor (er) mrna expression
and molecular subtype distribution in er-negative/progesterone receptor-positive breast
cancers. Breast cancer research and treatment, 143(2):403–409.
Ji, S. and Ye, J. (2009). An accelerated gradient method for trace norm minimization.
In Proceedings of the 26th annual international conference on machine learning, pages
457–464. ACM.
Jiang, J., Li, C., Paul, D., Yang, C., Zhao, H., et al. (2016). On high-dimensional misspecified mixed model analysis in genome-wide association study. The Annals of Statistics,
44(5):2127–2160.
Li, R., Liu, J., and Lou, L. (2017). Variable selection via partial correlation. Statistica
Sinica, 27(3):983.
Li, R., Zhong, W., and Zhu, L. (2012). Feature screening via distance correlation learning.
Journal of the American Statistical Association, 107(499):1129–1139.
Liang, F., Song, Q., and Qiu, P. (2015). An equivalent measure of partial correlation coefficients for high-dimensional gaussian graphical models. Journal of the American Statistical
Association, 110(511):1248–1265.
Liu, Y.-R., Jiang, Y.-Z., Xu, X.-E., Hu, X., Yu, K.-D., and Shao, Z.-M. (2016). Comprehensive transcriptome profiling reveals multigene signatures in triple-negative breast cancer.
Clinical cancer research, 22(7):1653–1662.
Luo, S., Song, R., and Witten, D. (2014). Sure screening for gaussian graphical models.
arXiv preprint arXiv:1407.7819.
Ma, S., Li, R., and Tsai, C.-L. (2017). Variable screening via quantile partial correlation.
Journal of the American Statistical Association, pages 1–14.
Meier, L., Van De Geer, S., and Bühlmann, P. (2008). The group lasso for logistic regression.
Journal of the Royal Statistical Society: Series B (Statistical Methodology), 70(1):53–71.
Nardi, Y., Rinaldo, A., et al. (2008). On the asymptotic properties of the group lasso
estimator for linear models. Electronic Journal of Statistics, 2:605–633.
24

Rajalingam, K., Schreck, R., Rapp, U. R., and Albert, S. (2007). Ras oncogenes and
their downstream targets. Biochimica et Biophysica Acta (BBA)-Molecular Cell Research,
1773(8):1177–1195.
Shao, Q.-M. (1999). A cramér type large deviation result for student’s t-statistic. Journal
of Theoretical Probability, 12(2):385–398.
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal
Statistical Society. Series B (Methodological), pages 267–288.
Tseng, G. C., Ghosh, D., and Feingold, E. (2012). Comprehensive literature review and
statistical considerations for microarray meta-analysis. Nucleic acids research, 40(9):3785–
3799.
Vershynin, R. (2010). Introduction to the non-asymptotic analysis of random matrices. arXiv
preprint arXiv:1011.3027.
Yadav, B. S., Chanana, P., and Jhamb, S. (2015). Biomarkers in triple negative breast
cancer: A review. World journal of clinical oncology, 6(6):252.
Yuan, M. and Lin, Y. (2006). Model selection and estimation in regression with grouped
variables. Journal of the Royal Statistical Society: Series B (Statistical Methodology),
68(1):49–67.
Zhou, J., Liu, J., Narayan, V. A., and Ye, J. (2012). Modeling disease progression via fused
sparse group lasso. In Proceedings of the 18th ACM SIGKDD international conference on
Knowledge discovery and data mining, pages 1095–1103. ACM.
Zhu, L.-P., Li, L., Li, R., and Zhu, L.-X. (2011). Model-free feature screening for ultrahighdimensional data. Journal of the American Statistical Association, 106(496):1464–1475.
Zou, H. (2006). The adaptive lasso and its oracle properties. Journal of the American
statistical association, 101(476):1418–1429.
Zou, H. and Hastie, T. (2005). Regularization and variable selection via the elastic net.
Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67(2):301–
320.

25

