Targeted Random Projection for Prediction

arXiv:1712.02445v1 [] 6 Dec 2017

from High-Dimensional Features
Minerva Mukhopadhyay

David B. Dunson

Abstract
We consider the problem of computationally-efficient prediction from high dimensional and highly correlated predictors in challenging settings where accurate variable selection is effectively impossible. Direct application of penalization or Bayesian
methods implemented with Markov chain Monte Carlo can be computationally daunting and unstable. Hence, some type of dimensionality reduction prior to statistical
analysis is in order. Common solutions include application of screening algorithms to
reduce the regressors, or dimension reduction using projections of the design matrix.
The former approach can be highly sensitive to threshold choice in finite samples,
while the later can have poor performance in very high-dimensional settings. We
propose a TArgeted Random Projection (TARP) approach that combines positive
aspects of both strategies to boost performance. In particular, we propose to use information from independent screening to order the inclusion probabilities of the features in the projection matrix used for dimension reduction, leading to data-informed

1

sparsity. We provide theoretical support for a Bayesian predictive algorithm based
on TARP, including both statistical and computational complexity guarantees. Examples for simulated and real data applications illustrate gains relative to a variety
of competitors.

Some key words: Bayesian; Dimension reduction; High-dimensional; Large p, small n;
Random projection; Screening.
Short title: Targeted Random Projection

1

Introduction

In many applications, the focus is on prediction of a response variable y given a massivedimensional vector of predictors x = (x1 , x2 , . . . xp )0 . Often enormous numbers of possibly
collinear predictors x are collected, and the sample size n is modest relative to p, so that
p  n. In such situations, it is common to assume that x can be replaced by a much
lower-dimensional feature vector comprised of sparse linear combinations of the original
predictors. However, accurate learning of the precise lower-dimensional structure is often
not possible, as the data simply do not contain sufficient information even putting aside
the intractable computational problem.
There is a large literature on variable selection in p  n settings. Much of the focus has
been on penalized optimization-based approaches, with some popular methods including
LASSO (Tibshirani, 1996), SCAD (Fan and Li, 2001)), elastic net (Zou and Hastie, 2005),
the Dantzig selector (Candes and Tao, 2007), and MCP (Zhang et al., 2010). There is
2

also a large literature on Bayesian approaches that attempt to characterize uncertainty
in variable selection. Most approaches use some variation on the spike and slab prior
(e.g. Ishwaran and Rao (2005)) or continuous shrinkage priors concentrated at zero with
heavy tails (e.g. Carvalho et al. (2009)). There is an increasingly rich theoretical literature
providing asymptotic support for such methods in the p → ∞ as n → ∞ setting. However,
positive results rely on strong assumptions in terms of a high signal-to-noise ratio and low
linear dependence in the columns of the design matrix.
We are interested in settings where practical performance of the above methods is poor,
due to a combination of statistical and computational intractability. In such settings, it
is common to use variable screening as a pre-processing step. In particular, independent
screening tests for association between y and each xj separately, and selects predictors with
the largest or most significant associations for second stage analysis. In general, screening
can be guaranteed asymptotically to select a superset of the ‘true’ predictors (Fan et al.,
2009). When the number of predictors is sufficiently reduced, one can apply a simple
maximum likelihood approach, penalized optimization, or Bayesian Markov chain Monte
Carlo (MCMC) algorithms in the second stage. However, when the predictors are highly
correlated and/or the true data generating process does not exhibit strong sparsity with
a high signal-to-noise ratio, it may be necessary to use a very conservative threshold for
the measure of marginal association, limiting the dimensionality reduction occurring in the
first stage.
As an alternative to variable screening, there is a rich literature on using random projections (RPs) to reduce data dimensionality prior to statistical analysis. For example,
3

compressive sensing uses RPs to reduce storage and communication costs in signal processing. By exploiting sparsity, the original signal can be reconstructed from the compressive
measurements with high accuracy (see, e.g., Candes and Tao (2005), Donoho (2006), Davenport et al. (2010)). Usual compressive sensing acts in a row-wise manner, reducing the
dimensionality of the design matrix from n × p to m × p, with m  n. This does not
solve the big p problem. There is a relatively smaller literature on column-wise compression, which instead reduces the design matrix from n × p to n × m, with m  p, while
providing bounds on predictive errors (see, e.g., Maillard and Munos (2009), Fard et al.
(2012), Kabán (2014), Thanei et al. (2017), Guhaniyogi and Dunson (2015), Pettenuzzo
et al. (2016)). Guhaniyogi and Dunson (2015) concentrate on approximating predictive
distributions in Bayesian regression. The above literature on RPs focuses primarily on
random matrices with i.i.d elements.
When predictors are very high-dimensional, existing RP methods can fail as they tend
to include many unimportant predictors in each linear combination, diluting the signal.
Potentially, one can attempt to improve performance by estimating the projection matrix,
but this results in a daunting computational and statistical problem. Alternatively, we
propose a TArgeted Random Projection (TARP) approach, which includes predictors in
the RP matrix with probability proportional to their marginal utilities. These utilities are
estimated quickly in a first stage using an independent screening-type approach. To reduce
sensitivity of the results to the different realizations of the RP matrices, we aggregate over
multiple realizations. TARP can be viewed as a type of rapid preconditioning, enabling
improved predictive performance in high-dimensional settings. Compared with applying
4

RPs after screening out predictors, TARP has the advantage of removing sensitivity to
threshold choice by using a soft probabilistic approach.
In Section 2, we propose the methodology including the computational algorithm, choice
of different tuning parameters and an analysis of computational complexity. We focus on
generalized linear models (GLMs) for ease in presentation and development of a strong
theoretical justification, although TARP can be applied directly in general settings to provide lower-dimensional features that can then be used in any predictive algorithm (random
forests, Gaussian processes, neural nets, etc). Section 3 provides theory on convergence
rates for the predictive distribution of y. Section 4 contains a simulation study comparing
performance with a variety of competitors. In Section 5, we apply TARP to a variety of
real data applications including a genomics dataset with millions of predictors. Section 6
contains a brief discussion, and proofs are included in an Appendix.

2

The Proposed Method

Let Dn = {(yn ; Xn ) : yn ∈ Rn , Xn ∈ Rn×pn } denote the dataset consisting of n observations
on pn predictors x1 , x2 , . . . , xpn and a response y, and (yi ; xi ) denote the ith data point,
i = 1, 2, . . . , n. Suppose that the data can be characterized by a generalized linear model
(GLM). The density of y is related to the predictors as

1
0
0
{yi a(xi β) + b(xi β) + c(yi )} ,
f (yi |β, σ ) = exp
d(σ 2 )
2



5

(1)

where a(·) and b(·) are continuously differentiable functions, a(·) has non-zero derivative
and d(·) is a non-zero function. The vector of regression coefficients β ∈ Rpn , and σ 2 ∈ R+ is
the scale parameter. We approximate the density of y in a compressed regression framework
as follows:

1
0
0
{yi a ((Rn xi ) θ) + b ((Rn xi ) θ) + c(yi )} .
f (yi |θ, Rn , σ ) = exp
d(σ 2 )
2



Here Rn ∈ Rmn ×pn is a random projection matrix, θ ∈ Rmn is the vector of compressed
regression coefficients, and mn  pn . We discuss the choice of the random matrix Rn in
Section 2.1, and illustrate the method in detail in Section 2.2.
Priors. We assume that the covariates are standardized. Taking a Bayesian approach
to inference, we assign priors to θ and σ 2 . The vector of compressed regression parameters
θ is assigned a Nmn (0, σ 2 I) prior given σ 2 , where 0 is a vector of zeros, and I is the identity
matrix. The scale parameter σ 2 is assigned an Inv-Gamma (aσ , bσ ) prior, with aσ , bσ > 0.
The Normal-Inverse Gamma (ING) prior is a common choice of prior for GLMs. In the
special case of a Gaussian likelihood, this prior is conjugate, and the posterior and predictive
distributions are available in analytic forms.

2.1

Choice of the projection matrix

The projection matrix Rn embeds Xn to a lower dimensional subspace. If pn is not large,
the best linear embedding can be estimated using a singular value decomposition (SVD)
of Xn . The projection of Xn to the space spanned by the singular vectors associated with

6

the first mn singular values is guaranteed to be closer to Xn in an appropriate sense than
any other mn -dimensional matrix. However, if pn is very large with pn  n, then it is
problematic to estimate the projection, both computational and statistically, and random
projection (RP) provides a practical alternative. If an appropriate RP matrix is chosen, due
to Johnson-Lindenstrauss (JL) type embedding results, distances between sample points
are maintained (see Dasgupta and Gupta (2003), Achlioptas (2003)).
Our focus is on modifying current approaches by constructing RP matrices that incorporate sparsity in a way that the predictors xj having relatively weak marginal relationships
with y are less likely to be included in the matrix. In particular, the TArgeted Random
Projection (TARP) matrices are constructed as follows:
i.i.d.

γ = (γ1 , γ2 , . . . , γpn )0
qj ∝ |rxj ,y |δ
Rγ = Omn ×(pn −pγ )

and γj ∼ Bernoulli (qj )
for some constant δ > 0,

and

where
(2)

Rγ = Rn∗ ,

where rxj ,y is the Pearson’s correlation coefficient of xj and yn , qj is the inclusion probability
of xj , Rγ and Rγ are the sub-matrices of Rn with columns corresponding to non-zero and
P
zero values of γ, respectively, and Rn∗ is the mn × pγ projection matrix where pγ = j γj .

We prioritize predictors based on their marginal utilities, q = (q1 , q2 , . . . , qpn )0 , and

consider a random subset of the predictors with inclusion probabilities q. This can be
viewed as a randomized version of independent screening. The selected subset is further
projected to a lower dimensional sub-space using Rn∗ . There are many possible choices of Rn∗
7

which can successfully reduce the dimension of the selected variables while having minimal
impact on prediction accuracy. Two predominant classes of such projection matrices are
based on partial SVD and random projections facilitating JL type embedding. We consider
both these choices, as described below.
∗
, of Rn∗ is sampled independently from a three
Random projection. Each element, Rk,j

point distribution as

∗
Rk,j



 ±1/√2ψ
=

 0

with probability ψ,

(3)

with probability 1 − 2ψ,

where ψ ∈ (0, 0.5) is a constant.

Projection matrices of this form are widely used due to their inter point distance preservation property. Incorporating zero values facilitates data reduction and improves computational efficiency. We refer to the method that generates Rn∗ in (2) from (3) as RIS-RP
(Randomized Independent Screening-Random Projection).
Remark 1. The choice of projection matrix in (3) can be replaced by a wide variety of
matrices having i.i.d. components with mean zero and finite fourth moments. One of the
√
∗
sparsest choices is of the form Rk,j
= ±nκ/2 / mn with probability 1/2nκ , 0 with probability
(1 − 1/nκ ), where mn ∼ nκ (see Li et al. (2006)). This choice of projection matrix reduces
the data to a great extent, and is useful in compressing extremely large dimensional data.
Our theoretical results would hold also if we consider a random matrix Rn∗ as above.
Principal component projection. Another alternative is to use the matrix of principal
8

component scores. Let Xγ be the sub-matrix of Xn with columns corresponding to non-zero
0
values of γ. Consider the partial spectral decomposition Xγ0 Xγ = Vγ,m
Dγ,mn Vγ,mn . The
n

projection matrix Rn∗ we consider is
0
.
Rn∗ = Vγ,m
n

(4)

The method corresponding to this choice of projection matrix combines a randomized
version of independence screening with principal component regression (PCR). Therefore,
we refer to this method as RIS-PCR (Randomized Independence Screening-PCR).
The performance of TARP depends on tuning parameters mn , δ and ψ. In addition, for
any given choice of tuning parameters, different realizations of the random projection matrix will vary, leading to some corresponding variation in the results. To limit dependence
of the results on the choice of tuning parameters and random variation in the projection,
we take the approach of generating multiple realizations of the matrix for different choices
of tuning parameters and aggregating these results together. Potentially, one could estimate weights for aggregation using Bayesian methods (see Hoeting et al. (1999)) or other
ensemble learning approaches, but we focus on simple averaging due to its computational
and conceptual simplicity.

2.2

Posteriors and Predictive Distribution

We illustrate the proposed method in normal linear regression for simplicity, although it is
applicable to more general settings. In the compressed regression framework, we replace
9

the normal linear model yi = x0i β + ei by yi = (Rn xi )0 θ + ei , where ei ∼ N (0, σ 2 ).
Given the NIG prior structure stated above, the posterior distribution of θ follows a scaled
mn -variate t distribution with degrees of freedom (d.f.) n + 2aσ , location vector µt , and
scale matrix Σt , where µt = Z(Xn Rn0 )0 yn , Σt = (yn0 yn − µ0t Z −1 µt + 2bσ )Z/(n + 2aσ ) and
Z = (Rn Xn0 Xn Rn0 + I)−1 . Moreover, the posterior distribution of σ 2 , given Dn and Rn , is
inverse gamma with parameters aσ + n/2 and (yn0 yn − µ0t Σ−1
t µt + bσ )/2.
Consider the problem of point prediction of y when nnew new data points on the predictors are obtained, given the dataset Dn . The predicted values of y, say ynew , can be
obtained using the Bayes estimator of θ under squared error loss as follows
ŷnew = Xnew Rn0 θ̂ Bayes where θ̂ Bayes = (Zn0 Zn + I)−1 Zn0 yn and Zn = Xn Rn .
Here Xnew is the new design matrix. Moreover, the posterior predictive distribution of ynew
is a nnew -variate t distribution with degrees of freedom n + 2aσ , location vector ŷnew and
0
scale parameter (yn0 yn − µ0t Z −1 µt + 2bσ )(I + Xnew ZXnew
)/(n + 2aσ ).

When the distribution yn is non-normal, analytical expressions of the posteriors of θ,
σ 2 and predictive distribution of ynew are not available. In such cases, it is common to rely
on a Laplace approximation (see Tierney and Kadane (1986)) or sampling algorithm, such
as MCMC. In the compressed regression framework, as the pn -dimensional variables are
projected to a much lower-dimensional mn -hyperplane, with mn  pn , we are no longer in
a high-dimensional setting. Hence, MCMC is computationally tractable.

2.3

Tuning Parameter Choice

Next, we describe the choices of tuning parameters that are involved in TARP.
10

Choice of mn . The parameter mn determines the number of linear combinations of
predictors we consider. Instead of choosing a fixed value of mn , we consider a range of
values. In particular, we suggest choosing values in the range (2 log pn , min{3n/4, pn }),
consistent with our theoretical results in Section 3 requiring that mn < n and with our
experiments assessing performance for different choices in a variety of scenarios.
Choice of δ and q. The parameter δ plays an important role in screening. Higher values
of δ lead to fewer variables selected in the screening step. If pn  n, one would like to select
a small proportion of predictors, and if pn ∼ n then selection of a moderate proportion is
desirable. We recommend δ = max{0, (1+log(pn /n))/2} as a default. This function selects
all the predictors if pn  n, and becomes more restrictive as pn becomes larger than n. The
selection probabilities in the RIS stage are then qj = |rxj ,y |δ / maxj |rxj ,y |δ , j = 1, 2, . . . , pn .
Hence, the variable with highest marginal correlation is definitely included in the model.
Choice of ψ. The value of ψ controls sparsity in the random matrix, and it is necessary
to let ψ ∈ (0, 0.5). Achlioptas (2003) suggest choosing ψ = 1/6 as a default value. To avoid
sensitivity of the results for a particular choice of ψ, we choose ψ ∈ (0.1, 0.4) avoiding the
very sparse and dense cases.

2.4

Computational Algorithm and Complexity

We now illustrate the time complexity of TARP along with the algorithm for computation
in normal linear models.
RIS-RP. For a specific choice of (mn , δ, ψ), calculation of ŷnew using RIS-RP involves the
following steps:
11

1: Calculate rxj ,y for j = 1, . . . , pn .
2: Generate γj ∼ Bernoulli(qj ) where qj = |rxj ,y |δ / max{|rxj ,y |δ }, j = 1, . . . , pn .
IF γj = 1, generate Rn with Ri,j as in (3).

ELSE set Ri,j = 0.

3: Post-multiply Rn with Xn . Set Zn = Xn Rn .
4: For a given Xnew , compute Znew = Xnew Rn and ŷnew = Znew θ̂.
The complexity of steps 1, 2-3, 4 and 5 are O(pn ), O(npγ mn ), O(nm2n ) and O (nnew pγ mn ),
P
respectively, where pγ =
γj . Thus, if nnew ≤ n, the total complexity for a single choice
of (mn , δ, ψ) is O(pn ) + 2O(nmn pγ ) + O(nm2n ) without using parallelization.

RIS-PCR. RIS-PCR differs from RIS-RP in step 2 of the algorithm. After generation
of γ RIS-PCR requires SVD of Xγ involving complexity O (npγ min{n, pγ }). Thus, total
complexity of RIS-PCR for a single choice of (mn , δ, ψ) can similarly be derived. Therefore,
the two methods have comparable time complexity unless either n or pγ is much larger
than mn . Although theoretically we do not impose any restriction on pγ , in practice when
pn = exp{o(n)} and δ ≥ 2, pγ is usually of order n.
Increment of complexity due to aggregation. Suppose N different choices of (mn , ψ, Rn )
are considered. Each choice yields a model Ml : y ∼ f (y|x, mn,l , ψl , Rn,l ) along with a corresponding estimate of ynew (say ŷnew,l ), where l ∈ {1, 2, . . . , N }. The proposed estimate
is the simple average of these N estimates of ynew .
Step 1 in the TARP algorithm is not repeated over the aggregration replicates, while
the remaining steps are repeated N times. In addition, the first step of screening and
12

aggregration are embarassingly parallelizable. Hence, given k processors, if nnew ≤ n, the
total complexity is O(pn /k) + 2O(N nmn pγ /k) + O(N nm2n /k) for RIS-RP and O(pn /k) +
O(N npγ min{n, pγ }/k) + O(N nmn pγ /k) + O(N nm2n /k) for RIS-PCR.

3

Theory on Predictive Accuracy

We study the asymptotic performance of the predictive distribution produced by TARP
for a single random projection matrix without considering the aggregation step. We focus
on weakly sparse and dense cases where the absolute sum of the true regression coefficients
is bounded. This condition also includes strong sparsity where only a few covariates have
non-zero coefficients.
The projection matrix in TARP depends on the random variable γ, and therefore is
denoted by Rγ . We denote a particular realization of the response variable as y, and a
particular realization of the variables (x1 , x2 , . . . , xpn )0 as x. Let f0 be the true density of
y given the predictors, and f (y|x, γ, Rγ , θ) be the conditional density of y given the model
induced by (γ, Rγ ), the corresponding vector of regression coefficients θ and x. We follow
Jiang (2007) in showing that the predictive density under our procedure is close to the true
predictive density in an appropriate sense.
We assume that each covariate xj is standardized so that |xj | < M , for j = 1, 2, . . . , pn ,
with M a constant. We also assume that the scale parameter σ 2 in (1) is known. We
require the following two assumptions on the design matrix.
Assumption (A1) Let rxj ,y denote the correlation coefficient of the observed values
13

of xj and y, 1 ≤ j ≤ pn . Then for each data point (y, x) and constant δ in (2), there exists
a positive constant αδ such that
pn
1 X 2
xj |rxj ,y |δ → αδ .
lim
pn →∞ pn
j=1

Assumption (A2) Let q(γ) =

Qn

γj
(1−γj )
,
i=1 qj (1 − qj )

with qj defined in (2), denote the

probability of obtaining a particular γ = (γ1 , . . . , γpn )0 in the random screening step. Let
Pn
γj = l, and let Ml ⊂ Γl
Γl ⊂ {0, 1}pn denote the set of γ vectors such that pγ = pj=1

denote the first pknn elements of Γl ordered in their q(γ) values. Let An denote the event

that γ ∈ Ml for some l. Then, P (Acn ) = P {γ : γ ∈
/ ∪l Ml } ≤ exp(−nε2n /4), for some
increasing sequence of integers {kn } and sequence {εn } satisfying 0 < ε2n < 1 and nε2n → ∞.

Remark 2. As the probability of selection in the random screening step depends on the
empirical correlation between the predictor and the response, assumption (A2) is on the
data generating process. If l ≤ kn or l ≥ pn − kn , then all models of dimension l belong to
Ml , and none of the corresponding γ vectors belong to Acn . If kn < l < pn − kn , there are
more than pknn models of dimension l, but as the models are ordered in terms of their q(γ)
values, the models not falling in Ml should have extremely small probabilities of selection,
hence satisfying (A2). Violations of (A2) would imply that large numbers of predictors have
empirical correlations that are not close to zero.
Measure of closeness: Let νx (dx) be the probability measure for x, and νy (dy) be the
dominating measure for conditional densities f and f0 . The dominating measure of (y, x)
is taken to be the product of νy (dy)νx (dx).
14

The Hellinger distance between f and f0 is given by
rZ 
p
p 2
f − f0 νx (dx)νy (dy).
d(f, f0 ) =

The Kullback-Leibler divergence between f and f0 is given by
 
Z
f0
νx (dx)νy (dy).
d0 (f, f0 ) = f0 ln
f
 t
Z

f0
−1
Define dt (f, f0 ) = t
f0
νx (dx)νy (dy) − 1 , for any t > 0.
f

Consider the following two facts: (i) d(f, f0 ) ≤ (d0 (f, f0 ))1/2 , and (ii) dt (f, f0 ) decreases
to d0 (f, f0 ) as t decreases to 0 (see Jiang (2007)).
Let Pn be a sequence of sets of probability densities, and εn be a sequence of positive
numbers. Let N (εn , Pn ) be the εn -covering number, i.e., the minimal number of Hellinger
balls of radius εn needed to cover Pn .
RIS-RP. The result showing asymptotic accuracy in approximating the predictive density using RIS-RP is stated below.
Theorem 1. Let θ ∼ N (0, σθ2 I), and f (y|x, γ, Rγ , θ) be the conditional density of y given
the model induced by (γ, Rγ ), where Rγ is as in (2) and (3). Let β 0 be the true regression
P
parameter with
j |β0,j | < K for some constant K, and assumptions (A1)-(A2) hold.

Consider the sequence {εn } as in assumption (A2) satisfying 0 < ε2n < 1 and nε2n → ∞,
and assume that the following statements hold for sufficiently large n:

(i) mn | log ε2n | < nε2n /4,
(ii) kn log pn < nε2n /4, and

 p

< nε2n /4, where
(iii) mn log 1 + D σθ 6nε2n pn mn
15

D(h∗ ) = h∗ suph≤h∗ |a0 (h)| suph≤h∗ |a0 (h)/b0 (h)|, b(·) as in (1). Then,
h
i
2
2
Pf0 π {d(f, f0 ) > 4εn |Dn } > 2e−nεn /4 ≤ 2e−nεn /5 ,
where π{·|Dn } is the posterior measure.
The proof of Theorem 1 is given in the Appendix.
Remark 3. If we consider the sparse choice of Rγ , as described in Remark 1, the same
line of proof (as that of Theorem 1) would go though. For both the choices of Rγ , each
component of the projection matrix has expectation zero and finite fourth moment. For the
random matrix in (3), the probability of choosing a non-zero element, P (Ri,j 6= 0), is fixed,
while the probability is decaying with n for the sparse choice of Rγ . However, the rate of
decay is such that distances are preserved between the sample-points, a critical property for
proving consistency.
RIS-PCR. Asymptotic guarantees on predictive approximation for RIS-PCR requires an
additional assumption.
Assumption (A3) Let Xγ be the sub-matrix of Xn with columns corresponding to nonzero values of γ, and xγ be a row of Xγ . Let Vγ be the mn × pγ matrix of mn eigenvectors
corresponding to the first mn eigenvalues of Xγ0 Xγ . Then, for each γ and data point xγ ,
kVγ xγ k2
≥ αn ,
kxγ k2

16

where αn ∼ (nε2n )−1 , where the sequence {ε2n } is as in assumption (A2).
Remark 4. If the matrix Xγ0 Xγ has rank less than mn , then αn = 1 by Perseval’s identity.
Consider the situation where rank of the gram matrix, say rn (≤ n), is bigger than mn . Then
the row space of Xγ , or that of Xγ0 Xγ , is spanned by a set of rn basis vectors v1 , v2 , . . . , vrn .
Therefore, any data point x can be written as a linear combination of these rn vectors as
x = a1 v1 + a2 v2 + · · · + arn vrn , where a1 , a2 , . . . , arn are constants not all equal to zero. As
the vectors vj are orthonormal, vj0 x = aj for all j = 1, 2, . . . , rn , which in turn implies that
Pn 2
x0 x = rj=1
aj . Also, note that the first mn among these rn vectors constitute Vγ0 , which
P n 2
Pmn 2 Prn 2
2
2
implies kVγ0 xk2 = m
j=1 aj . Thus kVγ xk /kxk =
j=1 aj /
j=1 aj , and magnitude of the
ratio depends on the part of x explained by the last few principal component directions.

The lower bound αn ∼ (nε2n )−1 is weaker than many real data scenarios where most of the
variation is explained by the first few principal components.
Theorem 2. Let θ ∼ N (0, σθ2 I), and f (y|x, γ, Rγ , θ) be the conditional density of y given
the model induced by (γ, Rγ ), where Rγ is as in (2) and (4). Let β 0 be the true regression
P
parameter with j |β0,j | < K for some constant K, and assumptions (A1)-(A3) hold. Assume that the conditions (i)-(iii) of Theorem 1 hold for the sequence {εn } as in assumption
(A2) satisfying 0 < ε2n < 1 and nε2n → ∞. Then,
h
i
2
2
Pf0 π {d(f, f0 ) > 4εn |Dn } > 2e−nεn /4 ≤ 2e−nεn /5 ,
where π{·|Dn } is the posterior measure.
The proof of Theorem 2 is given in the Appendix.
17

Remark 5. The conditions (i)-(iii) in Theorems 1 and 2 are related to the sizes of pn ,
mn and kn in comparison with nε2n . A sufficient condition for (i) is mn log n < nε2n /4,
providing an upper bound on the dimension of the subspace mn . Condition (ii) restricts the
permissible number of regressors pn , and the number of possible models of each dimension
kn . If there is a strict ordering in the marginal correlation coefficients |rxj ,y |, so that
kn ≤ κ for some large number κ (see assumption (A2)), then the condition reduces to
log pn < nε2n /4. To illustrate that condition (iii) tends to be weak, consider distributions of
y corresponding to Bernoulli, Poisson and normal. For these cases, the quantity D(h∗ ) is
at most order O(h∗ ). Therefore, condition (iii) does not impose much additional restriction
over (i)-(ii), except mn log pn < nε2n /4, inducing a stronger upper-bound to mn .

4

Simulation Study

In this section, we consider different simulation schemes (Scheme I – IV ) to compare
TARP with a variety of methods. We mainly focus on high-dimensional and weakly sparse
regression problems with a variety of correlation structures in the predictors. The sample
size is taken to be 200, while pn varies. Additional results for different choices of n are
provided in the Supplement.
Competitors. We compare with: SCAD screened by iterative SIS (ISIS), ISIS-SCAD;
minimax concave penalty (MCP) method screened by ISIS, ISIS-MCP ; LASSO screened
by sequential strong rule (SSR, Tibshirani et al. (2012)), SSR-LASSO; ridge regression
screened by SSR, SSR-Ridge; elastic net screened by SSR, SSR-EN ; principal component
18

regression (PCR); sparse PCR, SPCR (see Witten et al. (2009)); robust PCR, RPCR (see
Candès et al. (2011)); and Bayesian compressed regression (BCR). ISIS-SCAD and ISISMCP are available in the ‘SIS’ package, and LASSO, ridge and elastic net are available in
the ‘biglasso’ package (Zeng and Breheny (2017)). SPCR and RPCR are performed using
‘PMA’ and ‘rsvd’ packages in R, respectively. To estimate PC scores, we rely on approximate SVD using fast.svd in the ‘corpcor’ package. For BCR, we average over 100 different
random projections with varying mn values within the range [2 log pn , 3n/4]. We use the
qr function in R to apply QR factorization in place of Gram-Schmidt orthogonalization of
the random matrix, which is computationally prohibitive for large pn .
The proposed method. We select the tuning parameters of TARP as described in Section
2.3. The parameter mn is chosen in the range [2 log pn , 3n/4]. We assign δ = 2 as the
function max{0, (1 + log(pn /n))/2} is close to 2 for all the choices of (n, pn ). Further, the
hyperparameters of the inverse gamma priors are set to 0.02 to correspond to a minimally
informative prior.
Simulation Schemes. In the first three simulation schemes, the predictors were generated
from N (0, Σ), with different choices of pn , Σ and the regression coefficients. In Scheme IV
we consider a functional regression setup. Different methods are compared with respect to
their performance in out of sample prediction. We calculate mean square prediction error
(MSPE), empirical coverage probability (ECP) of a 50% prediction interval (PI) and the
width of the PI for each of a 100 replicate datasets in each simulation case.
Scheme I: First order autoregressive structure. Σi,j = (0.9)|i−j| , i, j = 1, . . . , pn , with

19

Figure 1: Box-plots of MSPE for pn = 2000 in Scheme I.
pn ∈ {2000, 3000}, and βj = 0 for all but a randomly selected set of 30 predictors having
βj = 1.
Scheme II: Block diagonal covariance structure. We choose (pn /100 − 2) blocks of 100
predictors each, along with 200 independent predictors, with pn ∈ {104 , 2 × 104 }. The
within-block correlation is ρ and the across-block correlation is zero, with ρ = 0.3 for half
of the blocks and ρ = 0.9 for the remaining half. There are 21 non-zero βj s having βj = 1,
with 20 of the corresponding predictors in the ρ = 0.9 blocks and the remaining in the
independent block.
Scheme III: Principal Component Regression. We first choose a matrix P with orthonormal
columns, and a 3 × 3 matrix D = diag(152 , 102 , 72 ). We set Σ = P DP 0 and choose β = P·,1 ,
where P·,1 is the first column of P . This produces an Xn with three dominant principal
components, with the response yn dependent on the first and pn ∈ {104 , 5 × 104 }.
Scheme IV: Functional Regression. Finally, we consider a functional regression setup, where
the covariates are generated from Brownian bridge Bt with t ∈ (0, 5) and values ranging
from (0, 10). A set of 20 covariates is randomly selected as active, each having regression
20

Figure 2: Box-plots of MSPE for pn = 10000 in Scheme II.

Figure 3: Box-plot of MSPEs for pn = 5000 in Scheme III.
parameters in the range (2, 2.5), and pn ∈ {104 , 2 × 104 }.
Results. For each of the simulation schemes, we present the results corresponding to the
first choice of pn , with the other results provided in the Supplement.
In Scheme I (see Figure 1) SSR-Ridge, PCR, RIS-RP and RIS-PCR show competitive
performance, with RIS-PCR showing the best overall result. Performance of RIS-PCR is
closely followed by RIS-RP, which in turn is followed by PCR and SSR-Ridge. Apart from
these four methods, ISIS based approaches and SSR-EN exhibit reasonable performance.
Although ISIS based methods have lower average MSPE, the variances of MSPEs are high
indicating less stability. SPCR, RPCR, BCR and SSR-LASSO fail to perform adequately.
21

(a) MSPE of all the methods.

(b) MSPE of selected methods.

Figure 4: Box-plot of MSPEs for pn = 10000 in Scheme IV.
For Scheme II, PCR, RIS-RP and RIS-PCR again have the best performance (see Figure
2). RIS-PCR yields lowest overall MSPE, closely followed by RIS-RP and PCR. SSR-EN,
SSR-Ridge and RPCR exhibit moderate performance. Performance of these three methods
are followed by LASSO, SPCR and BCR. Although SPCR enjoys better average MSPE, it
yields higher dispersion as well. ISIS based methods fail to perform adequately.
In Scheme III (see Figure 3) SSR based methods fail to exhibit competitive performance.
Among the SSR based methods, performance of SSR-LASSO is better than the others. All
the other methods perform well, although ISIS based methods occasionally show higher
MSPEs.
In Scheme IV (see Figure 4(a)), SPCR and ISIS based methods fail completely, making
the other methods indistinguishable in the figure. Hence, we separately show these methods
in Figure 4(b). Among the other methods, PCR, RIS-RP and RIS-PCR have the best
overall performance closely followed by BCR and then by SSR-EN. The other three methods
fail to exhibit competitive performance.

22

We next consider the empirical coverage probabilities (ECPs) of 50% prediction intervals
(PI), and the width of PIs. For the Bayesian methods, the PI is obtained from the highest
credible region of the predictive distribution of ynew given Dn , Xnew . The PIs for the freq

quentist methods can be obtained as ynew ±tn−pγ −1, α2 M SP E 1 + x0γ,new (Xγ0 Xγ )−1 xγ,new ,

where tn−pγ −1, α2 is the upper α/2 point of t distribution with (n − pγ − 1) degrees of freedom, the suffix γ indicates consideration of the regressors selected by the corresponding

method, and pγ is the number of selected regressors. The PIs for PCR based methods
can be obtained similarly, except xγ,new and Xγ are replaced by the principal component
scores. For SSR-Ridge and SSR-EN, pγ is much larger than n. Therefore, the PIs of these
methods could not be calculated using the above formula, instead we obtain the interval as
p
√
ynew ± tn, α2 M SP E + se(ŷ|Dn , xnew )2 , or ynew ± tn, α2 2M SP E, whichever gives better re-

sults, where se(ŷ|Dn , xnew )2 is the variance of the fitted values. The results are summarized
in Table 1.
Summary. Considering all the simulation schemes, RIS-RP, RPCR, SPCR and SSRLASSO have the best performance with respect to mean ECP, with RIS-RP having lowest
average width among these four methods. The average width of PI is much bigger for
SPCR, particularly in Scheme IV. ISIS based methods and RIS-PCR have relatively lower
coverage probabilities in general, although among these methods, RIS-PCR has higher coverage with much lower average width than others, especially in Scheme IV. The average
ECP for PCR is satisfactory, although the corresponding widths of PIs have large variances
in almost all the simulation schemes. This indicates instability in the overall performance

23

Table 1: Mean and standard deviation(sd) of empirical coverage probabilities of 50% prediction intervals, and mean and sd of width of 50% prediction intervals.
Methods→
Scheme, pn

ISISSACD

ISISSACD

SSRSSRLASSO EN

SSRRidge

PCR

SPCR

RPCR BCR

RISRP

RISPCR

Average and standard deviation (in braces) of empirical coverage probability
I

2 × 103

0.309
(.052)

0.311
(.058)

0.432
(.051)

0.549
(.040)

0.692
(.055)

0.498
(.051)

0.457
(.055)

0.433
(.290)

0.286
(.053)

0.493
(.059)

0.431
(.057)

II

104

0.327
(.056)

0.324
(.055)

0.429
(.049)

0.607
(.062)

0.702
(.075)

0.455
(.308)

0.499
(.058)

0.445
(.050)

0.278
(.049)

0.502
(.056)

0.362
(.050)

III

5 × 103

0.494
(.053)

0.494
(.053)

0.503
(.058)

0.657
(.031)

0.678
(.127)

0.503
(.058)

0.494
(.049)

0.494
(.049)

0.527
(.054)

0.494
(.055)

0.507
(.054)

0.425
(.053)

0.416
(.056)

0.477
(.065)

0.660
(.031)

0.665
0.448
(0.102) (.300)

0.488
(.048)

0.491
(.053)

0.487
(.058)

0.696
(.052)

0.418
(.055)

IV

104

Average and standard deviation (in braces) of width of the 50% prediction interval
I

2 × 103

3.894
(.500)

3.857
(.480)

5.939
(.424)

7.410
(.578)

8.192
(.528)

7.739
7.964
(8.309) (.700)

6.483
(.451)

5.260
(.416)

5.029
(.194)

4.253
(.365)

II

104

5.868
(.657)

5.836
(.684)

5.113
(.527)

7.335
(.615)

8.139
(.830)

8.407
6.066
4.635
(25.713) (1.103) (.407)

4.894
(.433)

4.204
(.198)

2.744
(.232)

III

5 × 103

1.424
(.174)

1.385
(.144)

1.362
(.072)

1.945
(.145)

5.755
1.362
(1.928) (.072)

1.366
(.069)

1.463
(.073)

1.351
(.079)

1.391
(.077)

13.63
13.34
2.170
(2.705) (2.540) (.422)

1.537
(.139)

3.582
(.501)

2.508
26.512 2.480
(3.169) (6.209) (.345)

1.792
(.115)

2.284
(.268)

1.154
(.123)

IV

104

1.366
(.069)

of PCR. BCR shows under-coverage in the first two simulation schemes, but performs well
with respect to both the measures in Schemes III and IV. Finally the other two methods,
viz., SSR-Ridge and SSR-EN have higher values of ECP, along with higher width of PIs.
SSR-Ridge has highest average width of PI in Schemes I and III. In all the simulation
schemes SSR-EN outperforms SSR-Ridge with respect to width of PI.

4.1

Computational Time

The computational time of a method may depend on the simulation scheme due to varying
level of complexity in the dataset. We only present the computational time for Scheme
24

IV as an example. Figure 5 presents the time (in minutes) taken by different methods to
compute ŷnew using a single core, as pn grows and n = nnew = 100. We run all the methods
in R 3.4.1 in a 64 bit Dell-Inspiron desktop with Ubuntu 16.04 LTS operating system, 15.6
GB random access memory and Intel R CoreTM i5-4460 CPU @ 3.20GHz processor.
Results. When pn is below 104 all the methods require comparable computational
time. When pn is increased, SSR based methods, except SSR-Ridge, RPCR, and fast.svd
based PCR continue to require low computational time. ISIS-SCAD, SPCR and RIS-RP
also have reasonable computational expense (approximately, 5 minutes for pn = 5 × 105 ).
Computational time of BCR, ISIS-MCP and RIS-PCR tends to increase rapidly after pn =
105 . Among these three methods, RIS-PCR requires highest system time (approximately
27 minutes for pn = 5 × 105 ). The computational time required by SSR-Ridge exceeds all
other methods for pn > 5×104 , and for pn = 5×105 it becomes computationally prohibitive
(it takes more than 2 hours).
The increment of computational time of RIS-PCR is due to the computation of exact
SVD of the screened design matrix Xγ . However, this burden would immediately be reduced
if one uses some approximation of the SVD. In that case the computational time would be
comparable to RIS-RP.

5

Real Data Analysis

In this section, we study the performance of TARP using three real datasets, viz., Golub
dataset, GTEx dataset and Eye dataset. The Golub dataset is available at GitHub https:

25

Figure 5: System time required by different methods to predict y as pn grows.

//github.com/ramhiser/datamicroarray/wiki, the GTEx dataset is available at GTEx
portal https://www.gtexportal.org/home/ and the eye data is available in the flare
package of R. In each case, we assess out-of-sample predictive performance averaging over
multiple training-test splits of the data.
Golub data. The Golub data consist of 47 patients with acute lymphoblastic leukemia
(ALL) and 25 patients with acute myeloid leukemia (AML). Each of the 72 (= n) patients
had bone marrow samples obtained at the time of diagnosis (see Golub et al. (1999)).
Expression levels of 7129 (= pn ) genes have been measured for each patient. We consider
a training set of size 60 with 20 AML patients, and 40 ALL patients. The test set consists
of the remaining 12 samples.
GTEx Data. To understand the functional consequences of genetic variation, Consortium
26

et al. (2015) presented an analysis of RNA sequencing data from 1641 samples across 43
tissues from 175 individuals, generated as part of the pilot phase of the Genotype-Tissue
Expression (GTEx) project. We selected RNA-seq data on two normal tissues, viz., ArteryAorta and Artery-Tibial. The dataset contains RNA-seq expressions on 36115 (= pn ) genes
and 556 (= n) samples, among which 224 are from Artery-Aorta, and 332 are from ArteryTibial. A training set of 100 samples from each of the tissue types is considered, and the
remaining 446 samples are used as test set.
Eye Data. The Eye dataset consists of gene expressions for 200 (= pn ) gene probes
from the microarray experiments of mammalian-eye tissue samples of 120 (= n) rats (see
Scheetz et al. (2006)). The response variable is the expression level of the TRIM32 gene.
We consider 100 sample points as the training set, and the remaining 20 samples as the
test set.
Golub and GTEx datasets have nominal response, and therefore the methods are evaluated by the misclassification rate (in %) and the area under the receiver operating characteristic (ROC) curve. Table 2 provides the average and standard deviation (sd) of percentages
of misclassification, and those for the area under the ROC curve over 100 random subsets
of the same size chosen from the dataset for the competing methods. We further compare
the predictive performance of the methods in terms of mean squared difference of predictive and empirical probabilities for these two datasets. Most methods (except SSR based
methods and SPCR) exhibit similar performance in this aspect. We provide the details of
the predictive calibration in the Supplement.
The eye dataset has continuous response, and therefore we evaluate the methods by
27

Table 2: Mean and standard deviation (in braces) of percentage of misclassification and
area under the ROC curve for Golub and GTEx datasets, and those of MSPE, ECP of 50%
PI and width of PI for Eye dataset for all the competing methods.
Methods→
Dataset ↓

ISISSACD

ISISSACD

SSRSSRLASSO EN

SSRRidge

PCR

SPCR

RPCR BCR

RISRP

RISPCR

Misclassification rate and Area under ROC curve for Datasets with Categorical response
Mean and SD of Misclassification Rate (in %)
Golub

11.82
(6.90)

11.50
(7.06)

45.45
(0.00)

45.45
(0.00)

45.45
(0.00)

7.09
(5.68)

41.36
9.73
(13.31) (7.28)

19.36
(9.79)

5.54
(4.36)

5.77
(4.52)

GTEx

0.00
(0.00)

0.00
(0.00)

34.83
(0.00)

34.83
(0.00)

34.83
(0.00)

0.06
(0.13)

3.53
(3.31)

13.28
(3.79)

0.39
(0.20)

0.49
(0.32)

Golub

0.876
(.073)

0.879
(.074)

0.500
(.000)

0.500
(.000)

0.500
(.000)

0.923
(.062)

0.582
(.134)

0.895
(.078)

0.816
(.093)

0.978
(.027)

0.943
(.044)

GTEx

1.00
(.000)

1.00
(.000)

0.500
(.000)

0.500
(.000)

0.500
(.000)

0.999
(.001)

0.964
(.033)

0.998
(.001)

0.877
(.041)

1.00
(.000)

.996
(.002)

10.01
(4.04)

8.54
(3.09)

8.29
(2.99)

0.22
(0.18)

Mean and SD of Area under Receiver Operating Characteristic curve

MSPE, ECP of 50% PI and Width of 50% PI
Mean and SD of Mean Square Prediction Error
Eye

11.66
(4.06)

11.66
(4.06)

20.92
20.92
7.31
(19.33) (19.33) (2.91)

13.84
(3.94)

8.65
(3.08)

7.67
(3.30)

Mean and SD of Empirical Coverage Probability and Width of the Prediction Interval
Empirical
coverage

0.502
(.138)

0.502
(.138)

0.634
(.130)

0.709
(.106)

0.700
(.076)

0.423
(.325)

0.508
(.123)

0.522
(.114)

0.564
(.117)

0.598
(.101)

0.507
(.107)

Width of
interval

1.208
(.057)

1.208
(.057)

1.970
(.190)

2.033
(.917)

1.539
(.303)

1.884
1.202
(1.612) (.079)

1.055
(.049)

1.249
(.056)

1.341
(.038)

1.056
(.036)

MSPE and empirical coverage probabilities (ECP) of 50% prediction intervals (PI) as in
Section 4. As variation in the expression levels of the TRIM32 gene is very small (the
range is 1.37) we multiply the MSPEs of different methods by 10 to increase the variability.
Table 2 provides the mean and sd of MSPEs, ECPs of 50% PIs, and widths of the PIs over
100 different training and test sets selected from the dataset, for the competing methods.

28

Results. For the Golub data set, both the lowest misclassification rate and the highest area
under ROC curve are achieved by RIS-RP, which is closely followed by RIS-PCR. TARP
based methods attain lower sd than other methods as well. PCR also yields reasonable
performance with 7% average misclassification rate and area under ROC more than 0.9.
RPCR and ISIS-based methods produce average rates of misclassification of at least 10%,
with area under the ROC of ∼ 0.9. BCR possesses high misclassification rate (about 19%),
although area under ROC is more than 0.8. Neither the MSPE, nor the area under ROC
curve, is satisfactory for SPCR. Finally, for all the 100 repetitions, SSR based methods
invariably select the intercept-only model. Thus, the MSPEs of these methods depend
entirely on the proportion of test samples obtained from the two classes.
For the GTEx dataset, perfect classification is achieved by ISIS based methods. These
methods along with RIS-RP have the highest area under the ROC curve. PCR, RPCR,
RIS-RP and RIS-PCR also yield satisfactory results, having less than 0.5% average misclassification rate and more than 99% area under the ROC curve. SPCR is comparable with an
average MSPE of less than 4%. BCR attains 13.3% average misclassification rate, with the
area under the ROC curve almost 0.9. SSR based methods fail to show any discriminatory
power in the GTEx dataset.
SSR-Ridge, RPCR, RIS-PCR, RIS-RP and SPCR yield excellent performance in terms
of MSPE in the eye dataset with an average MSPE of less than 0.9. SSR-Ridge has an
average ECP of about 0.7. RIS-PCR shows more stable performance in terms of ECP,
followed by SPCR. BCR and ISIS based methods have similar overall performance. In
terms of MSPE, BCR outperforms ISIS based methods but is outperformed by ISIS based
29

methods in terms of ECP. PCR is not quite as good in terms of either measure. SSR-LASSO
and SSR-EN again fail to perform adequately for the eye dataset.
The GEUVADIS cis-eQTL dataset

We conclude this section by illustrating the

TARP approach on a massive dataset. The GEUVADIS cis-eQTL dataset (Lappalainen
et al. (2013)) is publicly available at http://www.ebi.ac.uk/Tools/geuvadis-das/. This
dataset consists of messenger RNA and microRNA on lymphoblastoid cell line (LCL) samples from 462 individuals provided by the 1000 Genomes Project along with roughly 38
million SNPs. E2F2 plays a key role in the control of the cell cycle. Hence, as in Chen
and Dunson (2017), we choose the gene E2F2 (Ensemble ID: ENSG00000000003) as the
response. A total of 8.2 million (= pn ) SNPs are pre-selected as candidate predictors on
the basis of having at least 30 non-zero expressions. The total number of subjects included in the dataset is about 450 (= n). The genotype of each SNP is coded as 0, 1 or 2
corresponding to the number of copies of the minor allele.
TARP is applied on this dataset. We consider four different training sample sizes, viz.,
nt = 200, 250, 300 and 350, and test sample size 100 in each case. As pn is huge, we
applied three different values of δ, namely, 2, 5 and 8, to analyze the effect of a conservative
screening. The recommended choice of δ lies within (5, 6) when pn = 8.2 × 106 and n ∈
[200, 400]. To perform SVD for RIS-PCR, we use fast.svd instead of the usual svd to cope
with the massive number of regressors. Table 3 provides the MSPE, the ECP of 50% PI
and width of the PI, obtained by two different variants of TARP.
Results: The MSPEs of RIS-RP and RIS-PCR are comparable for all the choices on
n. However, RIS-RP yields much better empirical coverage probabilities than RIS-PCR,
30

Table 3: MSPE, ECP and width of PI (in order) obtained by RIS-RP and RIS-PCR for
three values of δ and different training sample sizes (nt ).

nt
200
250
300
350

δ=2
MSPE ECP Width
0.800 0.39 1.059
0.852 0.39 1.102
0.860 0.36 1.126
0.778 0.45 1.210

RIS-RP
δ=5
MSPE ECP Width
0.872 0.42 0.983
0.920 0.42 1.023
0.855 0.44 1.075
0.779 0.48 1.221

δ=8
MSPE
0.855
0.921
0.866
0.829

ECP
0.34
0.35
0.36
0.46

Width
0.928
1.013
1.069
1.219

nt
200
250
300
350

δ=2
MSPE
0.834
0.858
0.845
0.757

RIS-PCR
δ=5
MSPE ECP Width
0.838 0.12 0.192
0.882 0.12 0.289
0.867 0.20 0.511
0.786 0.36 0.886

δ=8
MSPE
0.831
0.896
0.865
0.826

ECP
0.10
0.19
0.20
0.41

Width
0.252
0.420
0.487
0.984

ECP
0.06
0.14
0.14
0.35

Width
0.177
0.355
0.399
0.893

especially when n ≤ 300. The three choices of δ yield comparable results in terms of all
the measures in general. For RIS-RP, δ = 5 results in higher ECP and for RIS-PCR higher
ECP is obtained using δ = 8. Moreover, the choice δ = 8 makes both the procedures
much faster compared to other choices of δ. When the training sample is 350, δ = 2, 5
and 8 select about 290800, 12600 and 7960 variables, respectively, on an average in the
screening stage out of 8.2 × 106 variables. In view of the results in this massive dimensional
dataset, it seems reasonable to use a higher value of δ for filtering out noisy regressors, and
computational convenience.

31

6

Appendix

This section contains proofs of the theorems stated in the paper. We use a generic notation c for
the constants, although all of them may not be equal.
Some Useful Results
Lemma 1. Let εn be a sequence of positive numbers such that nε2n  1. Then under conditions
(a) lnN (εn , Pn ) ≤ nε2n for all sufficiently large n.
2

(b) π (Pnc ) ≤ e−2nεn for all sufficiently large n.

(c) π f : dt (f, f0 ) ≤ ε2n /4 ≥ exp{−nε2n /4} for all sufficiently large n and for some t > 0.
Then

i
h
2
2
Pf0 π {d(f, f0 ) > 4εn |(yn , X)} > 2e−nεn (0.5∧(t/4)) ≤ 2e−nεn (0.5∧(t/4)) .

The proof is given in Jiang (2007).
Lemma 2. Suppose assumption (A1) holds. Let αδ be such that

P

j

x2j |rxj ,y |δ /pn → αδ as n → ∞

and x be a pn × 1 sample vector of the regressors. Then the following holds:
a. The random matrix Rγ described in (2) and (3) satisfies
p

kRγ xk2 /pn →
− cαδ .
b. Let kxγ k2 =

Ppn

2
j=1 xj I(γj

= 1) where γj is the j th element of the vector γ described in (2) and

I(·) is the indicator function, then
p

kxγ k2 /pn →
− cαδ ,

32

where c is the proportionality constant in (2).
The proof is given in the Supplement.
Lemma 3. Let θ ∼ N (0, σθ2 I), then for a given Rγ , x and y the following holds


24 ∆4
(x0 β )2 + ∆2
.
P (|(Rγ x)0 θ − x0 β 0 | < ∆) > exp − 2 0
2
σθ kRγ xk2
σθ kRγ xk2
The proof is given in Guhaniyogi and Dunson (2015).
Proofs of the Theorems. Without loss of generality we consider |xj | < M with M = 1,
j = 1, 2, . . . , pn and d(σ 2 ) = 1, although the proofs go through for any fixed value of M and σ 2 .
Proof of Theorem 1. Define the sequence of events
n
o
2
2
Bn = π {d(f, f0 ) > 4εn |(yn , X)} > 2e−nεn /4 and we need to show P (Bnc ) > 1 − 2e−nε /5 .
We first consider the sequence of events An in assumption (A2), and show that P (Bnc |An ) >
2 /4

1 − 2e−nε

. The proof then follows from assumption (A2) for moderately large n.

The proof of P (Bnc |An ) > 1 − 2e−nε

2 /4

hinges on showing the three conditions of Lemma 1

for the approximating distribution
f (y) = exp{ya(h) + b(h) + c(y)} with h = (Rγ x)0 θ,

(5)

and the true distribution f0 , where Rγ is as given in (4).
Checking condition (a). Let Pn be the set of densities f (y) stated above with parameter |θj | < cn ,
√
j = 1, 2, . . . , mn , where {cn } = {σθ 5nεn } and the model γ is such that γ ∈ Mk , for some
k ∈ {0, 1, . . . , pn }, given An . For any γ the corresponding set of regression parameters can be
n
covered by l∞ balls of the form B = (vj − , vj + )m
j=1 of radius  > 0 and center vj . It takes

33

(cn / + 1)mn balls to cover the parameter space for each of model γ in Pn . There are at most

min{ pkn , pknn } models for each γ under consideration as we are only concerned with models in

An (see assumption (A2)), and there are (pn + 1) possible choices of k. Hence it requires at most
N (, k) ≤ cpnkn +1 (cn / + 1)mn l∞ balls to cover the space of regression parameters Pn , for some
constant c.
Next we find the number of Hellinger balls required to cover Pn . We first consider the KL
distance between f and f0 , then use the fact d(f, f0 ) ≤ (d0 (f, f0 ))1/2 . Given any density in Pn , it
n
can be represented by a set of regression parameters (uj )m
j=1 falling in one of these N (, k) balls
n
B = (vj − , vj + )m
j=1 and pγ = k. More specifically, let fu and fv be two densities in Pn of the

form (5), where u = (Rγ x)0 θ 1 , v = (Rγ x)0 θ 2 with |θi,j | < cn , i = 1, 2 and pγ = k, then
d0 (fu , fv ) =

Z Z

Z Z

fv log



fv
fu



νy (dy)νx (dx)

{y(a(u) − a(v)) + (b(u) − b(v))} fv νy (dy)νx (dx)

 0


Z
b (v)
0
0
=
(u − v) a (uv ) − 0
+ b (uv ) νx (dx).
a (v)

=

The last expression is achieved by integrating with respect to y and using mean value theorem,
where uv is an intermediate point between u and v. Next consider
|u − v| = |(Rγ x)0 θ 1 − (Rγ x)0 θ 2 | ≤ kRγ xkkθ 1 − θ 2 k,

p

using the Cauchy-Schwartz inequality. Now, by Lemma 2 we have kRγ xk2 /pn →
− αδ as
n → ∞ for some constant 0 < αδ < 1. Therefore we can assume that for sufficiently large pn ,
√
√
√
kRγ xk ≤ pn . Also, kθ 1 − θ 2 k ≤ mn . Combining these facts we have |u − v| ≤  mn pn .
√
Similarly max{|u|, |v|} ≤ cn mn pn . These together imply that
(
)
0 (h)|
|b
√
d0 (fu , fv ) ≤  mn pn
sup
|a0 (h)|
sup
.
0
√
√
|h|≤cn mn pn
|h|≤cn mn pn |a (h)|

34

Therefore d(fu , fv ) ≤ εn if we choose
o
n√
mn pn sup|h|≤cn √mn pn |a0 (h)| sup|h|≤cn √mn pn (|b0 (h)|/|a0 (h)|) .
 = ε2n /

Therefore, density fu falls in a Hellinger ball of size εn , centered at fv . As shown earlier, there

are at most N (, k) such balls. Thus, the Hellinger covering number
c
mn
n
N (εn , Pn ) ≤ N (, k) = cpknn +1
+1

(
)
!#mn
"
cn √
|b0 (h)|
kn +1
0
mn pn
sup
+1
= cpn
|a (h)|
sup
0
√
√
ε2n
|h|≤cn mn pn
|h|≤cn mn pn |a (h)|

mn
1
√
kn +1
≤ cpn
D(cn mn pn ) + 1
,
ε2n
where D(R) = R suph≤R |a0 h| suph≤R |b0 (h)/a0 (h)|. The logarithm of the above quantity is no

√
more than log c + (kn + 1) log pn − mn log(ε2n ) + mn log 1 + D(cn mn pn ) , as 0 < ε2n < 1. Using

the assumptions in Theorem 1 condition (a) follows.

n
Checking condition (b) For the Pn defined in condition (a), π(Pnc ) ≤ π(∪m
j=1 |θj | > cn ).
q
Observe that π(|θj | > cn ) ≤ 2 exp{−c2n /(2σθ2 )}/ 2πc2n /σθ2 by Mills ratio. Now for the choice that
p
√
cn = σθ 5nεn the above quantity is 2 exp{−5nε2n /2}/ 10πnε2n . Therefore

π(Pnc ) ≤

mn
X
j=1

π(|θj | > cn ) ≤ 2mn exp{−5nε2n /2}/

p
2
10πnε2n ≤ e−2nεn

for sufficiently large n. Thus condition (b) follows.
Checking condition (c) Condition (c) is verified for t = 1. Observe that

dt=1 (f, f0 ) =

Integrating out y we would get

R

Z Z

f0




f0
− 1 νy (dy)νx (dx).
f

Ey|x [{(f0 /f )(Y ) − 1}] νx (dx). Note that under f and f0 we

35

have same function of y as given in (5) with h = x0 β 0 for f0 .

Therefore, the above can

be written as Ex [{(Rγ x)0 θ − x0 β 0 } g (u∗ )] using mean value theorem where g is a continuous
derivative function, and u∗ is an intermediate point between (Rγ x)0 θ and x0 β 0 . Therefore, if
|(Rγ x)0 θ − x0 β 0 | < ∆n , then |u∗ | < |x0 β 0 | + ∆n . This in turn implies that for sufficiently small
∆n , |g(u∗ )| will be bounded, say by M . Consider a positive constant ∆n . From Lemma 3 we have
P (|(Rγ x)0 θ − x0 β 0 | < ∆n ) =

X
γ

P (|(Rγ x)0 θ − x0 β 0 | < ∆n |γ)π(γ)





24 ∆4
(x0 β 0 )2 + ∆2n
≥ Eγ exp − 2
σθ kRγ xk2
σθ2 kRγ xk2



Zγ
Zγ
24 ∆4n
=
,
Eγ
exp −
(x0 β 0 )2 + ∆2n
pn
pn

(6)



where Zγ = (x0 β 0 )2 + ∆2n / σθ2 kRγ xk2 /pn . By part (a) of Lemma 2, and continuous mapping


p 
theorem Zγ →
− (x0 β 0 )2 + ∆2n / σθ2 cαδ > ∆2n / σθ2 cαδ .

For some non-negative random variable Z and non-random positive numbers p, a and b,

consider the following fact
E









Z
Z
Z
Z
exp −
≥ aP
exp −
>a
p
p
p
p




Z
a
Z
≥ aP
> , exp −
> ab
p
b
p


ap
= aP Z >
, Z < −p log(ab)
b

 ap
< Z < −p log(ab) .
= aP
b

(7)

Replacing Z by Zγ , p by pn and taking a = ∆2n exp{−nε2n /3}/(σθ2 cαδ ), and b = pn exp{−nε2n /3}.



Thus −pn log(ab) = −pn log ∆2n pn exp{−2nε2n /3}/(σθ2 cαδ ) > pn nε2n /2 and apn /b = ∆2n / σθ2 cαδ

36

for sufficiently large n. Therefore the expression in (7) is greater than
∆2n −nε2n /3
e
P
σθ2 cαδ
Note that (x0 β 0 )2 <

Ppn

j=1 |β0,j |



∆2n
1
≤ Zγ ≤ pn nε2n
2
2
σθ cαδ



< K, and the probability involved in the above expression can be

shown to be bigger than some positive constant p for sufficiently large n. Using these facts along
with equation (6), we have P (|(Rγ x)0 θ − x0 β 0 | < ∆n ) > exp{−nε2 /4}. Choosing ∆n < ε2 /(4M )
condition (c) follows.
Proof of Theorem 2. The outline of the proof of Theorem 2 closely follows the arguments given
in the proof of Theorem 1. Therefore we only present those parts of the proof which are different.
2 /4

As in Theorem 1, we show that P (Bnc |An ) > 1 − 2e−nε

by checking the three conditions of

Lemma 1.
The proof of Condition (a) is the same as for Theorem 1, except for the places involving
the projection matrix Rγ . Observe that given a dataset Dn and other tuning parameters we fix
a particular projection matrix Rγ . The only property of Rγ needed to prove condition (a) is
kRγ xk2 ≤ pn for sufficiently large n. To show this we use that fact that Rγ is a matrix with
orthonormal row vectors, and therefore kRγ xk2 ≤ kxγ k2 ≤ pn .
The proof of Condition (b) depends only on the prior assigned on θ, and therefore remains
the same under the conditions of Theorem 2.
The proof of Condition (c) differs from that of Theorem 1 in showing P (|(Rγ x)0 θ − x0 β 0 | <
∆n ) > exp{−nε2 /4} for some constant ∆n . To see this consider a positive constant ∆n . As
before, from Lemma 3 we have

37





(x0 β 0 )2 + ∆2n
24 ∆4
P (|(Rγ x) θ − x β 0 | < ∆n ) ≥ Eγ exp − 2
σθ kRγ xk2
σθ2 kRγ xk2




(x0 β 0 )2 + ∆2n
24 ∆4
≥ Eγ exp − 2
σθ αn kxγ k2
σθ2 kxγ k2



Zγ
Zγ
24 ∆4n
= 0
exp
−
,
E
γ
(x β 0 )2 + ∆2n
pn
αn pn
0

0

(8)



where Zγ = (x0 β 0 )2 + ∆2n / σθ2 kxγ k2 /pn , and αn is as in (A3). From part (b) of Lemma 2,


p 
and continuous mapping theorem we have Zγ →
− (x0 β 0 )2 + ∆2n / σθ2 cαδ > ∆2n / σθ2 cαδ .
For some positive random variable Z and non-random positive numbers p, a and b, consider

the following
E









Z
Z
Z
Z
exp −
≥ aP
exp −
≥a
p
αp
p
αp




Z
a
Z
≥ aP
≥ , exp −
≥ ab
p
b
αp
 ap

= aP
< Z < −αp log(ab) .
b

(9)

(10)

Replacing Z by Zγ , p by pn , α by αn and taking a = ∆2n exp{−nε2n /3}/(σθ2 cαδ ), and b =
pn exp{−nε2n /3}. Thus,



−αn pn log(ab) = −αn pn log ∆2n pn exp{−2nε2n /3}/(σθ2 cαδ ) ∼ 2pn log ∆2n pn /(σθ2 cαδ ) /3 > 2pn /3


for sufficiently large n and apn /b = ∆2n / σθ2 cαδ . Therefore the expression in (10) is greater than


∆2n
2
∆2n −nε2n /3
e
P
≤ Zγ ≤ p n .
3
σθ2 cαδ
σθ2 cαδ
P
n
Note that (x0 β 0 )2 < pj=1
|β0,j | < K, and the probability involved in the above expression

can be shown to be bigger than some positive constant p for sufficiently large n. Using these

facts along with equation (8), we have P (|(Rγ x)0 θ − x0 β 0 | < ∆n ) > exp{−nε2 /4}. Choosing

38

∆n < ε2 /(4M ) condition (c) follows.

References
Achlioptas, D. (2003). Database-friendly random projections: Johnson-Lindenstrauss with binary
coins. Journal of Computer and System Sciences, 66(4):671–687.
Candes, E. and Tao, T. (2007). The Dantzig selector: statistical estimation when p is much larger
than n. The Annals of Statistics, 35(6):2313–2351.
Candès, E. J., Li, X., Ma, Y., and Wright, J. (2011). Robust principal component analysis?
Journal of the ACM, 58(3):Art. 11, 37.
Candes, E. J. and Tao, T. (2005). Decoding by linear programming. IEEE Transactions on
Information Theory, 51(12):4203–4215.
Carvalho, C. M., Polson, N. G., and Scott, J. G. (2009). Handling sparsity via the horseshoe.
In International Conference on Artificial Intelligence and Statistics, volume 5 of Proceedings of
Machine Learning Research, pages 73–80.
Chen, Y. and Dunson, D. B. (2017). Modular Bayes screening for high-dimensional predictors.
ArXiv preprint arXiv:1703.09906.
Consortium, G. et al. (2015). The genotype-tissue expression (GTEx) pilot analysis: Multitissue
gene regulation in humans. Science, 348(6235):648–660.
Dasgupta, S. and Gupta, A. (2003). An elementary proof of a theorem of Johnson and Lindenstrauss. Random Structures & Algorithms, 22(1):60–65.

39

Davenport, M. A., Boufounos, P. T., Wakin, M. B., and Baraniuk, R. G. (2010). Signal processing with compressive measurements. IEEE Journal of Selected Topics in Signal Processing,
4(2):445–460.
Donoho, D. L. (2006).

Compressed sensing.

IEEE Transactions on Information Theory,

52(4):1289–1306.
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., Samworth, R., and Wu, Y. (2009). Ultrahigh dimensional feature selection: Beyond the
linear model. Journal of Machine Learning Research, 10(Sep):2013–2038.
Fard, M. M., Grinberg, Y., Pineau, J., and Precup, D. (2012). Compressed least-squares regression
on sparse spaces. In Association for the Advancement of Artificial Intelligence, pages 1055–1060.
Golub, T. R., Slonim, D. K., Tamayo, P., Huard, C., Gaasenbeek, M., Mesirov, J. P., Coller, H.,
Loh, M. L., Downing, J. R., Caligiuri, M. A., Bloomfield, C. D., and Lander, E. S. (1999).
Molecular classification of cancer: Class discovery and class prediction by gene expression
monitoring. Science, 286(5439):531–537.
Guhaniyogi, R. and Dunson, D. B. (2015). Bayesian compressed regression. Journal of the
American Statistical Association, 110(512):1500–1514.
Hoeting, J. A., Madigan, D., Raftery, A. E., and Volinsky, C. T. (1999). Bayesian model averaging:
A tutorial. Statistical Science. A Review Journal of the Institute of Mathematical Statistics,
14(4):382–417.

40

Ishwaran, H. and Rao, J. S. (2005). Spike and slab variable selection: Frequentist and Bayesian
strategies. The Annals of Statistics, 33(2):730–773.
Jiang, W. (2007). Bayesian variable selection for high dimensional generalized linear models:
Convergence rates of the fitted densities. The Annals of Statistics, 35(4):1487–1511.
Kabán, A. (2014). New bounds on compressive linear least squares regression. In Artificial
Intelligence and Statistics, pages 448–456.
Lappalainen, T., Sammeth, M., Friedländer, M. R., AC‘t Hoen, P., Monlong, J., Rivas, M. A.,
Gonzalez-Porta, M., Kurbatova, N., Griebel, T., Ferreira, P. G., et al. (2013). Transcriptome
and genome sequencing uncovers functional variation in humans. Nature, 501(7468):506.
Li, P., Hastie, T. J., and Church, K. W. (2006). Very sparse random projections. In Proceedings of
the 12th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining,
KDD ’06, pages 287–296.
Maillard, O. and Munos, R. (2009). Compressed least-squares regression. In Advances in Neural
Information Processing Systems, pages 1213–1221.
Pettenuzzo, D., Koop, G., Korobilis, D., et al. (2016). Bayesian compressed vector autoregressions. Working Papers 103, Brandeis University, Department of Economics and International
Businesss School.
Scheetz, T. E., Kim, K.-Y. A., Swiderski, R. E., Philp, A. R., Braun, T. A., Knudtson, K. L.,
Dorrance, A. M., DiBona, G. F., Huang, J., Casavant, T. L., et al. (2006). Regulation of gene
expression in the mammalian eye and its relevance to eye disease. Proceedings of the National
Academy of Sciences, 103(39):14429–14434.

41

Thanei, G.-A., Heinze, C., and Meinshausen, N. (2017). Random projections for large-scale
regression. In Big and Complex Data Analysis, pages 51–68.
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal
Statistical Society. Series B. Methodological, 58(1):267–288.
Tibshirani, R., Bien, J., Friedman, J., Hastie, T., Simon, N., Taylor, J., and Tibshirani, R. J.
(2012). Strong rules for discarding predictors in lasso-type problems. Journal of the Royal
Statistical Society: Series B (Statistical Methodology), 74(2):245–266.
Tierney, L. and Kadane, J. B. (1986). Accurate approximations for posterior moments and
marginal densities. Journal of the American Statistical Association, 81(393):82–86.
Witten, D. M., Tibshirani, R., and Hastie, T. (2009). A penalized matrix decomposition, with
applications to sparse principal components and canonical correlation analysis. Biostatistics,
10(3):515–534.
Zeng, Y. and Breheny, P. (2017). The biglasso package: A memory-and computation-efficient
solver for lasso model fitting with big data in R. ArXiv preprint arXiv:1701.05936.
Zhang, C.-H. et al. (2010). Nearly unbiased variable selection under minimax concave penalty.
The Annals of Statistics, 38(2):894–942.
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.

42

arXiv:1712.02445v1 [] 6 Dec 2017

Supplementary Material

1

Additional Simulation Results

1.1

Results for Larger pn

Here we present the performance of the different methods with respect to mean square prediction error (MSPE) and empirical coverage probability (ECP) of 50% prediction intervals
(PI), along with the width of PIs for a larger choice of pn in each scheme.
Mean Square Prediction Error (MSPE). In Scheme I (see Figure 1) for pn = 3000,
TARP based methods yield the best results, with RIS-PCR having marginally better performance than RIS-RP. These methods are followed by SSR-Ridge and PCR. ISIS based
methods and SSR-EN have overall comparable performance with SSR-EN showing more
stable results, and ISIS-MCP having minimum average MSPE among these 3 methods.
Performance of these six methods are followed by BCR, and that in turn is followed by
RPCR. SSR-LASSO and SPCR have very poor performance in this scenario.
As in the main paper simulation case, TARP based methods have the best overall
performance in Scheme II for pn = 2 ∗ 104 , immediately followed by PCR (see Figure 2).
These three methods are followed by SSR-EN and SSR-Ridge. Unlike for pn = 104 , here
1

Figure 1: Box-plot of MSPE of eleven competing methods for pn = 3000 in Scheme I.

SSR-EN outperforms SSR-Ridge. RPCR is substantially worse but is still better than SSRLASSO, BCR and SPCR (in that order). ISIS based methods had very poor performance.
In Scheme III (see Figure 3), the performance of all the methods are almost similar,
except SSR based methods. Among the SSR based methods, SSR-EN and SSR-LASSO
have very similar results. SSR-EN has lower average MSPE, however SSR-LASSO has a
more stable performance. Unlike for pn = 5000, ISIS based methods often yield large values
of MSPE for pn = 104 .
The relative performance of SPCR and ISIS based methods are unsatisfactory in Scheme
IV (see Figure 4(a)). Among the other methods (see Figure 4(b)), TARP based methods
and PCR have the best performance, followed by BCR and SSR-EN. The relative performance of RPCR and SSR-LASSO improve significantly for pn = 2 ∗ 104 . These two
methods outperform SSR-Ridge, and among them, RPCR yields much better results than
SSR-LASSO.
Empirical Coverage Probability (ECP) and Prediction Interval (PI) Widths.
The relative performance of all the methods is quite similar to that observed in the main
2

Figure 2: Box-plot of MSPE of eleven competing methods for pn = 20000 in Scheme II.

Figure 3: Box-plot of MSPE of eleven competing methods for pn = 10000 in Scheme III.

paper (see Table 1 of the paper and Table 4) for all the simulation schemes. The detailed
description is given below.
The average ECP of SPCR is closest to 0.5 in Scheme I, followed by RPCR, RIS-RP
and SSR-EN, with RIS-RP showing the smallest average width among these methods. The
sd of ECPs and widths of the PIs is highest for PCR, indicating a lack of stability in it’s
performance. RIS-PCR yields some under coverage, however the average width of PI for
RIS-PCR is lowest among all the competing methods. SSR-Ridge has an average ECP of

3

(a) MSPE of all the methods.

(b) MSPE of selected methods.

Figure 4: Box-plot of MSPE of the competing methods for pn = 20000 in Scheme IV.

4

Table 4: Mean and standard deviation(sd) of empirical coverage probabilities of 50% prediction intervals, and mean and sd of width of 50% prediction intervals.
Methods→
Scheme, pn

ISISISISSSR- SSRSACD SACD LASSO EN

SSRRidge

PCR

SPCR RPCR BCR

RISRP

RISPCR

Average and standard deviation (in braces) of empirical coverage probability
I

3∗

103

0.307
(.052)

0.306
(.050)

0.434
(.059)

0.536
(.043)

0.631
(.048)

0.437
(.288)

0.498
(.059)

0.471
(.055)

0.283
(.053)

0.464
(.052)

0.382
(.054)

II

2 ∗ 104

0.334
(.051)

0.333
(.048)

0.422
(.055)

0.577
(.045)

0.637
(.049)

0.365
(.260)

0.505
(.059)

0.454
(.052)

0.280
(.058)

0.488
(.056)

0.312
(.055)

III

104

0.495
(.053)

0.494
(.053)

0.502
(.058)

0.653
(.033)

0.671
(.144)

0.503
(.058)

0.498
(.052)

0.498
(.052)

0.550
(.054)

0.487
(.055)

0.499
(.054)

IV

2 ∗ 104

0.432
(.063)

0.433 0.485
(0.060) (.099)

0.658
(.034)

0.664
(.108)

0.489
(.308)

0.492
(.052)

0.483
(.058)

0.494
(.051)

0.675
(.060)

0.408
(.052)

Average and standard deviation (in braces) of width of the 50% prediction interval
I

3 ∗ 103

4.096
(.522)

3.942
(.531)

6.148
(.410)

7.388
(.543)

8.098
(.484)

8.082 8.441
(7.354) (.608)

7.023
(.492)

5.720
(.411)

5.192
(.152)

4.090
(.335)

6.589 7.649
(7.380) (.929)

6.270
(.565)

5.687
(.465)

4.594
(.204)

2.635
(.209)

1.373
(.070)

1.551
(.073)

1.343
(.079)

1.383
(.077)

2.895 32.430 1.925
(3.808) (7.647) (.234)

1.794
(.144)

2.066
(.142)

1.070
(.071)

II

2 ∗ 104

5.964
(.554)

5.954
(.554)

5.276
(.525)

7.039
(.479)

8.604
(.674)

III

104

1.483
(.174)

1.452
(.144)

1.36
(.072)

1.942
(.122)

5.439 1.362
(2.018) (.072)

IV

2 ∗ 104

18.00 17.98 3.428 1.621
(2.876) (2.789) (2.522) (.167)

4.087
(.634)

1.373
(.070)

about 0.63 and a high width of the PI. ISIS based methods and BCR have lower-coverage
probabilities on average.
For Scheme II, SSR-RP and SPCR have the best overall results, with SPCR having
higher PI widths than SSR-RP. The ECPs of SSR-LASSO, SSR-EN and RPCR are closer
to 0.5 than the remaining methods with SSR-LASSO having a lower average width of the
PIs. The average ECP of PCR is close to 0.5, however there are huge fluctuations in the
widths of PIs (sd= 7.4). Among the other methods, SSR-Ridge shows higher coverage with
a comparatively large width of PI. The remaining four methods show under coverage for
both the choices of pn . Among these four methods, RIS-PCR attains the lowest average
width of PI.

5

All the methods, except SSR-EN and SSR-Ridge have comparable performance in terms
of both the criteria in Scheme III. SSR-EN and SSR-Ridge tend to have higher ECP and
larger width of PI in this setup.
SSR-LASSO, RPCR and BCR yield the best performance in Scheme IV, and BCR
achieves the lowest average width of PI among them. Although the average ECP of SPCR
is close to 0.5, it possesses much higher width than all other methods. SSR-Ridge, SSR-EN
and RIS-RP show higher ECP in general. In terms of PI widths of these three methods,
performance of RIS-RP is the best, followed by SSR-EN. In Scheme IV as well, PCR has
a higher sd of the width of PIs indicating instability in its performance. ISIS-SCAD, ISISMCP and RIS-PCR show some under coverage, and averages of ECPs for all of them are
around 0.4. RIS-RCR yields the lowest overall width of PI among all the methods. Average
width of PI for ISIS based methods and SPCR are much higher than all the other methods
in Scheme IV.

1.2

Effect of sample size on different competing methods

We consider three different choices of n, viz., 150, 200, 250. The results for n = 200 are
presented in the paper and previous sections of Supplementary file. The results corresponding two other choices of n are presented in this section. The relative performance of
the competing methods for different values of n remains similar for the two choices of pn .
Therefore, we present results corresponding to higher values of pn only.
Comparison with respect to MSPE. In Scheme I (see Figures 5(a)-5(b)), SSR-Ridge,
SSR-EN, PCR and TARP based methods work well for lower sample size. Among these
methods, RIS-RP has the best overall performance, immediately followed by RIS-PCR.
6

(a) n = 150

(b) n = 250

Figure 5: Box-plot of MSPE of the competing methods for pn = 3 × 103 in Scheme I.
Among the other methods, BCR and SSR-LASSO yield better performance than others.
The overall scenario changes reasonably for higher sample size, when n = 250. ISIS-MCP
based methods outperform other methods with respect to average MSPE, closely followed
by RIS-RP, ISIS-SCAD and RIS-PCR which also have satisfactory performance. However,
the ISIS based methods have much larger dispersion than TARP based methods, indicating
lack of stability in their results. The other methods do not show much improvement with
increment of sample size. Among these methods, SSR-Ridge yields relatively better results
followed by SSR-EN and PCR.
The relative performance of all these methods does not vary much for different sample
sizes in Scheme II (see Figures 6(a)-6(b)). For both the choices of n, TARP outperforms
all other methods. SSR based methods show reasonable results, although the differences of
their performance with TARP increases with increment of sample size. The other methods
do not perform well for either choices of n. Notably, PCR yields a low median MSPE in
general, however, it results in very high values of MSPE occasionally. For n = 150, the

7

(a) n = 150

(b) n = 250

Figure 6: Box-plot of MSPE of the competing methods for pn = 2 × 104 in Scheme II.
MSPE is as large as 1038, while for n = 250 the MSPE increases up to 3089.
Projection based methods work extremely well in Scheme III for both choices of n (see
Figures 7(a)-7(b)). SSR based methods fail to perform well in this scheme, with SSR-Ridge
showing largest overall MSPE. ISIS based methods show low MSPE in general, however
they occasionally yield higher values of MSPEs.
As the ISIS based method and SPCR yield extremely large values of MSPE in Scheme
IV, we avoid presenting further results on those methods (see Figures 8(a)-8(b)). TARP
yields the best performance for both choices of n in this scheme. SSR-EN also gives comparable performance. Among the other methods, SSR-LASSO and BCR yield reasonable
results. The other three methods tend to have higher values of MSPE in this scheme. As
in Scheme II, PCR tends to have extremely large MSPEs. For n = 250, the highest MSPE
obtained by PCR is around 125, whereas the average MSPE of all other methods is below
7.

8

(a) n = 150

(b) n = 250

Figure 7: Box-plot of MSPE of the competing methods for pn = 104 in Scheme III.
Table 5: Mean and standard deviation (in small font) of empirical coverage probabilities
of 50% prediction intervals.
n→
Schemes →

n = 150

n = 250

I

II

III

VI

I

II

III

VI

Methods↓ pn →

3 × 103

2 × 104

104

2 × 104

3 × 103

2 × 104

104

2 × 104

ISIS-SCAD

0.290.060

0.340.068

0.488.064

0.429.057

0.344.060

0.349.057

0.507.061

0.447.050

ISIS-MCP

0.294.052

0.336.070

0.489.063

0.427.058

0.345.065

0.354.059

0.503.061

0.443.053

SSR-LASSO

0.406.059

0.386.053

0.494.064

0.660.031

0.441.052

0.436.050

0.502.061

0.656.049

ISIS-EN

0.658.035

0.660.035

0.655.031

0.657.030

0.656.031

0.661.035

0.657.033

0.656.030

ISIS-Ridge

0.660.031

0.652.034

0.652.032

0.656.030

0.661.031

0.660.034

0.659.032

0.657.033

PCR

0.407.200

0.450.285

0.494.064

0.484.287

0.421.193

0.416.275

0.502.061

0.504.316

SPCR

0.507.064

0.499.055

0.494.064

0.490.056

0.491.054

0.497.050

0.502.061

0.501.054

RPCR

0.487.052

0.475.056

0.494.065

0.490.053

0.489.057

0.484.051

0.502.061

0.494.058

BCR

0.519.056

0.353.062

0.626.061

0.525.058

0.555.057

0.353.054

0.557.059

0.510.052

RIS-RP

0.508.063

0.460.059

0.492.065

0.723.050

0.598.058

0.504.057

0.499.059

0.704.047

RIS-PCR

0.434.062

0.270.049

0.508.066

0.395.058

0.549.054

0.338.051

0.510.060

0.420.052

Empirical Coverage Probability (ECP) and Prediction Interval (PI) Widths.
Tables 5 and 6 provide the means and standard deviations (sd) of empirical coverage
probabilities (ECP) of 50% prediction intervals (PI) and the widths of the PIs, respectively,
9

(a) n = 150

(b) n = 250

Figure 8: Box-plot of MSPE of the competing methods for pn = 2 × 104 in Scheme IV.
for two different choices of n.
As n increases, ISIS based methods tend to improve rapidly with respect to both the
measures. The other methods show little difference with an increment of n. ISIS based
methods show under coverage for all the simulation schemes except Scheme III. However,
the PI widths are low in general for these methods, except for Scheme IV, where they
result in the highest average width. SSR-EN and SSR-Ridge tend to have higher coverage
probabilities, and higher average PI width. For most of the schemes, these two methods
possess the highest average width. PCR tends to show under coverage for Schemes I and II,
and yield reasonable ECPs for the other two schemes. However, the widths of PIs for PCR
attain the highest variance in all the schemes, except Scheme III, indicating lack of stability
in the results. SPCR and RPCR yield the best performance in terms of ECP. The width of
PI is higher on an average for SPCR than RPCR. For Scheme IV, SPCR yields extremely
large widths of PI. Except for Scheme II, BCR and SSR-LASSO perform reasonably well
in terms of ECP. However, these two methods result in some over-coverage in Scheme III

10

and Scheme IV. With respect to PI width, SSR-LASSO outperforms BCR. TARP based
methods yield the lowest PI width in general. RIS-RP yields an average ECP close to 0.5
in most of the cases, and tend to show high coverage probabilities occasionally. However,
RIS-PCR tends to show low coverage probabilities in general, especially for n = 150.

2

Predictive Calibration in Binary Response Datasets

Apart from measuring the misclassification rates and the area under ROC curve, we validate TARP in terms of it’s ability to quantify uncertainly in real datasets with binary
responses. To this end, we partition the interval [0, 1] into ten equal sub-intervals, viz.,
[0, 0.1), [0.1, 0.2) and so on, and classify the test data points (x, y)i,new to the k th class if
predictive probability of yi,new falls in that class. Next, we consider the squared difference
of the empirical proportion of yi,new = 1 among the data points classified in a given interval
with the middle point of the interval, and consider the mean of these squared differences
(MSD) of all the intervals. If a method is well calibrated, then the MSD would be small.
The following table shows means and standard deviations of MSDs of eleven competing
methods for Golub and GTEx datasets.
Table 7 indicates that the competing methods do not show much differences in MSD
of empirical and predictive probabilities. ISIS based methods, PCR, RPCR and TARP
perform relatively well compared to the other methods. Among these methods, RIS-PCR
and RIS-RP have lowest MSD for Golub and GTEx, respectively. Relative performance of
SSR based methods are deficient in terms of MSD as well. SPCR results in high MSD for
both the datasets too.

11

Table 6: Mean and standard deviation (in small font) of width of 50% prediction intervals.
n→

n = 150

Schemes →

n = 250

I

II

III

VI

I

II

III

VI

Methods↓ pn →

3 × 103

2 × 104

104

2 × 104

3 × 103

2 × 104

104

2 × 104

ISIS-SCAD

4.390.50

6.130.85

7.8744.45

13.762.60

2.990.87

5.870.54

1.490.25

9.674.11

ISIS-MCP

4.410.54

6.120.86

1.460.30

14.202.24

2.780.95

5.900.55

1.400.14

9.714.01

SSR-LASSO

5.800.41

4.950.38

1.360.08

3.920.40

6.450.42

5.490.50

1.360.06

3.760.47

ISIS-EN

9.660.77

8.870.79

3.350.69

2.650.34

9.350.78

8.110.74

3.270.53

2.670.38

ISIS-Ridge

9.310.75

9.641.19

5.392.05

6.320.99

8.120.64

8.830.95

5.731.97

4.490.63

PCR

5.642.89

7.978.33

1.360.08

2.602.82

5.912.97

6.776.41

1.360.06

2.823.53

SPCR

8.440.61

8.220.80

1.360.08

31.089.79

8.310.53

7.930.70

1.360.06

30.2110.41

RPCR

7.760.50

7.730.64

1.360.08

2.820.49

7.720.45

7.540.54

1.360.06

2.750.46

BCR

7.440.48

6.280.52

1.810.07

2.740.23

8.180.42

5.970.38

1.550.05

2.050.14

RIS-RP

6.120.26

4.910.22

1.350.08

2.460.16

6.260.21

4.350.17

1.360.06

1.970.13

RIS-PCR

5.810.33

2.530.20

1.40.08

1.090.08

6.100.37

2.630.20

1.390.06

0.980.06

Table 7: Mean and standard deviation (in bracket) of mean square differences (MSD) of
empirical and predictive probabilities.

3

Methods→
Dataset ↓

ISISSACD

ISISSACD

SSRSSRLASSO EN

SSRRidge

PCR

SPCR

RPCR BCR

RISRP

RISPCR

Golub

2.605
(.036)

2.606
(.035)

2.804
(.000)

2.804
(.000)

2.742
(.028)

2.589
(.012)

3.429
(.073)

2.587
(.017)

2.886
(.032)

2.611
(.045)

2.555
(.044)

GTEx

2.784
(.000)

2.784
(.000)

3.325
(.000)

3.325
(.000)

3.325
(.000)

2.784
(.000)

3.216
(.102)

2.784
(.000)

2.873
(.007)

2.782
(.001)

2.783
(.001)

Mathematical Details

Proof of Lemma 2 a. Consider the conditional expectation and variance of kRγ xk2 given
(γ, x) as follows:
E kRγ xk2 |γ
var kRγ xk2 |γ




= kxγ k2
= 4kxγ k2

(

12

1+



 Ppγ 4 )
1
j=1 xγ,j
−1
,
ψ
2kxγ k2

where xγ includes the regressors j for which γj = 1. The detailed proof is given in Result
1 below.
Next consider the conditional expectation of kRγ xk2 given x is given by


Eγ E kRγ xk2 |γ = Eγ

X

!

x2j I(γj = 1)

j

=c

X

x2j |rxj ,yn |δ ,

j

(1)

where c > 0 is the proportionality constant. Also the conditional variance of kRγ xk2 given
x is given by




varγ E kRγ xk2 |γ + Eγ var kRγ xk2 |γ .

(2)

Considering both the terms in (2) separately as follows:


2

varγ E kRγ xk |γ



= varγ
= c

X
j

X

x2j I(γj

j

!

= 1)


x4j |rxj ,yn |δ 1 − c|rxj ,yn |δ ≤ pn ,

(3)

as given x, γj s are independent, and each |xj | ≤ 1, and qj = c|rj |δ < 1. Again


Eγ var kRγ xk2 |γ



"

= Eγ 4kxγ k2

(

1+



≤ cEγ kxγ k2
X
≤ c
x2j |rxj ,yn |δ ,



 Ppγ 4 )#
1
j=1 xγ,j
−1
ψ
2kxγ k2
(4)

j

for some constant c, as

Ppγ

j=1

x4γ,j < kxγ k2 .

Therefore, from (1), (3) and (4) it can be shown that the expectation of kRγ xk2 /pn
13

converges to the limit cαδ , and variance of the same converges to 0.
Proof of Lemma 2 b. The proof follows from observing that
X

Eγ kxγ k2 = c
x2j |rxj ,yn |δ
j

and


 X 4
varγ E kxγ k2 |γ c
xj |rxj ,yn |δ 1 − c|rxj ,yn |δ ≤ pn .
j

Therefore it can be shown that the expectation of kxγ k2 /pn converges to the limit cαδ , and
variance of the same converges to 0.
Result 1. Consider a random matrix Rγ which depends on another random vector γ
distributed as in (2). Then the conditional distribution of Rγ satisfies the following:
a. E (kRγ xk2 |γ) = kxγ k2 , and

Ppγ 4
b. var (kRγ xk2 |γ) = 4kxγ k2 1 + (ψ −1 − 1) j=1
xγ,j / (2kxγ k2 ) .
Proof of part a. Observe that

2

kRγ xk

=

X

r1,j γj xj ,

j

=

X

r1,j γj xj

j

X

r2,j γj xj , . . . ,

+

r2,j γj xj

!2

j

X
j

rmn ,j γj xj

!0

+ ... +

X

X

j
!2

2

rmn ,j γj xj

j

!2

.

(5)

Now

E

X
j

r1,j γj xj

!2

=E

(
X
j

2
r1,j
γj x2j +

X

r1,j r1,j 0 γj γj 0 xj xj 0

j6=j 0

14

)

=2

X
j

γj x2j = 2kxγ k2 ,

2
as E(ri,j
) = 1 and E(ri,j ri,j 0 ) = 0 as i = 1, 2, . . . , mn , j, j 0 = 1, 2, . . . , pn , and j 6= j 0 .

Proof of part b. From (5) we have

2

var kRγ xk |γ



= var


X X

+

i

ri,j γj xj

j

X

cov

i6=i0


 X


!2 



=

X

var

,

X

i

ri,j γj xj

j

!2

X

ri,j γj xj

j

ri0 ,j γj xj

j

!2 



!2
.

We will consider each term of (6) one by one. Consider the first term. Note that

var

X

ri,j γj xj

j

!2

= var

(
X

2
ri,j
γj x2j +

j

= var

(
X

X

2
ri,j
γj x2j

j

+cov

(
X

)

ri,j ri,j 0 γj γk xj xj 0

j6=k

+ var

(
X

)

ri,j ri,j 0 γj γk xj xj 0

j6=j 0

2
ri,j
γj x2j ,

X

ri,j ri,j 0 γj γj 0 xj xj 0

j6=j 0

j

)

)

.

Consider the first term in (6).

var

(
X
j

2
ri,j
γj x2j

)

=

X
j

=

X


 X
2
2
2
2
cov ri,j
γj x2j , ri,j
var ri,j
γj x2j +
0 γj 0 xj 0


2
γj x4j var ri,j
+

j

=

X
j

j6=j 0

γj x4j



E

4
ri,j

15



X

2
2
γj x2j γj 0 x2j 0 cov ri,j
, ri,j
0

j6=j 0

−E

2

2
ri,j



=2





X
1
−1
γj x4j .
ψ
j

(6)

Again,

var
=

X
j6=j 0

(
X

ri,j ri,j 0 γj γk xj xj 0

j6=j 0

)


2 2
γj γk x2j x2j 0 E ri,j
ri,j 0 +
=4

X

=E
X

X

ri,j ri,j 0 γj γk xj xj 0

j6=j 0

!2

,

γj γk γj 0 γk0 x2j x2j 0 x2k x2k0 E (ri,j ri,j 0 ri,k ri,k0 )

(j, j 0 ) 6= (k, k0 )
j 6= j 0 , k 6= k0

γj γk x2j x2j 0 ,

j6=j 0

as the other term will be zero. Next
o
P
2
2
0 γj γj 0 xj xj 0
,
r
r
γ
x
r
0
i,j
i,j
j
j
j6=j
j i,j

P P
2
= j k6=k0 γj x2j , γk γk0 xk xk0 cov ri,j
, ri,k ri,k0 = 0.

cov

nP

Therefore the first term in (6) is
X

var

i

X

ri,j γj xj

j

!2

#
"
X
X
1
κ
4
2 2
=
γj γk xj xj 0
(n − 1)
γj xj +
mn
j
j6=j 0

!2 

X
X
1
= 2
−1
γj x4j + 2
γj x2j  .
ψ
j
j

(7)

Now consider the last term in (6).

cov


 X


= cov

j

(
X
j

ri,j γj xj

!2

2
ri,j
γj x2j +

,

X
j

X

ri0 ,j γj xj

!2 



ri,j ri,j 0 γj γj 0 x2j x2j 0 ,

j6=j 0

X
k

16

ri20 ,k γk x2k +

X

k6=k0

ri0 ,k ri0 ,k0 γk γk0 x2k x2k0

)

= cov

(
X

2
γj x2j ,
ri,j

j

+cov

(
X

X

ri20 ,k γk x2k

k

ri,j ri,j 0 γj γj 0 x2j x2j 0 ,

j6=j 0

+cov

(
X

X

)

+ cov

ri20 ,k γk x2k

k

ri,j ri,j 0 γj γj 0 x2j x2j 0 ,

j6=j 0

X

(
X

)

2
γj x2j ,
ri,j

k6=k0

ri0 ,k ri0 ,k0 γk γk0 x2k x2k0

k6=k0

j

ri0 ,k ri0 ,k0 γk γk0 x2k x2k0

X

)

.

)

(8)

Consider the first term of (8).

cov

(
X

2
ri,j
γj x2j ,

j

X
k

ri20 ,k γk x2k

)

=

XX
j

k


2
γj x2j γk x2k cov ri,j
, ri20 ,k = 0.

Similarly, it can be shown that all other terms in (8) are zero. Combining the above result
and (7) the proof follows.

17

