Scalable Bayesian Variable Selection Using Nonlocal

arXiv:1507.07106v4 [stat.ME] 18 Jan 2017

Prior Densities in Ultrahigh-dimensional Settings
Minsuk Shin∗, Anirban Bhattacharya†, and Valen E. Johnson‡
Department of Statistics, Texas A&M University, Texas, U.S.A

Abstract
Bayesian model selection procedures based on nonlocal alternative prior densities are extended to ultrahigh dimensional settings and compared to other variable selection procedures
using precision-recall curves. Variable selection procedures included in these comparisons include methods based on g-priors, reciprocal lasso, adaptive lasso, scad, and minimax concave
penalty criteria. The use of precision-recall curves eliminates the sensitivity of our conclusions
to the choice of tuning parameters. We find that Bayesian variable selection procedures based
on nonlocal priors are competitive to all other procedures in a range of simulation scenarios,
and we subsequently explain this favorable performance through a theoretical examination of
their consistency properties. When certain regularity conditions apply, we demonstrate that
the nonlocal procedures are consistent for linear models even when the number of covariates
p increases sub-exponentially with the sample size n. A model selection procedure based on
Zellner’s g-prior is also found to be competitive with penalized likelihood methods in identifying the true model, but the posterior distribution on the model space induced by this method
is much more dispersed than the posterior distribution induced on the model space by the
∗

Electronic address: minsuk@stat.tamu.edu
Electronic address: anirbanb@stat.tamu.edu
‡
Electronic address: vjohnson@stat.tamu.edu
†

1

nonlocal prior methods. We investigate the asymptotic form of the marginal likelihood based
on the nonlocal priors and show that it attains a unique term that cannot be derived from the
other Bayesian model selection procedures. We also propose a scalable and efficient algorithm
called Simplified Shotgun Stochastic Search with Screening (S5) to explore the enormous
model space, and we show that S5 dramatically reduces the computing time without losing the
capacity to search the interesting region in the model space, at least in the simulation settings
considered. The S5 algorithm is available in an R package BayesS5 on CRAN.

Key words: Bayesian variable selection; Nonlocal prior; Precision-recall curve; Strong model
consistency; Ultrahigh-dimensional data.

1

Introduction

In the context of hypothesis testing, Johnson and Rossell (2010) defined nonlocal (alternative) priors as densities that are exactly zero whenever a model parameter equals its null value. Nonlocal
priors were extended to model selection problems in Johnson and Rossell (2012), where product
moment (pMoM) prior and product inverse moment (piMoM) prior densities were introduced as
priors on a vector of regression coefficients. In p ≤ n settings, model selection procedures based
on these priors were demonstrated to have a strong model selection property: the posterior probability of the true model converges to 1 as the sample size n increases. More recently, Rossell et al.
(2013) and Rossell and Telesca (2017) proposed product exponential moment (peMoM) prior densities that have similar behavior to piMoM densities near the origin. However, the behavior of
nonlocal priors in p  n settings remains understudied to date (particularly in comparison to other
commonly used variables selection procedures), which serves as the motivation for this article.
We undertook a detailed simulation study to compare the performance of nonlocal priors in p 
n settings under sparsity with a host of penalization methods including the least absolute shrinkage
and selection operator (lasso; Tibshirani (1996)), smoothly clipped absolute deviation (scad; Fan
and Li (2001)), adaptive lasso (Zou, 2006), minimum convex penalty (mcp; Zhang (2010)), and the

2

reciprocal lasso (rlasso), recently been proposed by Song and Liang (2015). The penalty function
of the rlasso is equivalent to the negative log-kernel of nonlocal prior densities; further connections
are described in Section 5. As a natural Bayesian competitor, we also considered the widely used gprior (Zellner, 1986; Liang et al., 2008), which is a local prior in the sense of Johnson and Rossell
(2010). We used precision-recall curves (Davis and Goadrich, 2006) as a basis for comparison
between methods. These curves eliminate the effect of the choice of tuning parameters for each
method so that the comparison across different methods can be transparent. It has been argued
(Davis and Goadrich, 2006) that in cases where only a tiny proportion of variables are significant,
precision-recall curves are more appropriate tools for comparison than are the more widely used
receiver operating characteristic curves. While the ROC curves present a trade-off between the
type I error and the power of a decision procedure, precision-recall curves examine the trade-off
between the power and the false discovery rate.
Our studies indicate that Bayesian procedures based on nonlocal priors and the g-prior perform
better than penalized likelihood approaches in a sense that they achieve a lower false discovery
rate, while maintaining the same power of the decision procedure. Posterior distributions on the
model space based on nonlocal priors were found to be more tightly concentrated around the
maximum a posteriori model than the posterior based on g-priors, implying that they had a faster
rate of posterior concentration. We also identified the oracle hyperparameter that maximizes the
posterior probability of the true model for the Bayesian procedures. The growth-rate of these oracle
hyperparameters with p also offers an interesting contrast between nonlocal and local priors. In
the case of g-priors, the oracle value of g varied between 7.83 × 108 and 4.29 × 1013 as p ranged
between 1000 and 20000. For the same range of p, the oracle value of τ varied between 1.97
and 3.60, where τ is the tuning parameter for nonlocal priors described in Section 2. George and
Foster (2000) argued from a minimax perspective that the g parameter should satisfy g  p2 ,
which explains the large values of the optimal g. However, using asymptotic arguments to obtain
default hyperparameters is difficult because the constant of proportionality is typically unknown.
Moreover, when g is very large, the g-prior assigns negligible prior mass at the origin, essentially

3

resulting in a nonlocal like prior. A similar point can be made about the recently proposed Bayesian
shrinking and diffusing (BASAD) priors (Narisetty and He, 2014). On the other hand, the optimal
hyperparameter value for the nonlocal priors is stable with increasing p, growing at a very slow
rate.
Motivated by this empirical finding, we studied properties of two classes of nonlocal priors
allowing the hyperparameter τ to scale with p. Using a fixed value of τ , it seems that strong model
selection consistency is possible only when p ≤ n (Johnson and Rossell, 2012). In this article, we
establish that nonlocal priors can achieve strong model selection consistency even when the number
of variables p increases sub-exponentially in the sample size n, provided that the hyperparameter τ
is asymptotically larger than log p. This theoretical result is consistent with our empirical finding.

2

Nonlocal prior densities for regression coefficients

We consider the standard setup of a Gaussian linear regression model with a univariate response
and p candidate predictors. Let y = (y1 , . . . , yn )T denote a vector of responses for n individuals
and X an n × p matrix of covariates. We denote a model by k = {k1 , . . . , k|k| }, with 1 ≤ k1 <
. . . < k|k| ≤ p. Given a model k, let Xk denote the design matrix formed from the columns of
Xn corresponding to model k and βk = (βk,1 , . . . , βk,|k| )T the regression coefficient for model k.
Under each model k, the linear regression model for the data is

y = Xk βk + ,

(1)

where  ∼ Nn (0, σ 2 In ). Let t denote the true, or data-generating model and let βt0 be the true
regression coefficient under model t. We assume that the true model is fixed but unknown.
Given a model k, the product exponential moment (peMoM) prior density (Rossell et al., 2013;

4

0.10
0.00

0.05

density

0.15

peMoM
piMoM

−10

−5

0

5

10

parameter

Figure 1: Nonlocal prior density functions for a single regression coefficient with τ = 5; for the
piMoM prior, r = 1.
Rossell and Telesca, 2017) for the vector of regression coefficients βk is defined as

2

π(βk | σ , τ, k) = C

−|k|

|k|
Y

2
2
exp{−βk,j
/(2σ 2 τ ) − τ /βk,j
}.

(2)

j=1

The normalizing constant C can be explicitly calculated as
Z

∞

exp{−t2 /(2σ 2 τ ) − τ /t2 }dt = (2πσ 2 τ )1/2 exp{−(2/σ 2 )1/2 },

C=

(3)

−∞

since

R

exp{−µ/t2 − ζt2 }dt = (π/ζ)1/2 exp{−2(µζ)1/2 }.

Second, for a fixed positive integer r, the product inverse-moment (piMoM) prior density
(Johnson and Rossell, 2012) for βk is given by

2

π(βk | σ , τ, k) = C

∗−|k|

|k|
Y
2
[(βk,j )−2r exp{−τ /βk,j
}],

(4)

j=1

where C ∗ = τ −r+1/2 Γ(r − 1/2) for r > 1/2, where Γ(·) is the gamma function.
The piMoM and peMoM prior densities are nonlocal in the sense that the density value at the
origin is exactly zero. This feature of the densities for a single regression coefficient is illustrated
in Figure 1. Since the piMoM prior densities and the peMoM prior densities have the same term

5

exp{−τ /β 2 } that controls the behavior of the density function around the origin, they attain almost
the same shape of the density function at the origin, which yields the same theoretical results in an
asymptotic sense. Further details regarding this point are discussed in Section 4.
We focus on these two classes of nonlocal priors in the sequel. Note that in both (2) and (4),
π(βk ) = 0 when βk = 0; a defining feature of nonlocal priors. The distinction between the peMoM
and the piMoM priors mainly involves their tail behavior. Whereas peMoM priors possess Gaussian tails, the piMoM prior densities have inverse polynomial tails. For example, piMoM densities
with r = 1 have Cauchy-like tails, which has implications for their finite sample consistency and
asymptotic bias in posterior mean estimates of regression coefficients. Since similar conditions are
later imposed on the hyperparameter τ appearing in (2) and (4), at the risk of some ambiguity we
use the same notation for the two hyperparameters in these equations.
In addition to imposing priors on the regression parameters given a model, we need to place a
prior on the space of models to complete the prior specification. We consider a uniform prior on
the model space restricted to models having size less than or equal to qn , with qn < n, i.e.,

π(k) ∝ I(|k| ≤ qn ),

(5)

where I(·) denotes the indicator function and with a slight abuse of notation, we denote the prior
on the space of models by π as well. Similar priors have been considered in the literature by
Jiang (2007) and Liang et al. (2013). Since the peMoM and piMoM priors already induce a strong
penalty on the size of the model space (see Section 4), we do not need to additionally penalize
larger models using, for example, model space priors of the type discussed in Scott and Berger
(2010).
Under a peMoM prior (2) on the regression coefficients, the marginal likelihood mk (y) under
model k given σ 2 can be obtained by integrating out βk , resulting in
n

ek /(2σ 2 )},
mk (y) = (2πσ 2 )− 2 C −|k| Qk exp{−R

6

where

ek = y T (In − P
e k )y, P
e k = Xk (XkT Xk + 1/τ Ik )−1 XkT ,
R
Z
|k|
X
T e −1
2
2
e
e
Qk =
exp{−(βk − βk ) Σk (βk − βk )/(2σ ) −
}dβk ,
τ /βk,j

(6)

j=1

βek = (XkT Xk + 1/τ Ik )−1 XkT y,

e k = (X T Xk + 1/τ Ik )−1 .
Σ
k

Similarly, the marginal likelihood using the piMoM prior densities (4) can be expressed as
n

mk (y) = (2πσ 2 )− 2 C ∗−|k| Q∗k exp{−Rk∗ /(2σ 2 )}, where
Rk∗ = y T (In − Pk )y, Pk = Xk (XkT Xk )−1 XkT ,
Z Y
|k|
|k|
X
−2r
∗−1
∗
T
2
2
Qk =
βk,j exp{−(βk − βbk ) Σk (βk − βbk )/(2σ ) −
τ /βk,j
}dβk ,
j=1
T

(7)

j=1
−1

T

βbk = (Xk Xk ) Xk y,

Σ∗k

−1

T

= (Xk Xk ) .

The integrals for Qk and Q∗k cannot be obtained in closed forms, so for computational purposes
we make Laplace approximations to mk (y). The expressions for the marginal likelihood derived
here is nevertheless important for our theoretical study in Section 4.

3
3.1

Numerical results
Simulation studies using precision-recall curves

To illustrate the performance of nonlocal priors in ultrahigh-dimensional settings and to compare
their performance with other methods, we calculated precision-recall curves (Davis and Goadrich,
2006) for all selection procedures. A precision-recall curve plots the precision = TP/(TP + FP),
versus recall (or sensitivity) = TP/(TP + FN), where TP, FP and FN respectively denote the number
of true positives, false positives, and false negatives, as the tuning parameter is varied. The efficacy
of a procedure can be measured by the area under the precision-recall curve; the greater the area,

7

the more accurate the method. Since both precision and recall take values in [0, 1], the area under
the curve for an ideal precision-recall curve is 1. We used two (n, p) combinations, namely (n, p) =
(400, 10000) and (n, p) = (400, 20000), and plotted the average of the precision-recall curves
obtained from 100 independent replicates of each procedure.
We compared the performance of peMoM and piMoM priors to a number of frequentist penalized likelihood methods: lasso (Tibshirani, 1996), adaptive lasso (Zou, 2006), scad (Fan and Li,
2001), and minimax concave penalty (Zhang, 2010). We used the R package ncvreg to fit these penalized likelihood methods. We also included reciprocal lasso in our simulation studies. However,
due to computational constraints involved in implementing the full rlasso procedure, we followed
the recommendation in Song and Liang (2015) and instead implemented the reduced rlasso. The
reduced rlasso procedure is a simplified version of rlasso that uses the least square estimators of β
when minimizing the rlasso objective function.
We considered Zellner’s g-prior (Zellner, 1986; Liang et al., 2008) as a competing Bayesian
method, with βk | k, σ 2 ∼ N (0, gσ 2 (XkT Xk )−1 ) and g is the tuning parameter. With the prior
π(σ 2 ) ∝ 1/σ 2 , the marginal likelihood mk (y) ∝ (1 + g)−|k|/2 {1 + g(1 − Dk2 )}−(n−1)/2 can be
obtained in a closed form; see for example, Liang et al. (2008, pp 412), where Dk2 is the ordinary
coefficient of determination for the model k.
A uniform model prior (5) was considered for all Bayesian procedures. This prior was chosen
for several reasons. First, construction of the PR curves requires maximization over model hyperparameters, which is most easily achieved if there is only one unknown hyperparameter. We also
wished to avoid providing an advantage to the Bayesian methods by introducing additional tuning
parameters into these methods that were not present in the penalized likelihood methods. Furthermore, the use of non-uniform priors on the model space introduces (at least) one more degree of
freedom into the comparisons between methods, and our intent was to compare the effects of the
penalties imposed on regression coefficients by both penalized likelihood and Bayesian methods.
At first blush, this might appear to put Bayesian methods like those based on the g-prior at a disadvantage, since such methods do not yield consistent variable selection even in p < n settings

8

without prior sparsity penalties on the model space (when g is held fixed as n increases). However,
in the construction of our PR curves, we allowed prior hyperparameters to increase with n, which
effectively allowed the Bayesian methods to impose additional sparseness penalties through the
introduction of large hyperparameter values.
We arbitrarily fixed r = 1 for the piMoM prior (4) and used an inverse-gamma prior on σ 2
with parameters (0.1, 0.1) for the peMoM, piMoM priors, and g-priors. Posterior computations
for the peMoM, piMoM and g-priors were implemented using the Simplified Shotgun Stochastic
Search with Screening (S5) algorithm described in Section 7. The maximum a posteriori model
was used in each case to summarize the model selection performance. The precision-recall curves
are drawn by varying the hyperparameters (τ for the nonlocal priors and g for the g-priors), so the
comparison between the model selection based on the nonlocal priors and the g-prior is free of
the choice of hyperparameters. Because of their high computational burden, we could not include
BASAD (Narisetty and He, 2014) in the comparisons.
For each simulation setting, we simulated data according to a Gaussian linear model as in (1)
with the fixed true model t = {1, 2, 3, 4, 5} with the true regression coefficient βt0 = (0.50, 0.75, 1.00, 1.25, 1.50)T
and σ = 1.5. Also, the signs of the true regression coefficients were randomly determined with
probability one-half. Each row of X was independently generated from a N (0, Σ) distribution with
one of the following covariance structures:
Case (1): compound symmetry design; Σjj 0 = 0.5, if j 6= j 0 and Σjj = 1, 1 ≤ j, j 0 ≤ p.
0

Case (2): autoregressive correlated design; Σjj 0 = 0.5|j−j | , 1 ≤ j, j 0 ≤ p.
Case (3): isotropic design; Σ = Ip .
Figure 2 plots the precision-recall curves averaged over 100 simulation replicates for the different methods across the two (n,p) pairs and the three covariate designs. From Figure 2, it is evident
that the precision-recall curves for the peMoM and piMoM priors have an overall better performance than the penalized likelihood methods lasso, adaptive lasso, scad, and mcp. For decision
procedures having the same power, this implies that the nonlocal priors achieve lower false discovery rates. As discussed in Section 5, since the reduced rlasso shares the same nonlocal kernel
9

as the nonlocal priors, it has a similar selection performance. The figure also shows that Zellner’s
g-prior attains comparable performance with the nonlocal priors in terms of the precision-recall
curves.

3.2

Further comparison with Zellner’s g-prior

The similarity of the performances of the g-prior and the nonlocal priors in terms of precisionrecall curves begs for closer comparisons of these procedures. For this reason, we also investigated
the concentration of the posterior densities around their maximum models. To this end, we fixed
p = 20, 000 and varied n from 150 to 400; the data generating mechanism was exactly the same
as in Section 3.1. The left column of Figure 3 displays the posterior probability of the true model
under the peMoM, piMoM and g-prior models versus n for the three covariate designs in Section
3.1. The plot shows that the posterior probability of the true model increases with n for all three
methods, with the peMoM and piMoM priors almost uniformly dominating the g-prior, implying
a higher concentration of the posterior around the true model for the nonlocal priors.
This tendency is confirmed in the right panel of Figure 3, where we plot the number of models
b | y) > 0.001, where k
b is the maximum
k which achieve a posterior odds ratio π(k | y)/π(k
a posteriori model. This plot clearly shows that the posterior distribution on the model space
from the g-priors is more diffuse than those obtained using the nonlocal prior methods. These
comparisons were based on fitting the hyperparameters g and τ at their oracle value, i.e., the value
which maximized the posterior probability of the true model for a given value of n.
The magnitudes of the oracle hyperparameters under each model also present an interesting
contrast between the local and nonlocal priors. We observed that the oracle value of g increased
rapidly with p, whereas the oracle value of τ was much more stable. This phenomenon is illustrated
in Table 1 that shows the oracle hyperparameter value averaged over 100 replicates for the three
different covariate designs in Section 3.1. For this comparison, we fixed n = 400 and varied p
between 1000 and 20, 000; five representative values are displayed. The oracle values for g are
on a completely different scale from the oracle values τ , and they vary more with p. This table
10

1.0

1.0

0.8

0.8

peMoM
piMoM

Precision
0.4
0.6

Precision
0.4
0.6

g−prior
rrLASSO
AdapLASSO

0.2

0.2

LASSO
SCAD

0.0

0.0

MCP
0.4
0.6
Recall

0.8

1.0

0.0

0.2

0.4
0.6
Recall

0.8

1.0

0.0

0.2

0.4
0.6
Recall

0.8

1.0

0.0

0.2

0.4
0.6
Recall

0.8

1.0

0.0

0.2

0.4
0.6
Recall

0.8

1.0

0.0

0.2

0.4
0.6
Recall

0.8

1.0

0.8
Precision
0.4
0.6
0.2
0.0
1.0
0.8
Precision
0.4
0.6
0.2
0.0

0.0

0.2

Precision
0.4
0.6

0.8

1.0

0.0

0.2

Precision
0.4
0.6

0.8

1.0

0.2

1.0

0.0

Figure 2: Plot of the mean precision-precision curves over 100 datasets with (n, p) =
(400, 10000)(first column) and (n, p) = (400, 20000)(second column). Top: case (1); middle:
case (2); bottom: case (3).
11

80 100
60
40

# of models

20
0

posterior probability

0.0 0.2 0.4 0.6 0.8 1.0

peMoM
piMoM
g−prior

150

200

250

300

400

150

200

300

400

300

400

300

400

80 100
60
40

# of models

20
0
150

200

250

300

400

150

200

250

n

60
40
0

20

# of models

0.0 0.2 0.4 0.6 0.8 1.0

80 100

n

posterior probability

250

n

0.0 0.2 0.4 0.6 0.8 1.0

posterior probability

n

150

200

250

300

400

150

n

200

250

n

Figure 3: Averaged posterior true model probability and the number of models which attain the
posterior odds ratio, with respect to the maximum a posteriori model, larger than 0.001 with the
fixed p = 20000 and varying n. Top: case (1); middle: case (2); bottom: case (3).

12

Table 1: Optimal hyperparameters for Bayesian model selection methods

Case (1)

Case (2)

Case (3)

p = 1000
peMoM
2.24
piMoM
2.16
g-prior 7.83 × 108
peMoM
1.97
piMoM
1.97
g-prior 8.56 × 109
peMoM
2.66
piMoM
2.61
g-prior 1.26 × 1012

The number of predictors
p = 2000
p = 5000
p = 10000
2.72
2.88
3.32
2.59
2.70
3.04
2.87 × 109 3.05 × 109 9.66 × 109
2.29
2.34
2.75
2.20
2.32
2.66
2.55 × 1010 2.62 × 1010 6.58 × 1010
3.00
3.00
3.10
2.94
2.94
2.94
8.84 × 1012 9.67 × 1012 6.81 × 1012

p = 20000
3.60
3.26
1.70 × 1010
3.00
2.86
1.25 × 1011
3.60
3.46
4.29 × 1013

confirms the recommendations in George and Foster (2000) for setting g = p2 based on minimax
arguments. However, the finite sample behavior of the optimal choice of g is unclear, which means
that the large variance of the optimal hyperparameter value is likely to hinder the selection of g in
real applications. Finally, we note that such large values of g effectively convert the local g-priors
into nonlocal priors by effectively collapsing the g-prior density to 0 at the origin.

4

Model selection consistency

The empirical performance of the peMoM and piMoM priors suggests that the hyperparameter τ
should be increased slowly with p. While Johnson and Rossell (2012) were able to show strong
selection consistency with a fixed value of τ , it is not clear whether their proof can be extended to
p  n cases. Motivated by the empirical findings of the last section, we next investigated the strong
consistency properties of peMoM and piMoM priors when τ was allowed to grow at a logarithmic
rate in p. We found that in such cases, both peMoM and piMoM priors achieve strong model
selection consistency under standard regularity assumptions when p increases sub-exponentially
with n, i.e., log p = O(nα ) for α ∈ (0, 1).
Henceforth, we use τn,p instead of τ to denote the hyperparameter in the peMoM and piMoM
priors in (2) and (4) respectively. The normalizing constants for these priors is now denoted by Cn,p

13

∗
and Cn,p
respectively. Before providing our theoretical results, we first state a number of regularity

conditions. Let νj (A) denote the j-th largest nonzero eigenvalue of an arbitrary matrix A, and let

νk∗ =

min
1≤j≤min(n,|k|)

νj (XkT Xk /n),

νk∗ =

max
1≤j≤min(n,|k|)

νj (XkT Xk /n).

(8)

For sequences an and bn , an  bn indicates bn = O(an ), and an  bn indicates bn = o(an ).
With this notation, we assume that the following regularity conditions apply.
Assumption 1. There exists α ∈ (0, 1) such that log p = O(nα ).
Assumption 2. log p ≺ τn,p ≺ n.
Assumption 3. |k| ≤ qn , where qn ≺
Assumption 4. min νk∗ 
k:|k|≤qn

τn,p
.
log p

τn,p
.
n

Assumption 5. C1 < νt∗ ≤ νt∗ < C2 for some positive constants C1 and C2 .

Several comments regarding these conditions are worth making. Assumption 1 allows p to
grow sub-exponentially with n. Our theoretical results continue to hold when p grows polynomially in n, i.e., at the rate O(nγ ) for some γ > 1. Assumption 2 reflects our empirical findings about
the oracle τ ≡ τn,p in Section 3.1, which was observed to grow slowly with p. We need the bound
on qn in Assumption 3 to ensure that the least square estimator of a model is consistent when a
model contains the true model. In the p ≤ n setting, Johnson and Rossell (2012) assumed that
all eigenvalues of the Gram matrix (XkT Xk )/n are bounded above and below by global constants
for all k. However, this assumption is no longer viable when p  n and we replace that by Assumption 4, where the minimum of the minimum eigenvalue of (XkT Xk )/n over all submodels k
with |k| ≤ qn is allowed to decrease with increasing n and p. Assumption 4 is called the sparse
Riesz condition and is also used in Chen and Chen (2008) and Kim et al. (2012). Narisetty and
He (2014) showed that Assumption 4 holds with overwhelmingly large probability when the rows

14

of the design matrix are independent with an isotropic sub-Gaussian distribution. Even though the
assumption of sub-Gaussian tails on the covariates is difficult to verify, the results in Narisetty and
He (2014) show that Assumption 4 can be satisfied for some sequence of design matrices.
We now state a Theorem that demonstrates that model selection procedures based on the peMoM and piMoM nonlocal prior densities achieve strong consistency under the proposed regularity conditions. A proof of the Theorem is provided in the Supplemental Materials.
Theorem 1. Suppose σ 2 is known and that Assumptions 1 – 5 hold. Let π(t | y) denote the
posterior probability of the true model obtained under a peMoM prior (2). Also, assume a uniform
prior on all models of size less than or equal to qn , i.e., π(k) ∝ I(|k| ≤ qn ). Then, π(t | y)
converges to one in probability as n goes to ∞.
Corollary 2. Assume the conditions of the preceding Theorem apply. Let π(t | y) denote the
posterior probability of the true model obtained under a piMoM prior density (4). Then, π(t | y)
converges to one in probability as n goes to ∞.
We note that these results apply also if a beta-Bernoulli prior is imposed on the model space
as in Scott and Berger (2010), because the effect of that prior is asymptotically negligible when
|k| ≤ qn ≺ n.
In most applications, σ 2 is unknown, and it is thus necessary to specify a prior density on it. By
imposing a proper inverse gamma prior density on σ 2 , we can obtain the strong model consistency
result stated in the Theorem below. The proof is again deferred to the Supplemental Materials.
Theorem 3. Suppose σ 2 is unknown and a proper inverse gamma density with parameters (a0 , b0 )
is assumed for σ 2 . Also, let π(t | y) denote the posterior probability of the true model evaluated using peMoM priors. Then if Assumptions 1 – 5 are satisfied, π(t | y) converges to one in
probability as n goes to ∞.
Corollary 4. Suppose the conditions of the preceding Theorem apply, but that π(t | y) now denotes
the posterior probability of the true model obtained under a piMoM prior density. Then π(t | y)
converges to one in probability as n goes to ∞.
15

5

Connections between nonlocal priors and reciprocal lasso

In this section, we highlight the connection between the rlasso of Song and Liang (2015) and
Bayesian variable selection procedures based on our nonlocal priors. We begin by noting that the
objective function g(βk ; k) of rlasso on a model k can be expressed as follows:

g(βk ; k) = ky −

Xk βk k22

+

|k|
X

τn,p /|βk,j |.

(9)

j=1

The optimal model is selected by minimizing this objective function with respect to βk and k.
P
It is clear that the penalty function |k|
j=1 τn,p /|βk,j | in (9) is similar to the negative log-density
of piMoM nonlocal priors as proposed in Johnson and Rossell (2012, pp 659) and Johnson and
Rossell (2010, pp 149). The main difference between the nonlocal prior version of rlasso and the
piMoM-type prior densities proposed in the previous section is the power of β in the exponential
kernels. For the rlasso penalty this power is 1, while for piMoM-type prior densities it is 2. The
implications of this difference are apparent from the following proposition.
Proposition 5. For a given model k, suppose that βek∗ is the minimizer of the objective function (9),
and again let βbk denote the least square estimator of β under model k. Assume that τn,p ≺ n, and
there exist strictly positive contants CL and CU such that CL < νk∗ ≤ νk∗ < CU . Then, for any
∗n  (τn,p /n)1/3 ,
h
i
∗
∗
e
b
P βk ∈
/ R βk ; n → 0,
where R(u; ) = {x ∈ R|k| : |xj − uj | ≤ , j = 1, . . . , |k|}.
The proposition shows that under standard conditions on the eigenvalues of the Gram matrix
XkT Xk /n, the estimator derived from (9) is asymptotically within (τn,p /n)1/3 distance of the least
squares estimator βbk . On the other hand, results cited in the previous section show that maximum
a posteriori estimators obtained from the piMoM-type prior densities reside at an asymptotic distance of (τn,p /n)1/4 from the least squares estimator. Variable selection procedures based on both
forms of piMoM priors thus achieve adaptive penalties on the regression coefficients in the sense
16

described in Song and Liang (2015).
Although rlasso is proposed as a penalized likelihood approach, the computational procedure
to optimize its objective function is quite different from the other penalized likelihood methods.
The resulting computational complexity of this optimization procedure, which contains a discontinuous penalty function, is NP-hard. This suggests that the formulation of this nonlocal penalty
in a penalized likelihood framework is unlikely to provide significant computational advantages
over related Bayesian model selection procedures, even though the inferential advantages of the
Bayesian framework are lost.

6

Asymptotic behavior of marginal likelihoods based on nonlocal priors

From Lemma 1 in the Supplemental Materials, it follows that the asymptotic log-marginal likelihood of a model k based on a peMoM or piMoM prior density can be expressed as

log π(k | y) = l(βbk ) + log Qk − |k| log Cn,p
≈ l(βbk ) −

|k|
X


pτn,p βbk,j + C,

j=1

for some constant C, βbk is the maximum likelihood estimator under the model k, i.e. βbk =
(XkT Xk )−1 XkT y, and



pτn,p βbk,j ≈




(nτ

n,p uk )

1/2

, if |βbk,j | < c



2
τn,p /βbk,j
,

if |βbk,j | ≥ c


nuk −1/4
τn,p

(10)


nuk −1/4
,
τn,p

for some constant c and some arbitrary sequence uk with νk∗ ≤ uk ≤ νk∗ . We note that the
strength of the correlation between the variables in the model k affects the behavior of uk , and
(nuk /τn,p )−1/4 converges to zero as n tends to infinity due to Assumption 4 described in Section

17

4.
On the other hand, the penalty term in the other Bayesian model selection approaches is quite
different from that of the nonlocal priors as in (10). The marginal likelihood based on the g-prior
when σ 2 is known can be expressed as

l(βbk ) − |k| log(1 + g)/2.

Narisetty and He(2014) demonstrated that BASAD achieves the strong model selection consistency. This consistency follows from that the fact that the BASAD “penalty” is asymptotically
equivalent to

l(βbk ) − c|k| log(p),

(11)

where c is some constant. Yang et al. (2016) and Castillo et al. (2012) also considered a similar
penalty term on the model space, which implies that the posterior probability for their procedures
can be expressed in the same form as (11). When g = p2c , the marginal likelihood based on a
g-prior is asymptotically equivalent to (11).
The asymptotic term of the marginal likelihoods is quite different from that of the nonlocal priors, since the penalty terms in the other Bayesian approaches only focus on the model
size without considering the different weights on variables in the model. The marginal likelihoods based on nonlocal priors, however, impose different penalties on each predictor in the given
model. When the MLE of the regression coefficient in the model is asymptotically close to zero
(|βbk,j | < c(nuk /τn,p )−1/4 ), the model that contains the corresponding variable would be strongly
penalized by (nτn,p uk )1/2 . In contrast, when the MLE is asymptotically significant (|βbk,j | ≥
2
c(nuk /τn,p )−1/4 ), the penalty attains a different weight based on the MLE (pτn,p (βbk,j ) ≈ τn,p /βbk,j
).

This analysis highlights the fact that the nonlocal priors are able to adapt their penalty for the
inclusion of covariates based on the observed data, whereas the local priors must instead rely on a
prior penalty on non-sparse models.
18

7

Computational strategy

In p  n settings, full posterior sampling using existing Markov chain Monte Carlo (MCMC) algorithms is highly inefficient and often not feasible from a practical perspective. To overcome this
problem, we propose a scalable stochastic search algorithm aimed at rapidly identifying regions of
high posterior probability and finding the maximum a posteriori (MAP) model. Our main innovation is to develop a stochastic search algorithm combining isis-like screening techniques (Fan and
Lv, 2008) and temperature control that is commonly used in global optimization algorithms such
as simulated annealing (Kirkpatrick and Vecchi, 1983).
b that can be expressed as
To describe our proposed algorithm, consider the MAP model k

b = argmax{π(k | y)},
k

(12)

k∈Γ∗

where Γ∗ is the set of all models assigned non-zero prior probability.

7.1

Shotgun stochastic search algorithm

Hans et al. (2007) proposed the shotgun stochastic search (SSS) algorithm in an attempt to efficiently navigate through very large model spaces and identify global maxima. Letting nbd(k) =
{Γ+ , Γ− , Γ0 }, where Γ+ = {k∪{j} : j ∈ kc }, Γ− = {k\{j} : j ∈ k}, and Γ0 = {[k\{j}]∪{l} :
l ∈ kc , j ∈ k}, the SSS procedure is described in Algorithm 1.
Algorithm 1 Shotgun Stochastic Search (SSS)
Choose an initial model k(1)
For i = 1 to i = N − 1
Compute π(k | y) for all k ∈ nbd(k(i) )
Sample k+ , k− , and k0 , from Γ+ , Γ− , and Γ0 , with probabilities proportional to π(k | y)
Sample k(i+1) from {k+ , k− , k0 }, with probability proportional to
{π(k+ | y), π(k− | y), π(k0 | y)}

The MAP model can be identified by the model that achieves the largest (unnormalized) pos19

terior probability among those models searched by SSS.

7.2

Simplified shotgun stochastic search algorithm with screening (S5)

SSS is effective in exploring regions of high posterior model probability, but its computational cost
is still expensive because it requires the evaluation of marginal probabilities for models in Γ+ , Γ− ,
and Γ0 at each iteration. The largest computational burden occurs for the evaluation of marginal
likelihood for models in Γ0 , since |Γ0 | = |k|(p − |k|). To improve the computational efficiency
of SSS, we propose a modified version which only examines models in Γ+ and Γ− , which have
cardinality p − |k| and |k|, respectively. However, by ignoring Γ0 in the sampling updates we
make the algorithm less likely to explore “interesting” regions of high posterior model probability,
and therefore more likely to get stuck in local maxima. To counter this problem, we introduce a
“temperature parameter” analogous to simulated annealing which allows our algorithm to explore
a broader spectrum of models.
Even though ignoring models in Γ0 reduces the computational burden of the SSS algorithm, the
calculation of p posterior model probabilities in every iteration is still computationally prohibitive
when p is very large. To further reduce the computational burden, we borrow ideas from the
Iterative Sure Independence Screening (isis; Fan and Lv (2008)) and consider only those variables
which have a large correlation with the residuals of the current model. More precisely, we examine
the products |rkT Xj |, where rk is the residual of the model k, for j = 1, . . . , p, after every iteration
of the modified shotgun stochastic search algorithm, and then restrict attention to variables for
which {|rkT Xj | : j = 1, . . . , p} is large (we assume that the columns of X have been standardized).
This yields a scalable algorithm even when the number of variables p is large.
With these ingredients, we propose a new stochastic model search algorithm called Simplified
Shotgun Stochastic Search with Screening (S5), which is described in Algorithm 2.
In S5, Sk is the union of variables in k and the top Mn variables, obtained by screening using
the residuals from the model k. The screened neighborhood of the model k can be defined as
−
+
c
nbdscr (k) = {Γ+
scr , Γ }, where Γscr = {k ∪ {j} : j ∈ k ∩ Sk }.

20

Algorithm 2 Simplified Shotgun Stochastic Search with Screening (S5)
Set a temperature schedule t1 > t2 > . . . > tL > 0
Choose an initial model k(1,1) and a set of variables after screening Sk(1,1) based on k(1,1)
For l = 1 in l = L
For i in 1, . . . , J − 1
Compute all π(k | y) for all k ∈ nbdscr (k(i,l) )
−
1/tl
Sample k+ and k− , from Γ+
scr and Γ , with probabilities proportional to π(k | y)
Sample k(i+1,l) from {k+ , k− }, with probability proportional to {π(k+ | y)1/tl , π(k− | y)1/tl }
Update the set of considered variables Sk(i+1,l) to be the union of variables in k(i+1,l) and
the top Mn variables according to {|rkT(i+1,l) Xj | : j = 1, . . . , p}

Even though this algorithm is designed to identify the MAP model, it also provides an approximation to the posterior model probability of each model. The uncertainty of the model space
can be measured by approximating the normalizing constant from the (unnormalized) posterior
probabilities of the models explored by the algorithm.
Denoting the computational complexity of the evaluation of the unnormalized posterior model
probability of the largest model among searched models by En , the computational complexity
of the SSS algorithm can be expressed as the product of the number of explored models by the
algorithm and En , which is [O{N p}+O{N qn }+O{N (p−qn )qn }]×En , where qn is the maximum
size of model among searched models and qn < n  p.
On the other hand, the S5 only considers Mn variables after the screening step in each iteration,
which dramatically reduces the number of models to be considered in constructing the neighborhood, O{JL(Mn − qn )} + O(JLMn ). Therefore, the resulting computational complexity is

[O{JL(Mn − qn )} + O(JLMn )] × En + O(JLnp),

where qn < Mn . When the computational complexity for screening steps, O(JLnp), is dominated
by the other terms, the computational complexity is almost independent of p. As a result, the
proposed algorithm is scalable in the sense that the resulting computational complexity is typically
robust to the size of p.

21

7.3

Performance comparisons between S5 and SSS

We examined the computational performance of S5 to SSS in identifying the MAP model under
a piMOM prior with τn,p = log n log p and r = 1. We generated data according to Case (1) in
Section 3, with a fixed sample size (n = 200), and a varying number of covariates p. We set
Mn = 20, L = 20, and J = 20 for S5. To match the total number of iterations between S5 and
SSS, we set N = 400 for SSS. All computations were implemented in R.
Figure 4 shows the average computation time and the number of models searched before hitting
the MAP model for the first time for the S5 and SSS algorithms. All averages were based on 100
simulated datasets, and both algorithms obtained the same MAP model for all data sets. Panel (a)
shows that the computation time of SSS increases roughly at a p2 rate, but that the computation time
for S5 was nearly independent of the number of covariates p (about 4 seconds). For example when
p = 2, 000, SSS first found the MAP model in an average of 1,360 seconds (about 23 minutes),
whereas S5 hit the MAP model after about only 4 seconds. Interestingly, panel (b) of Figure 4
also illustrates that the S5 algorithms explored only 181 models on average to hit the MAP model,
whereas SSS typically visited slightly more than 38,000 models. Thus, not only is S5 much faster
than SSS in identifying the MAP model, but it also visited far fewer models before visiting the
MAP model.

8
8.1

Real data analysis
Analysis of polymerase chain reaction (PCR) data

Lan et al. (2006) studied coordinated regulation of gene expression levels on 31 female and 29
male mice (n = 60). A number of psychological phenotypes, including numbers of stearoyl-CoA
desaturase 1 (SCD1), glycerol-3-phosphate acyltransferase (GPAT) and phos- phoenopyruvate carboxykinase (PEPCK), were measured by quantitative real-time RT-PCR, along with 22,575 gene
expression values. The resulting data set is publicly available at http://www.ncbi.nlm.

22

38561.1

11

100

200

500

1000

4447.1

# of models

1636

8
7

log (# of models)

181.3

5

1

3.9

2

6

9

9

10

1359.7
60.3 148.4

seconds

5
4
3

log (seconds)

6

403.4

7

S5
SSS

2000

100

200

500

p

p

(a)

(b)

1000

2000

Figure 4: (a) the average computation time to first hit the MAP model; (b) the average number of
models searched before hitting the MAP model. The left y-axis is in a logarithmic scale and the
right y-axis is in the raw scale.
nih.gov/geo (accession number GSE3330).
Zhang et al. (2009) used penalized orthogonal components regression to predict the three phenotypes mentioned above based on the high-dimensional gene expression data. Bondell and Reich
(2012) also used the same data set to examine their model selection procedure based on penalizing regression coefficients within a (marginal or joint) credible interval obtained from a ridge-type
prior. For brevity, we restrict attention here to SCD1 as the response variable.
Since the ground truth regarding the true significant variables is not known for this data, we
compared our approach with a host of competitors on predictive accuracy and parsimony of the
selected model.
Prior to analyses, we standardized the covariates and randomly split the data set into 5 test samples and 55 training samples to evaluate the out-of-sample mean square prediction error (MSPE)

MSPE =

X

(yi − XiT βbkbtr )2 /|Ttest |,

i∈Ttest

where Ttest is the index set of the test samples and βbkbtr is the least square estimator under the
23

b based on the training samples. To avoid sensitivity to a particular split, we
estimated model k
considered 100 replications of the training and test sample generation. To measure the stability of
model selection, we considered the number of variables that were (i) selected at least 95 times, and
(ii) at least once, out of the 100 replicates.
Due to the high-computational burden of the penalized credible interval (Bondell and Reich,
2012) approach, we followed the pre-processing step suggested in their article to marginally screen
variables to reduce to 2000 variables (1999 genes and gender). For all the other approaches, all
22,575 genes were used. For the nonlocal priors, we considered both the MAP estimator and the
least squares (LS) estimator from the MAP model. For the g-prior, we set g = p2 as recommended
in George and Foster (2000). For the penalized likelihood procedures, we used ten-fold cross
validation to choose the tuning parameter.
To choose the hyperparameter τn,p for the nonlocal priors, we used a procedure proposed by
Nikooienejad et al. (2016). That procedure sets the hyperparameter so that the L1 distance between
the posterior distribution on the regression parameters under the null distribution (i.e., β = 0)
and the nonlocal prior distributions on these parameters is constrained to be less than a specified
value (e.g., p−1/2 ). The average value of the hyperparameter values chosen by this procedure were
τn,p = 1.12 and τn,p = 1.16 for piMoM and peMoM priors, respectively.
To make the comparison between the nonlocal priors and the g-prior more transparent, we used
the same beta-binomial prior on the model space in both models, rather than the uniform prior on
the model space described previously. The form of the beta-binomial prior was given by

π(k) ∝ ρ|k| (1 − ρ)p−|k| I(|k| ≤ qn ),

(13)

with a uniform prior on ρ and qn = 40. We note that this prior does not strongly induce sparsity as
does, for example, the prior obtained by imposing a Beta(1, pu ), u > 1 prior on ρ, as suggested in
Castillo et al. (2015).

24

Method
piMoM(MAP)
piMoM(LS)
peMoM(MAP)
peMoM(LS)
g-prior
lasso
scad
mcp
Marginal(p = 2000)
Joint(p = 2000)

MSPE
0.283 (0.17)
0.282 (0.17)
0.291 (0.18)
0.287 (0.17)
0.368 (0.20)
0.542 (0.39)
0.308 (0.23)
0.308 (0.21)
0.456 (0.40)
0.440 (0.40)

MS
FS
1.00 (0.00)
1
1.00 (0.00)
1
1.02 (0.14)
1
1.02 (0.14)
1
4.07 (0.56)
1
17.97 (8.62)
1
12.66 (7.62)
2
2.20 (0.94)
0
17.47 (11.16) 0
16.42 (11.06) 1

TS
1
1
2
2
133
211
163
29
273
185

Table 2: Analysis of the PCR data. Marginal and Joint refer to the variable selection procedures
(Bondell and Reich, 2012) based on Bayesian marginal credible set and Bayesian joint credible
set, respectively. MS is the average size of the selected model. FS is the number of frequently
selected variables, i.e., that were selected at least 95 times in 100 repetitions. TS refers to the total
number of variables selected at least once from 100 repetitions. Standard errors are provided in
parenthesis.
Table 2 summarizes the results from the analysis of the gene expression data set. On average,
the nonlocal priors simultaneously produced the lowest MSPE and the most parsimonious model.
The other model selection methods selected a wide array of different variables for different splits
of the data set. In particular, lasso and the penalized credible region approach selected more than
180 different variables from 100 repeated splits, while the average size of the selected model was
less than 20 and the number of frequently selected variables was only zero or one, indicating a
potentially large number of false positives picked up by these methods.

8.2

A simulation study based on the Boston housing data

We next examined the Boston housing data set that contains the median value of owner-occupied
homes in the Boston area, together with several variables that might be associated with their median
value. There were n = 506 median values in the data set, and we considered 10 continuous
variables as the predictor variables: crim, indus, nox, rm, age, dis, tax, ptratio, b, and
lstat. This data set has been used to validate a variety of approaches; some recent examples
relevant to variable selection include Radchenko et al. (2011), Yuan and Lin (2005), and Rockova

25

and George (2014).
To examine the model selection performance in high-dimensional settings, we added 1,000
noise variables that were generated independently from a standard Gaussian distribution (p =
1, 010). The same competitors from the previous subsection were used with the aforementioned
choice of hyperparameters. For nonlocal priors, the hyperparameter was chosen by the aforementioned procedure (Nikooienejad et al., 2016); the average of the chosen hyperparameter values
were τn,p = 2.01 and τn,p = 0.47 for piMoM and peMoM priors, respectively. Prior to analyses,
we standardized the covariates and considered a simulation test size of 100 samples.
Methods
MSPE
MS-O
MS-N
FS-O
piMoM(MAP) 24.281 (9.01) 5.05 (0.22) 0.01 (0.10)
5
piMoM(LS)
24.265 (9.04) 5.05 (0.22) 0.01 (0.10)
5
peMoM(MAP) 24.156 (9.02) 5.02 (0.14) 0.00 (0.00)
5
peMoM(LS)
24.165 (9.00) 5.02 (0.14) 0.00 (0.00)
5
g-prior
26.314 (9.87) 3.10 (0.44) 0.00 (0.00)
3
lasso
30.243 (11.82) 5.07 (0.87) 7.77 (11.16)
4
scad
33.993 (10.66) 5.39 (0.57) 31.60 (28.28)
5
mcp
26.191 (9.87) 4.66 (0.74) 0.54 (1.04)
3
Marginal
26.612 (10.16) 3.74 (0.88) 0.41 (0.72)
3
Joint
26.385 (10.25) 3.77 (0.94) 0.02 (0.20)
3

TS-O
6
6
6
6
5
8
7
6
7
6

Table 3: The Boston Housing data set: MS-O and MS-N refer to the average number of selected
original variables and selected noise variables, respectively. FS-O is the number of original variables that are frequently selected at least 95 times out of 100 repetitions. TS-O refers to the number
of original variables selected at least once from 100 repetitions.
The results of are analysis are summarized in Table 3. The conclusions are similar to those
reported in Section 8.1; the nonlocal priors consistently choose more parsimonious models and
had better predictive performance. The model selection procedure resulting from the nonlocal prior
selects almost the same variables across the 100 repetitions. The average number of the original
variables selected more than 95 times over 100 repetitions is 5, which is close to the average model
size. It is also reliable in the sense that the average number of the original variables that are selected
at least once across the repetitions is only 6. This means that model selection based on the nonlocal
prior selects the same model in most data splits. On the other hand, penalized likelihood methods
such as lasso and scad tend to select a large number of noise variables.
26

9

Conclusion

This article describes theoretical properties of peMoM and piMoM priors for variable selection
in ultrahigh-dimensional linear model settings. In terms of identifying a “true” model, selection
procedures based on peMoM priors are asymptotically equivalent to piMoM priors in Johnson and
Rossell (2012) because they share the same kernel, exp{−τn,p /β 2 }. We demonstrated that model
selection procedures based on peMoM priors and piMoM priors achieve strong model selection
consistency in p  n settings.
In Section 3.1, precision-recall curves were used to show that the model selection procedure
based on a g-prior can achieve nearly the same performance in identifying the MAP model as
nonlocal priors when an optimal value for the hyperparameter g is chosen. However, as shown
in Section 3.2, the value of the hyperparameter that maximizes the posterior probability of the
true model is very large and has high variability, which may limit the practical application of
this method. To overcome this problem, one can consider mixtures of g-prior as in Liang et al.
(2008), but the asymptotic behavior of Bayes factor and model selection consistency in ultrahighdimensional settings have not been examined for hyper-g priors, and they are difficult to implement
computationally.
In Section 7, we proposed an efficient and scalable model selection algorithm called S5. By
incorporating the SSS with a screening idea and a temperature control, S5 was able to accelerate
the computation speed without losing the capacity to explore the interesting region in the model
space. Under some simulation settings, it outperformed the SSS in a sense that not only did S5
search the MAP model much faster than the SSS, but it also found exactly the same MAP model
that was searched by the SSS.
Because the explicit form of the marginal likelihood of the nonlocal priors is not available,
we used the Laplace approximation throughout the paper, and Barber et al. (2016) studied the
accuracy of the approximation in Bayesian high-dimensional variable selection, especially when
the dimension of the approximation (which is qn ) and n are both increasing. However, their results
do not apply to the case of the nonlocal priors, since the nonlocal priors violate their regularity
27

condition (nonzero density at the origin). While empirical results in this paper and Johnson and
Rossell (2012) suggest that the use of the Laplace approximation is reasonable, in future work
it is still worth paying attention to the approximation error of the Laplace approximation to the
marginal likelihood of the nonlocal priors.
The close connection between our methods and the reduced rlasso procedures provides a useful contrast between Bayesian and penalized likelihood methods for variable selection procedures.
According to the evaluation criteria proposed in Section 5, the two classes of methods appear to
perform quite similarly. A potential advantage of the reduced rlasso procedure, and to the lesser
extent the rlasso procedure, is reduced computation cost. This advantage accrues primarily because
the reduced rlasso can be computed from the least squares estimate of each model’s regression parameter, whereas the Bayesian procedures require numerical optimization to obtain the maximum
a posteriori estimate used in the evaluation of the Laplace approximation to the marginal density of
each model visited. However, the procedures used to search the model space, given the value of a
marginal density or objective function, are approximately equally complex for both classes of procedures. There are also potential advantages of the Bayesian methods. For example, it is possible
to approximate the normalizing constant of the posterior model probability from the models visited
by S5 algorithm, and to use this normalizing constant to obtain an approximation to the posterior
probability assigned to each model. In so doing, the Bayesian procedures provide a natural estimate of uncertainty associated with model selection. These posterior model probabilities can also
be used in Bayesian modeling averaging procedures, which have been demonstrated to improve
prediction accuracy (e.g., Raftery et al. (1997)) over prediction procedures based on maximum a
posteriori estimates. Finally, the availability of prior densities may prove useful in setting model
hyperparameters (i.e., τn,p ) in actual applications, where scientific knowledge is typically available
to guide the definition of the magnitude of substantively important regression parameters.
We also developed an R package BayesS5 that provides all computational functions used in
this paper, including a support of parallel computing environments. It is available on the author’s
website and on CRAN.

28

10

Supplemental Materials

The supplemental material contains proofs of the technical results stated in the paper and the
Laplace approximations to evaluate the marginal likelihoods based on the nonlocal priors.

References
Barber, R. F., Drton, M., and Tan, K. M. (2016). Laplace approximation in high-dimensional
Bayesian regression. In Statistical Analysis for High-Dimensional Data, pages 15–36. Springer.
Bondell, H. and Reich, B. (2012). Consistent high-dimensional Bayesian variable selection via
penalized credible regions. Journal of the American Statistical Association 107, 1610–1624.
Castillo, I., Schmidt-Hieber, J., Van der Vaart, A., et al. (2015). Bayesian linear regression with
sparse priors. Annals of Statistics 43, 1986–2018.
Castillo, I., van der Vaart, A., et al. (2012). Needles and straw in a haystack: Posterior concentration for possibly sparse sequences. Annals of Statistics 40, 2069–2101.
Chen, J. and Chen, Z. (2008). Extended Bayesian information criteria for model selection with
large model spaces. Biometrika 95, 759–771.
Davis, J. and Goadrich, M. (2006). The relationship between Precision-Recall and ROC curves. In
Proceedings of the 23rd international conference on Machine learning, pages 233–240. ACM.
Fan, J. and Li, R. (2001). Variable selection via nonconcave penalized likelihood and its oracle
properties. Journal of the American Statistical Association 96, 1348–1360.
Fan, J. and Lv, J. (2008). Sure independence screening for ultrahigh dimensional feature space.
Journal of the Royal Statistical Society: Series B 70, 849–911.
George, E. and Foster, D. P. (2000). Calibration and empirical Bayes variable selection. Biometrika
87, 731–747.
29

Hans, C., Dobra, A., and West, M. (2007). Shotgun stochastic search for “large p regression.
Journal of the American Statistical Association 102, 507–516.
Jiang, W. (2007). Bayesian variable selection for high dimensional generalized linear models:
Convergence rates of the fitted densities. Annals of Statistics 35, 1487–1511.
Johnson, V. E. and Rossell, D. (2010). On the use of non-local prior densities in Bayesian hypothesis tests. Journal of the Royal Statistical Society: Series B 72, 143–170.
Johnson, V. E. and Rossell, D. (2012). Bayesian model selection in high-dimensional settings.
Journal of the American Statistical Association 107, 649–660.
Kim, Y., Kwon, S., and Choi, H. (2012). Consistent model selection criteria on high dimensions.
Journal of Machine Learning Research 13, 1037–1057.
Kirkpatrick, S. and Vecchi, M. (1983). Optimization by simulated annealing. Science 220, 671–
680.
Lan, H., Chen, M., Flowers, J. B., Yandell, B. S., Stapleton, D. S., Mata, C. M., Mui, E. T.K., Flowers, M. T., Schueler, K. L., Manly, K. F., et al. (2006). Combined expression trait
correlations and expression quantitative trait locus mapping. PLoS Genet 2, e6.
Liang, F., Paulo, R., Molina, G., Clyde, M., and Berger, J. (2008). Mixtures of g priors for Bayesian
variable selection. Journal of the American Statistical Association 103,.
Liang, F., Song, Q., and Yu, K. (2013). Bayesian Subset Modeling for High-Dimensional Generalized Linear Models. Journal of the American Statistical Association 108, 589–606.
Narisetty, N. N. and He, X. (2014). Bayesian variable selection with shrinking and diffusing priors.
Annals of Statistics 42, 789–817.
Nikooienejad, A., Wang, W., and Johnson, V. E. (2016). Bayesian variable selection for binary
outcomes in high dimensional genomic studies using non-local priors. Bioinformatics 32, 1338–
1345.
30

Radchenko, P., James, G. M., et al. (2011). Improved variable selection with forward-lasso adaptive
shrinkage. Annals of Applied Statistics 5, 427–448.
Raftery, A. E., Madigan, D., and Hoeting, J. A. (1997). Bayesian model averaging for linear
regression models. Journal of the American Statistical Association 92, 179–191.
Rockova, V. and George, E. I. (2014). EMVS: The EM approach to Bayesian variable selection.
Journal of the American Statistical Association 109, 828–846.
Rossell, D. and Telesca, D. (2017+). Non-local priors for high-dimensional estimation. Journal of
the American Statistical Association .
Rossell, D., Telesca, D., and Johnson, V. E. (2013). High-dimensional Bayesian classifiers using
non-local priors. In Statistical Models for Data Analysis, pages 305–313. Springer.
Scott, J. and Berger, J. (2010). Bayes and empirical-Bayes multiplicity adjustment in the variableselection problem. Annals of Statistics 38, 2587–2619.
Song, Q. and Liang, F. (2015).

High dimensional variable selection with reciprocal L1 -

regularization. Journal of the American Statistical Association 110, 1602–1620.
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal
Statistical Society: Series B 58, 267–288.
Yang, Y., Wainwright, M. J., and Jordan, M. I. (2016). On the computational complexity of highdimensional bayesian variable selection. Annals of Statistics 44, 2497–2532.
Yuan, M. and Lin, Y. (2005). Efficient empirical Bayes variable selection and estimation in linear
models. Journal of the American Statistical Association 100, 1215–1225.
Zellner, A. (1986). On assessing prior distributions and Bayesian regression analysis with g-prior
distributions. In Bayesian inference and decision techniques: Essays in Honor of Bruno de
Finetti, pages 233–243. North Holland, Amsterdam.
31

Zhang, C.-H. (2010). Nearly unbiased variable selection under minimax concave penalty. Annals
of Statistics 38, 894–942.
Zhang, D., Lin, Y., Zhang, M., et al. (2009). Penalized orthogonal-components regression for large
p small n data. Electronic Journal of Statistics 3, 781–796.
Zou, H. (2006). The adaptive lasso and its oracle properties. Journal of the American Statistical
Association 101, 1418–1429.

32

Supplemental Materials to “Scalable Bayesian Variable
Selection Using Nonlocal Prior Densities in
Ultrahigh-dimensional Settings”

1

Preliminary Results

Lemma 1. For Qk defined in (6),

Qk

j=1

QLk,j ≤ Qk ≤

Qk

j=1

QUk,j ,

where

∗2
QLk,j = c1 (σ 2 )1/2 (nνk∗ + 1/τn,p )−1/2 exp{−τn,p /βek,j
},

QUk,j = c2 (σ 2 )1/2 (nνk∗ + 1/τn,p )−1/2 exp{−τn,p /(|βek,j | + e
n )2 },
∗
∈ [βek,j − e
n , βek,j + e
n ] \ (−e
n , e
n )c for some positive constants
and e
n  (nνk∗ /τn,p )−1/4 , with βek,j

c1 and c2 .
e k = (X T Xk + 1/τn,p Ik )−1 . From (8), all eigenvalues of (Σ
e k )−1 are bounded
Proof. Recall Σ
k
between nνk∗ + 1/τn,p and nνk∗ + 1/τn,p , which implies for all x ∈ R|k| , (nνk∗ + 1/τn,p )xT x ≤
e k )−1 x ≤ (nν ∗ + 1/τn,p )xT x. Let T1n = {(nν ∗ + 1/τn,p )/σ 2 }1/2 and T2n = {(nνk∗ +
xT (Σ
k
k
1/τn,p )/σ 2 }1/2 . Substituting the above inequality in the expression for Qk , we have
|k|
Y

g1 (βek,j ) ≤ Qk ≤

j=1

|k|
Y

g2 (βek,j ),

(S1)

j=1

where
Z

∞

gi (βek,j ) =

2
2
}dβk,j ,
exp{−Tin
(βk,j − βek,j )2 /2 − τn,p /βk,j

(S2)

−∞

for i = 1, 2. We establish the lower bound first by showing that g1 (βek,j ) ≥ QLk,j for all j =

1

1, . . . , |k|. Recall e
n  (nνk∗ /τn,p )−1/4 from the statement of the Lemma. We have
Z
g1 (βek,j ) ≥
[βek,j −e
n ,βek,j +e
n ]\(−e
n ,e
n )c
∗2
≥ exp{−τn,p /βek,j
}

2
2
}dβk,j
(βk,j − βek,j )2 /2 − τn,p /βk,j
exp{−T1n

Z
[βek,j −e
n ,βek,j +e
n ]\(−e
n ,e
n )c

2
exp{−T1n
(βk,j − βek,j )2 /2}dβk,j ,

∗
for some βek,j
∈ [βek,j − e
n , βek,j + e
n ] \ (−e
n , e
n )c . Then, the integral in the last line of the above

display is equivalent to
Z

2 t2 /2
−T1n

e

dt ≥

[−e
n ,e
n ]\(−βek,j −e
n ,−βek,j +e
n )c

−1
c1 T1n

Z

T1n e
n

e−z

2 /2

−1
dz ≥ c2 T1n
,

0

where c1 and c2 are some positive constants and the last inequality in the above display follows
−1
∗2
since T1ne
n ≥ 1 for large n. Substituting back in the previous display, g1 (βek,j ) ≥ c1 T1n
exp{−τn,p /βek,j
}

for some constant c1 > 0, completing the proof of the lower bound.
We now establish the upper bound by showing that g2 (βek,j ) ≤ QUk,j for all j = 1, . . . , |k|. It
is straightforward to see that g2 is a symmetric function (i.e, g2 (βek,j ) = g2 (|βek,j |)), so that it is
enough to establish the bound for βek,j > 0; without loss of generality we assume that βek,j > 0.
We have
Z

∞
2
2
exp{−T2n
(βk,j − βek,j )2 /2 − τn,p βk,j
}dβk,j

−∞
0

Z
=

2
2
}dβk,j
exp{−T2n
(βk,j − βek,j )2 /2 − τn,p /βk,j

−∞
βek,j +e
n

Z
+

Z0 ∞
+

2
2
exp{−T2n
(βk,j − βek,j )2 /2 − τn,p /βk,j
}dβk,j

2
2
exp{−T2n
(βk,j − βek,j )2 /2 − τn,p /βk,j
}dβk,j .

βek,j +e
n

Define the first term of the above as W1 , the second as W2 , and the third term as W3 . First, we
−1
shall show that W1 ≤ cT2n
exp{−T2n (2τn,p )1/2 } for some positive constant c. By transforming

2

the variable t = βk,j − βek,j ,
Z

0
2 e2
2 e
2 2
βk,j /2 − τn,p /t2 }dt
tβk,j − T2n
t /2 + T2n
exp{−T2n

W1 =
−∞
0

Z

2 2
exp{−T2n
t /2 − τn,p /t2 }dt

≤

−∞
−1
c3 T2n

≤

for some constant c3 , since

exp{−T2n (2τn,p )1/2 },

exp{−µ/t2 − ζt2 }dt = (π/ζ)−1/2 exp{−2(µζ)1/2 } for µ > 0 and

R

ζ > 0.
Second, by changing the variable z = t − e
,
Z

βek,j

2
exp{−T2n
(z − βek,j + e
n )2 /2 − τn,p /(z + e
n )2 }dz
−e
n
Z ∞
2
2
≤ exp{−τn,p /(βek,j + e
n ) }
exp{−T2n
(z − βek,j + e
n )2 /2}

W2 =

−∞

≤

−1
c4 T2n

exp{−τn,p /(βek,j + e
n )2 },

for some positive constant c4 .
Third, by changing the variable z = t − βek,j , there exists some positive constant c such that
Z

∞

2 2
exp{−T2n
z /2 − τn,p /(z + βek,j )2 }dz
e
n
Z ∞
2 2
2 2
≤ exp{−T2ne
n /4}
exp{−T2n
z /4}dz

W3 =

≤

−1
c5 T2n

−∞
1/2
exp{−c6 T2n τn,p
},

for some constants c5 and c6 , since e
n  (nνk∗ /τn,p )−1/4 . Then,
−1
−1
g2 (βek,j ) ≤ c3 T2n
exp{−T2n (2τn,p )1/2 } + c4 T2n
exp{−τn,p /(βek,j + e
n )2 }
−1
1/2
+c5 T2n
exp{−c6 T2n τn,p
}.

1/2
Since e
n  (nνk∗ /τn,p )−1/4 , when βek,j < e
n , τn,p /(βek,j + e
n )2 < τn,p /(4e
2n )  T2n τn,p , and

3

1/2
2
when βek,j ≥ e
n , τn,p /(βek,j + e
n )2 ≤ τn,p /(4βek,j
) < T2n τn,p . In overall, the right-hand side

of the above display would be dominated by the second term, which shows that g2 (βek,j ) ≤
−1
cT2n
exp{−τn,p /(βek,j + e
n )2 } for some constant c. When βek,j < 0, we can show the same re-

sult by following exactly the same steps explained above.
We now present some auxiliary results that are used to prove Theorems 1 and 2. We make use of
the following simple union bound multiple times: for non-negative random variables V1 , . . . , Vm
and a > 0,

P(

m
X
l=1

Vl > a) ≤

m
X

P (Vl > a/m) ≤ m max P (Vl > a/m).
1≤l≤m

l=1

(S3)

We define some notations that are used in the subsequent proofs. Let t denote the true data
generating model, and let βt0 denote the true regression coefficient corresponding to t. Let ct =
t \ k, ck = k \ t, and u = k ∪ t. Also, we define the cardinality of a model k as k and in the
same spirit, denote ck = |ck |, ct = |ct |, and t = |t|. {x}j denotes the j-th element of the vector
x, and diag{A}j refers to the j-th diagonal element in the square matrix A. We denote χ2m (λ) a
non-central chi-square distribution with the degrees of freedom m and non-centrality parameter λ;
a central chi-square distribution is simply denoted by χ2m .
An important property that is used in the subsequent proofs concerns the distribution of the
marginal ridge estimator. Let βek = (XkT X + 1/τn,p Ik )−1 XkT y and βek,j = {βek }j . Then,
2∗
∗
),
βek,j ∼ N (βk,j
, σk,j

(S4)

∗
2∗
where βk,j
= {(XkT X + 1/τn,p Ik )−1 XkT Xt βt∗ }j and σk,j
= σ 2 diag{(XkT Xk + 1/τn,p Ik )−1 }j . It is
∗ 2
2∗
also evident that (βek,j − βk,j
) /σk,j
∼ χ21 .

4

A set of technical results follow that are used in the proof of the main results. Define

H1n =

X mk (y)π(k)
X π(k | y)
=
,
m
π(t
|
y)
t (y)π(t)
k:t(k,
k:t(k,
|k|≤qn

H2n =

X mk (y)π(k)
X π(k | y)
=
. (S5)
mt (y)π(t)
π(t | y)

k:t*k,
|k|≤qn

|k|≤qn

k:t*k,
|k|≤qn

Lemma 2. Fix  > 0. Let Γd = {k : |k| ≤ qn , t ( k, |k| − |t| = d} for d = 1, . . . , qn − |t|.

Suppose there exist constants c, δ > 0 such that maxk∈Γd P π(k | y)/π(t | y) > p−d /qn ≤
cp−d(1+δ) for d = 1, . . . , qn − |t|. Then, H1n converges to zero in probability as n tends to ∞,
where H1n is as in (S5).
Proof. Clearly, |Γd | =

p−|t|
d


. Using (S3), we bound

n −|t|
n X π(k | y)
o
n qX
o
X π(k | y)
P
>
= P
>
π(t | y)
π(t | y)
k:t(k
d=1 k∈Γ
d

qn −|t|

≤

X

P

d=1
qn −|t| 

≤

X
d=1

Finally,

Pqn −|t|
d=1

n X π(k | y)
o
> /qn
π(t | y)
k∈Γ
d


n −|t|
n π(k | y)
o qX
p − |t|
−d
maxP
cp−dδ .
> p /qn ≤
k∈Γd
d
π(t | y)
d=1

cp−dδ ≤ cqn p−δ → 0 as n → ∞.

Lemma 3. Fix  > 0 and let t = |t|. Define Γk,ck ,ct = {k : |k| ≤ qn , |k| = k, |k\t| = ck , |t\k| =
ct } for k = 0, . . . , qn ; ck = 0, . . . , k; ct = 1, . . . , t. Suppose

max P

k∈Γk,ck ,ct

h π(k | y)
π(t | y)

i
> n−3 p−k n−ck t−t ≤ cp−k(1+δ) ,

with some postive constants c and δ. Then, H2n converges to zero as n tends to ∞, where H2n is
as in (S5).

5

Proof. Clearly, |Γk,ck ,ct | =

P

p
k



k
ck



t
ct


.

qn
k X
t
n X π(k | y)
o
nX
X
X
> ≤ P
π(t | y)
k=1 c =0 c =1 k∈Γ
k:t*k

k

≤ P

t

k,ck ,ct

qn
k X
t
nX
X

X

k=1 ck =0 ct =1 k∈Γk,ck ,ct
qn
k X
t
n
X
X
P
≤

X

k=1 ck =0 ct =1

k∈Γk,ck ,ct

≤

qn
k X
t
X
X

≤

o
π(k | y)
>
π(t | y)
o
π(k | y)
−3
> n
π(t | y)

pk nck tt max P
k∈Γk,ck ,ct

k=1 ck =0 ct =1
qn
k X
t
X
X

o
π(k | y)
>
π(t | y)

n π(k | y)
π(t | y)

> n−3 p−k n−ck t−t

o

pk nck tt p−k(1+δ) → 0,

k=1 ck =0 ct =1

as n → ∞.
Lemma 4. Suppose W follows a non-central chi-square distribution with the degree of freedom
mn that is a positive integer and the non-central parameter λn ≥ 0, i.e, W ∼ χ2mn (λn ). Also,
consider wn and tn such that wn → 0 and tn → ∞ as n tends to ∞. Also, assume that mn ≺ tn .
Then,
2
P (W ≤ λn wn ) ≤ c1 λ−1
n exp{−λn (1 − wn ) },

(S6)

And

P (W > λn + tn ) ≤ c2

tn
2mn

mn /2
exp {mn /2 − tn /2} +

c3 λn1/2 t−1
n

t2n
exp −
32λn



,

(S7)

where c1 , c2 , and c3 are some positive constants.
Proof. W can be expressed as W =

Pmn

i=1 {Zi

i.i.d

+ (λn /mn )1/2 }2 , where Zi ∼ N (0, 1) for i =

1, . . . , m. Then, by the fact that P (Z > a) ≤ (2π)−1/2 a−1 exp{−a2 /2} for any a > 0, we can

6

show that there exist some positive constants c1 such that

P (W ≤ λn wn ) = P

mn
X

Zi2 + 2(λn /mn )1/2
mn
X


≤ P mn−1/2
= P

Zi + λn ≤ λn wn

i=1

i=1



mn
X

Zi ≤ −λn1/2 (1 − wn )/2

i=1
|Z1 | ≥ λn1/2 (1

− wn )/2 /2

2
≤ c1 λ−1
n exp{−λn (1 − wn ) /2},

since Z1 follows a standard normal distribution.
Also, by using Chernoffs’s bound and the fact that P (Z > a) ≤ (2π)−1/2 a−1 exp{−a2 /2} for
any a > 0, one can show that

P (W > λn + tn ) = P

(m
n
X

Zi2 + 2(λn /mn )1/2

mn
X

i=1
mn
X

≤ P

!
Zi2 > tn /2

mn−1/2

i=1


≤ c2

tn
2mn

Zi > tn

i=1

(
+P

)

mn
X

)
Zi > λn−1/2 tn /4

i=1

mn /2
exp {mn /2 − tn /2} +

c3 λn1/2 t−1
n



t2n
exp −
,
32λn

where c2 and c3 are some positive constants.
Lemma 5. Consider Qk defined in (6) for an arbitrary model k. Fix any δ > 0. For any k with
t ( k,



2/3
1−δ/8
P Qk /Qt > exp −|k \ t|τn,p
(nνk∗ )1/3 + |t|τn,p
(nνk∗ )δ/8 ≤ p−|k\t|(1+δ) ,

(S8)

and for k such that t * k,



P Qk /Qt > exp kβt0 k22 nνu∗ /{2 log(τn,p / log p)} ≤ p−|k|(1+δ) .

7

(S9)

Proof. By Lemma 1, it is sufficient to show that




#

"

Y
Y
2/3
1−δ/8
(nνk∗ )1/3 }
(nνk∗ )δ/8 } + P 
QUk,j > exp{−|k \ t|τn,p
(QUk,j /QLt,j ) > exp{|t|τn,p

P

j∈t

j∈k\t

≤ p−|k\t|(1+δ) .

(S10)

We first shall show that the first term in the left-hand side of (S10) is bounded above by
exp{−cnνk∗ } for some constant c.
"
P

=

X

≤

X

Y QUk,j
QLt,j
j∈t
" 
P c0

j∈t

#
1−δ/8
> exp |t|τn,p
(nνk∗ )δ/8



"

≤

X
j∈t

nνk∗ + 1/τn,p
nνt∗ + 1/τn,p

−1/2

∗
P [|βek,j − βk,j
| > 0 ] +

j∈t

X

 1−δ/8
QUk,j
P
>
exp
τn,p (nνk∗ )δ/8
QLt,j

#

n

o
 1−δ/8
∗2
exp −τn,p 1/(|βek,j | + e
n )2 − 1/βek,j
> exp τn,p
(nνk∗ )δ/8

∗
P [|βet,j − βt,j
| > 0 ],

j∈t

∗
for some small enough 0 > 0 and some positive constant c0 and βek,j
∈ [βek,j − e
n , βek,j + e
n ] \
∗
defined in (S4). The last inequality in the
(−e
n , e
n )c as defined in Lemma 1, and βek,j and βk,j

above display asymptotically holds, since

1−δ/8
∗
τn,p
(nνk∗ )δ/8  τn,p /(|βk,j
| − 0 − e
n )2 ,

for any δ > 0.
∗ 2
∗2
∗2
Since (βek,j − βk,j
) /σk,j
∼ χ21 and σk,j
≥ (nνk∗ + 1/τn,p )−1 , by using Lemma 4, one can show

that the first term in (S11) bounded above by exp{−c1 02 nνk∗ } for some constant c1 . Similarly, the
second term in (S11) is bounded above by exp{−c2 02 n} for some constant c2 , since Assumption
5 states that XtT Xt /n is asymptically isotropic. Therefore, (S11) is asymptotically bounded by
p−qn (1+δ) by Assumption 3.
Next, we shall show that the second term in the left-hand side of (S10) is bounded above by

8

#

(S11)

1/3

∗
exp{−cτn,p (nνk∗ )2/3 } for some positive constant c. Since when j ∈ k \ t and t ( k, βk,j
 n−1 ,


P


Y

2/3
QUk,j > exp{−|k \ t|τn,p
(nνk∗ )1/3 }

j∈k\t

"
X

≤

(

P c0 (nνk,j + 1/τn,p )−1/2 exp −

j∈k\t

X

=


P

2
βek,j

τn,p
(|βek,j | + e
n )2

#

)
2/3
(nνk∗ )1/3 }
> exp{−τn,p

n
o2 

1/2
1/3 2/3
0 −1/2
> τn,p (nνk∗ ) τn.p − log(nνk∗ + 1/τn,p )/2 + log c
−e
n

j∈k\t

"
X

≤

∗ 2
∗
P (βek,j − βk,j
) /σk,j
> c00

j∈k\t



τn,p
nνk∗

1/3

#
(nνk∗ + 1/τn,p )/σ 2 ,

∗ 2
∗
for some positive contant c0 and c00 . Since (βek,j − βk,j
) /σk,j
∼ χ21 , by Lemma 4 the last quantity
1/3

in the above display can be bounded by exp{−cτn,p (nνk∗ )2/3 } for some contant c. By Assumption
1/3

3, exp{−cτn,p (nνk∗ )2/3 } ≺ p−qn (1+δ) ≤ p|k\t|(1+δ)| , which proves the statement (S10).
We now shall show that the equation (S9) holds for any δ > 0. The left-hand side of (S9) can
be bounded above by
"
Y

P

QUk,j

j∈k

≤

X

Y

QLt,j

−1

#
> exp



kβt0 k22 nνu∗ /{2 log(τn,p / log p)}

j∈t

i
h
n
o

P c(nνk∗ + 1/τn,p )−1/2 exp −τn,p /(|βek,j | + e
n )2 > exp kβt0 k22 nνu∗ /{4|k| log(τn,p / log p)}

j∈k

+

X

h
n
o
i

∗2
P c0 (nνt∗ + 1/τn,p )1/2 exp τn,p /(βet,j
) > exp kβt0 k22 nνu∗ /{4|t| log(τn,p / log p)}

j∈t

"
≤

X

+

X

#
> kβt0 k22 nνu∗ /{4|k| log(τn,p / log p)} + log c

(S12)

h
i
∗
−1/2
1/2 1/2
P |βet,j
| < c00 kβt0 k−1
(nν
)
{4|t|
log(τ
/
log
p)}
τ
u∗
n,p
2
n,p ,

(S13)

P −

j∈k

τn,p
(|βek,j | + e
n )2

j∈t

where c, c0 , and c00 are some positive constants.
(S12) is always zero since the left-hand side in the probability is always negative and the righthand side in the probability operator is always positive. So, we focus on (S13) as below:
∗
∗
Since βet,j − e
n ≤ βet,j
≤ βet,j + e
n implies |βet,j | − e
n ≤ |βet,j
| ≤ |βet,j | + e
n , (S13) can be

9

bounded above by
X

h
i
−1/2
1/2 1/2
∗
(nν
)
{4|t|
log(τ
/
log
p)}
τ
| < c00 kβt0 k−1
P |βet,j
u∗
n,p
n,p
2

j∈t

≤

X

i
h
−1/2
1/2 1/2
+
e

,
(nν
)
{4|t|
log(τ
/
log
p)}
τ
P |βet,j | < c00 kβt0 k−1
n
u∗
n,p
n,p
2

j∈t

∗
2
2
∗2
2
2
where βt,j
is defined in (S4). Since βet,j
/σt,j
∼ χ21 (βt,j
/σt,j
) and σt,j
 σ 2 /n for j ∈ t, by

using Lemma 4 and Assumption 5, one can show that the probability is bounded by exp{−cn} for
some constant c, and it is evident that exp{−cn} ≺ p−|k|(1+δ) , which completes the proof of the
Lemma.

2

Proofs of Main Results

Proof of Theorem 1. We have π(t | y) = mt (y)π(t)/{

P

k:|k|≤qn

mk (y)π(k)}, since π(k) = 0

for any k with |k| > qn . Recall H1n and H2n from (S5) and note that π(t | y) = (1+H1n +H2n )−1 .
Hence to show that π(t | y) converges to one in probability, it is sufficient to establish that H1n
and H2n both converge in probability to zero as n tends to ∞. We shall prove the Theorem by
showing:
For any δ ∈ (0, 8/3) and any model k ∈ Γd (defined in Lemma 2),



π(k | y)
−d −1
P
> p qn ≤ p−d(1+δ) ,
π(t | y)

(S14)

and for any model k ∈ Γk,ck ,ct (defined in Lemma 3),



π(k | y)
−3 −k −ck −t
> n p n t
≤ cp−k(1+δ) .
P
π(t | y)

(S15)

Then, it is evident that H1n and H2n both converge to zero in probability by Lemma 2 and 3
respectively.

10

First, we shall show that (S14) holds. For any k ∈ Γd , recall that

P

h π(k | y)

i
h
n
i
o
1 e
−d Qk
−d
e
> p−d qn−1 ≤ P Cn,p
exp − 2 R
−
R
>
p
/q
k
t
n .
π(t | y)
Qt
2σ

ek > R∗ and R
et < R∗ + η, where η = d1 βbT βbt /τn,p for some constant d1 and βbt is the
Since R
t
t
k
ordinary least square estimator of βt in the true model t, by using (S3), the term in the last display
can be bounded above by
i
h


−d Qk
exp − Rk∗ − Rt∗ /(2σ 2 ) + η/(2σ 2 ) > p−d /qn
P Cn,p
Qt
h
i
−d Qk d(1+δ)+δ
≤ P Cn,p
p
> p−d /qn
Qt
 ∗

+P Rt − Rk∗ > 2σ 2 d(1 + δ) log p


+P exp{η/(2σ 2 )} > pδ .

(S16)
(S17)
(S18)

By using Lemma 5, (S16) is less than p−d(1+δ) when δ < 8/3. Since (Rt∗ − Rk∗ )/σ 2 ∼ χ2|k\t| , by
using (S6) in Lemma 4, we can show that (S17) is bounded by cp−d(1+δ) for some positive constant

c. Since τn,p nνt∗ η/d1 σ 2 ≤ βbtT XtT Xt βbt /σ 2 ∼ χ2|t| βt0T XtT Xt βt0 , by using the inequality (S7) in
Lemma 4, (S18) can be expressed as






P exp η/2σ 2 > pδ ≤ P τn,p nνt∗ η/d1 σ 2 > 2τn,p nνt∗ (log  + δ log p)/d1


≤ P βbtT XtT Xt βbt /σ 2 > 2τn,p nνt∗ (log  + δ log p)/d1
≤ (nδ log p)|t|/2 exp{−c1 δ(n log p)} + n−1/2 (δ log p)−1 exp{−c2 (n log p)2 /n}
≤ c3 p−|k|(1+δ) ,

(S19)

for some positive constant c1 , c2 , and c3 , which proves that (S14) holds.

11

Next, we consider (S15). Recall that u = k ∪ t. By using (S3), it can be shown that

P

h π(k | y)

> n−3 p−|k| n−|k\t| |t|−|t|

i

π(t | y)
h
n
o
i
−(|k|−|t|) Qk
ek − R
et )/(2σ 2 ) > n−3 p−|k| n−|k\t| |t|−|t|
≤ P Cn,p
exp −(R
Qt
h
i

−|k|−|t|) Qk
≤ P Cn,p
exp −(Rk∗ − Ru∗ )/(2σ 2 ) > n−3−|k\t| |t|−|t| p−|k|(2+δ)+δ
Qt


 ∗




+P exp Rt − Ru∗ /(2σ 2 ) ≥ p|k|(1+δ) + P exp η/(2σ 2 ) > pδ




≤ P exp Rt∗ − Ru∗ /2σ 2 > p|k|(1+δ)



+P exp η/2σ 2 > pδ


+P Rk∗ − Ru∗ < 2σ 2 kβt0 k22 nνu∗ / log(τn,p / log p)



+P Qk /Qt > exp kβt0 k22 nνu∗ /{2 log(τn,p / log p)} .

(S20)
(S21)
(S22)
(S23)

Since (Rt∗ − Ru∗ )/σ 2 follows a χ2|u\t| distribution, (S20) is also bounded by c1 p−|k|(1+δ) with
some constant c1 . By following the same steps regarding (S19), one can show that (S21) is
bounded by c2 p−|k|(1+δ) for some constant c2 . We note that (Rk∗ − Ru∗ )/σ 2 ∼ χ2|u\k| (λn ) with
λn = βt0T XtT (Pu − Pk )Xt βt0 , where Pk is the projection matrix of Xk . As discussed in Narisetty
and He (2014), λn ≥ nνu∗ kβt0 k22 . Hence, by using Lemma 4, one can show that (S22) is bounded
by exp{−c3 kβt0 k22 nνu∗ / log(τn,p / log p)} for some constant c3 . Lemma 5 states that (S23) is
bounded by p−|k|(1+δ) . In summary, since qn ≺ τn,p / log p by Assumption 3, there exists some
positive constant c4 such that P [π(k | y)/π(t | y) > n−3 p−|k| n−|k\t| |t|−|t| ] ≤ c4 p−|k|(1+δ) . which
completes the proof of Theorem 1.
Proof of Corollary 2. Recall the penalty term of a model k, Q∗k , based on the piMoM priors is

Q∗k

Z
=

2
b
exp − (βk − βbk )T Σ∗−1
k (βk − βk )/(2σ ) −



|k|
X
j=1

2
τn,p /βk,j

−r

|k|
X

2
log(βk,j
) dβk ,

j=1

 P

2
2
in (7). Since, for any  > 0, exp − |k|
j=1 {τn,p /βk,j + r log(βk,j )} is bounded above with respect
R
P|k|
2
2
b
to βk,j , Q∗k ≤ C exp{−(βk − βbk )T Σ∗−1
k (βk − βk )/(2σ ) −
j=1 (1 − )τn,p /βk,j }dβk for some
12

constant C. Following the exactly same steps in Lemma 1, Q∗k ≤ C 0 (nνk∗ )−1/2

Q|k|

j=1

exp{−(1 −

)τn,p /(|βbk,j | + e
n )2 } for some constant C 0 > 0.
We shall show that the model selection procedure based on piMoM priors as in (4) assures
consistency by proving that Q∗k and Qk are asymptotically equivalent.
Q
b∗2
Next, we shall show that Q∗k is bounded below by C(nνk∗ )−1/2 |k|
j=1 exp{−(1 − )τn,p /βk,j }

∗
2
2
for some constant C > 0 and βbk,j
∈ [βbk,j − e
n , βbk,j + e
n ]. Since exp − τn,p /βk,j
+ r log(βk,j
)
can be minimized in [βbk,j − e
n , βbk,j + e
n ], by following the proof of Lemma 1,
Z

∞

exp{−nνk∗ (β − βbk,j )2 /(2σ 2 ) − τn,p /β 2 − r log(β 2 )}dβ

−∞
βbk,j +e
n

Z

exp{−nνk∗ (β − βbk,j )2 /(2σ 2 ) − (1 − )τn,p /β 2 } exp{−τn,p /β 2 − r log(β 2 )}dβ
βbk,j −e
n
n
o
∗2
≥ C(nνk∗ )−1/2 exp −(1 − )τn,p /βbk,j
,

≥

∗
where C is some constant and βbk,j
∈ [βbk,j − e
n , βbk,j + e
n ] \ (−e
n , e
n )c .

Therefore, due to the asymptotic similarity between the ridge estimator and the least square
estimator, the lower and upper bounds of Q∗k are asymptotically equivalent to those of Qk with the
penalty parameter (1 − )τn,p , which assures the strong consistency of the model selection based
on the piMoM priors.
Proof of Theorem 3. Under a situation where σ 2 is unknown, it is clear that
|k|
− 2

mk (y) = τn,p

Z

(2πσ 2 )−

n+|k|
2

Z



|k|
  2 1/2 (β − βe )T Σ
−1
X
e
e
τn,p 
k
k
k (βk − βk )
exp |k|
−
−
π(σ 2 )dβk dσ 2 ,
2 

σ2
2σ 2
βk,j
j=1

where π(σ 2 ) is the prior for σ 2 (Inverse-gamma density with hyperparameters a0 and b0 ).
First, we shall show that the ratio between marginal likelihoods of a model k and the true model

13

t can be bounded as
|k|−|t|
mk (y)
≤c 2
mt (y)

ek + 2b0
R
et + 2b0
R

!−n/2−a0



|k|
|t|
 X
X
τn,p
τn,p  (nνk∗ τn,p + 1)−|k|/2
exp −
+
,
∗2  (nν ∗ τn,p + 1)−|t|/2

βet,j
(|βek,j | + e
n )2
t
j=1

j=1

(S24)

∗
∈ [βet,j − e
n , βet,j + e
n ] \ (−e
n , e
n )c for j ∈ 1, . . . , |t| and c is some constant. Next, we
where βet,j

ek − R
et )/(2σ02 (1 + un ))}, where σ02
ek + 2b0 )/(R
et + 2b0 )}−n/2−a0 ≤ exp{−(R
shall show that {(R
is the true regression variance that involves in the data-generating process, and un is some random
variable that is concentrated around a finite value with at least probability 1 − exp{−cn} for some
constant c. Then, by following the same steps in the proof of Theorem 1, the proof of Corollary 2
is completed.
By Lemma 1, the marginal likelihood of a model k can be bounded by

mk (y) ≤ {c1 (nνk∗ τn,p + 1)}−

|k|
2

≤ {c1 (nνk∗ τn,p + 1)}−

|k|
2

Z






1/2



1/2

|k|
X


e
Rk + 2b0  2
τn,p
−
dσ
2σ 2 
(|βek,j | + e
n )2

2
−

σ2
j=1


|k|
 X


− n+2a
0
τn,p
2
e
exp −
(1 + exp{2|k|}) Rk + 2b0
,
e

n )2 
j=1 (|βk,j | + e
(σ 2 )−

n+2a0
−1
2

exp

|k|

for some constant c1 .
Also, by using Lemma 1, one can show that

mk (y) ≥ {c2 (nνk∗ τn,p + 1)}−

≥ {c2 (nνk∗ τn,p + 1)}−

|k|
2

|k|
2

Z

(σ 2 )−

n+2a0
−1
2

exp




|k|

2
σ2


|k|
X
e
τn,p
Rk + 2b0  2
−
−
dσ
e∗2
2σ 2 
j=1 βk,j




|k|
 X
− n+2a
0
τn,p   e
2
Rk + 2b0
,
exp −

βe∗2 
j=1

k,j

∗
where c2 is some constant and βek,j
∈ [βek,j − e
n , βek,j + e
n ] \ (−e
n , e
n )c for j ∈ 1, . . . , |k|. These

results shows that (S24) holds.
ek +2b0 )/(R
et +2b0 )}−n/2−a0 in (S24). Define
Next, we consider the asymptotic behavior of {(R

14

ρn as the follows:
et + 2b0 )/(nσ 2 ) − 1.
ρn = (R
0
Since − log(1 − u) < u/(1 − u) for u ∈ R,
ek + 2b0 )/(R
et + 2b0 )} = − log[1 + (R
ek − R
et )/{n(1 + ρn )σ02 }]
− log{(R
et − R
ek )/{nσ 2 (1 + un )},
≤ (R
0
ek − R
et )/(nσ 2 ).
where un = ρn + (R
0
Since (Rk∗ − Ru∗ )/σ02 ∼ χ|u\k| (λn ) with λn = βt0T XtT (Pu − Pk )Xt βt0 /σ02 , by using Lemma 4
one can show that


P (|un − λn /n| > ) ≤ P (|ρn | > /4) + P (Rt∗ − Ru∗ )/(nσ02 ) > /4


+P (Rk∗ − Ru∗ )/(nσ02 ) − λn /n > /4 + P η/2nσ02 > /4

≤ exp{−c0 n} + P (Rk∗ − Ru∗ )/(nσ02 ) − λn /n > /4
≤ exp{−c00 n},

for some constant c0 and c00 , and η is defined in the proof of Theorem 1. Also, by Assumption 5,
λn /n will be bounded below and above.
Proof of Corollary 4. Since we showed that the asymptotic equivalence between Qk and Q∗k in
the proof of Corollary 2, by following exactly same steps in the proof of Theorem 3 we can prove
the model selection consistency under piMoM prior densities.
Proof of Proposition 5. We shall show that for any αk = βbk + n with n = {n,j }j=1,...,|k| and
|n,j |  ∗n for at least one j ∈ {1, . . . , |k|}, P {g(αk ; k) < g(βek∗ ; k)} → 0 as n tends to ∞, where
∗
βek∗ ∈ B(βbk ; ∗n ) with ∗n  (τn,p /n)1/3 . More specifically, we set βek,j
= βbk,j + ∗n for j ∈ t and
∗
βek,j
= βbk,j for j ∈ tc . Without loss of generality, we assume that XjT Xj = n for j = 1, . . . , p.

15

Note that

g(αk ; k) = ||Xk αk −

Xk βbk ||22

+

|k|
X

τn,p /|αk,j | + Dn

j=1
|k|

=

X

{cj n2n,j + τn,p /|βbk,j + n,j |} + Dn ,

j=1

for some constants cj such that CL < cj < CU for j = 1, . . . , |k|, and some randome variable Dn
that are not relevant to αk . Then,
P {g(αk ; k) < g(βek∗ ; k)}

(
)
(
)
|k|
|k|
X
X
τn,p
τn,p 
≤ P
cj n2n,j +
<
cj n∗2
n +
∗
|βbk,j + n,j |
|βek,j
|
j=1
j=1

)
)
(
(
X
X
τn,p 
τn,p
≤ P
− tn,j <
(S25)
cj n∗2
cj n2n,j +
n +
bk,j | + |n,j |
e∗ |
|
β
|
β
∗
∗
k,j
j∈S ∩Sk,n
j∈S ∩Sk,n

)
)
(
(
X
X
τn,p
τn,p 
+P 
− tn,j <
(S26)
cj n2n,j +
cj n∗2
n +
bk,j | + |n,j |
e∗ |
|
β
|
β
c
c
∗
∗
k,j
j∈S ∩Sk,n
j∈S ∩Sk,n
(
)
"
(
)#
X
X
X
τ
t
τ
n,p
n,j
n,p
cj n∗2
<
+P
cj n2n,j +
+
, (S27)
n +
∗c |
b
e
|S
|βk,j | + |n,j |
|β ∗ |
∗c
∗
∗c
j∈S

j∈S

j∈S

k,j

1/3

where tn is an arbitrary sequence such that tn,j = n2/3 τn,p n,j , and S ∗ = {j ∈ {1, . . . , p} : |n,j | 
∗n }, and Sk,n = {j ∈ k : |βbk,j | < ∗n }. Then, to complete the proof, it is sufficient to show that
each of (S25), (S26), and (S27) converges to zero.
0 2
Since n(βbk,j − βt,j
) /σ 2 ∼ χ21 for j = 1, . . . , |k|,

0
P (|βbt,j − βt,j
| > ζn ) ≤ (πnζn2 /2)−1/2 exp{−nζn2 /(2σ 2 )},

−1/2
2
for any ζn > 0. This implies that Sk,n = t at least probability 1−|tc |(πn∗2
exp{−n∗2
n /2)
n /(2σ )}.

16

Therefore, the equation (S25) can be asymptotically bounded by
"
X

P cj n2n,j

j∈S ∗ ∩t

X

≤

τn,p
τn,p
+
− tn,j < cj n∗2
n +
2|n,j |
|βbk,j + ∗n |

#

h
i
P |βbk,j + ∗n | < cτn,p (n2n,j − tn,j + τn,p /|n,j |)−1 ,

j∈S ∗ ∩t

2
2 2
∗2
2
for some positive constant c. Consider Lemma 4 with λn = n∗2
n /σ and wn = c τn,p /{n (nn,j −
2/3

tn,j + τn,p /|n,j |)2 } for j ∈ S ∗ ∩ t. Since n2n,j  n1/3 τn,p for j ∈ S ∗ implies wn → 0, Lemma 4
2
guarantees that the last display is bounded by c0 |S ∗ ∩ t|λ−1
n exp{−λn (1 − wn ) } for some constant

c0 , which means that (S25) converges to zero as n tends to 0. By following the same steps, one can
show that (S26) converges to zero.
Also, (S27) can be asymptotically bounded by
"
X

P cj n2n,j +

j∈S ∗c ∩t

τn,p
τn,p
+ c min∗ tn,j < cj n∗2 +
b
j∈S
2|n,j |
|βk,j + ∗n |

"
+

X

P cj n2n,j +

j∈S ∗c ∩tc

≤

X

τn,p
2|βbk,j + ∗n |

+ c min∗ tn,j < cj n∗2 +

P |βbk,j + ∗n | < c0 τn,p (n2n,j − n∗2
tn,j
n + c min
∗
j∈S

j∈S ∗c ∩t

+

j∈S



X

#

τn,p

#

|βbk,j + ∗n |

−1
+ τn,p /|n,j |)




∗
00
2
∗2
−1
b
P |βk,j + n | < c τn,p (nn,j − nn + c min∗ tn,j + τn,p /|n,j |) /2 ,
j∈S

j∈S ∗c ∩tc

where c, c0 , and c00 are some positive constants. For the first term in the last line of the above display,
2
2
∗
2
by setting λn = n∗2 /σ 2 and wn = c2 τn,p
/{∗2
n (nn,j − nn + c minj∈S ∗ tn,j + τn,p /|n,j |) }, we can
2
apply Lemma 4. Since wn ≺ τn,p
(∗n minj∈S ∗ tn,j )−2 implies wn → 0, the first term in the above

display converges to zero by Lemma 4. Similarly, the second term also converges to zero.

17

3

Laplace Approximations of Marginal Likelihoods

In this section, we provide the Laplace approximation of the marginal likelihoods based on the
nonlocal priors. Because closed form expressions for posterior model probabilities based on modified peMoM priors and modified piMoM priors are not available, we estimate the posterior model
probabilities using Laplace approximations. For posterior probabilities based on the peMoM priors, an inverse-Gamma density with parameters (a0 , b0 ) on σ 2 the Laplace approximation to the
marginal density of the data for model k can be expressed as

π(k | y) ∝ (2π)|k|/2 V (βk∗ , σ 2∗ )

−1/2

exp{f (βk∗ , σ 2∗ )}p(k),

(S28)

where

(βk∗ , σ 2∗ ) = argmaxf (βk , σ 2 )
(βk ,σ 2 )

f (βk , σ 2 ) = − (n/2 + |k|/2 + a0 + 1) log σ 2 − (y − Xk βk )T (y − Xk βk )/(2σ 2 ) − βkT βk /(2σ 2 τn,p )
−

|k|
X

2
τn,p /βk,j
+ |k|(2/σ 2 )1/2 − b0 /σ 2 + |k|(log τn,p )/2,

j=1

and V (βk , σ 2 ) is a (|k| + 1) × (|k| + 1) matrix with the following blocks:

4
V11 = XkT Xk /σ 2 + Ik /σ 2 τn,p + diag 6τn,p /βk,j

j=1,...,|k|

V12 = XkT (y − Xk βk )/σ 4 − βk /{σ 4 τn,p }
V22 = −(n/2 + |k|/2 + a0 + 1)/σ 4 + (y − Xk βk )T (y − Xk βk )/σ 6 − βkT βk /τn,p
−3|k|21/2 σ −5 /4 + 2b0 /σ 6 .

18

For the piMoM priors on βk , the Laplace approximation of the posterior model probability can be
expressed as in (S28), but with

f (βk , σ 2 ) = − (n/2 + a0 + 1) log σ 2 − (y − Xk βk )T (y − Xk βk )/(2σ 2 )
−

|k|
X



2
2
r log(βk,j
) + τn,p /βk,j
+ |k| (r − 1/2) log τn,p − log Γ(r − 1/2) − b0 /σ 2 ,

j=1

and V (βk , σ 2 ) a (|k| + 1) × (|k| + 1) matrix with the following blocks:

4
2
V11 = XkT Xk /σ 2 + diag 6τn,p /βk,j
− 2r/βk,j

j=1,...,|k|

V12 = XkT (y − Xk βk )/σ 4
V22 = −(n/2 + a0 + 1)/σ 4 + (y − Xk βk )T (y − Xk βk )/σ 6 + 2b0 /σ 6 .

19

