arXiv:1602.09100v2 [stat.ME] 3 Jul 2017

Bayesian Variable Selection for Skewed
Heteroscedastic Response
Libo Wang1 , Yuanyuan Tang1 , Debajyoti Sinha1 , Debdeep Pati1 , and Stuart Lipsitz2
1 Department of Statistics, Florida State University, Tallahassee, FL 32306
2 Brigham and Women’s Hospital

Abstract
In this article, we propose new Bayesian methods with proper theoretical justification for selecting and estimating a sparse regression coefficient vector for skewed
heteroscedastic response. Our novel Bayesian procedures effectively estimate the median and other quantile functions, accommodate non-local prior for regression effects
without compromising ease of implementation via sampling based tools. First time for
skewed and heteroscedastic response, this Bayesian method asymptotically selects the
true set of predictors even when the number of covariates increases in the same order
of the sample size. We also extend our method to deal with some observations with
very large errors. Via simulation studies and a re-analysis of a medical cost study with
large number of potential predictors, we illustrate the ease of implementation and other
practical advantages of our approach compared to existing methods for such studies.

Keywords: Bayesian consistency; median regression; sparsity

1

1

Introduction

Large number of possible predictors and highly skewed heteroscedastic response are often
major challenges for many biomedical and econometric applications. Selection of an optimal
set of covariates and subsequent estimation of the regression function are important steps
for scientific conclusions and policy decisions based on such studies. For example, previous
analyses of Medical Expenditure Panel Study (Natarajan et al., 2008; Cohen, 2003) testify
to the highly skewed and heteroscedastic nature of the main response of interest, total health
care expenditure in a year. Also, it is common in such studies to have a small proportion of
patients with either very high or very low medical costs. Popular classical sparse-regression
methods such as Lasso (Least absolute shrinkage operator) by Tibshirani (1996) and Efron
et al. (2004), and later related methods of Fan and Li (2001), Zou and Hastie (2005), Zou
(2006) and MCP (Zhang, 2010) assume Gaussian (or, at least symmetric) response density
with common variance. Limited recent literature on consistent variable selection for nonGaussian response includes Zhao and Yu (2006) under common variance assumption, Bach
(2008) under weak conditions on covariate structures, and Chen et al. (2014) under skew-t
errors. However, none of these methods deal with estimation of quantile function for heteroscedastic response frequently encountered in complex biomedical studies. Many authors
including Koenker (2005) argue effectively against focusing on mean regression for skewed
heteroscedastic response. Our simulation studies demonstrate that directly modeling skewness and heteroscedasticity, particularly in presence of analogous empirical evidence, leads
to better estimators of quantile functions for finite samples compared to existing methods
which ignore skewness and heteroscedasticity.
Bayesian methods for variable selection have some important practical advantages including incorporation of prior information about sparsity, evaluation of uncertainty about
the final model, interval estimate for any coefficient of interest, and evaluation of the relative
importance of different coefficients. Asymptotic properties of Bayesian variable selection
methods when the number of potential predictors, p, increases as a function of sample size n
have received lot of attention recently in the literature. Traditionally, to select the important
variables out of (X1 , . . . , Xp ), a two component mixture prior, also referred to as “spike and
slab” prior, (Mitchell and Beauchamp, 1988; George and McCulloch, 1993, 1997) is placed
on the coefficients β = (β1 , . . . , βp ). These mixture priors include a discrete mass, called a
“spike”, at zero to characterize the prior probability of a coefficient being exactly zero (that
is, not including the corresponding predictor in the model) and a continuous density called
a “slab”, usually centered at null-value zero, representing the prior opinion when the coefficient is non-zero. Following Johnson and Rossell (2010a,b), when the continuous density of
the slab part of a spike and slab prior has value 0 at null-value 0, we will call it a non-local
mixture prior. Continuous analogues of local mixture priors are being proposed recently
by Park and Casella (2008); Carvalho et al. (2010); Bhattacharya et al. (2014) among others. Bondell and Reich (2012) presented the selection consistency of joint Bayesian credible
sets. However, current Bayesian variable selection methods usually focus on mean regression
function for models with symmetric error density and common variance.
Johnson and Rossell (2010b) recently showed a startling selection inconsistency phenomenon for using several commonly used mixture priors, including local mixture (spike and
slab prior with non-zero value at null-value 0 of the slab density) priors, when p is larger than
2

√
the order of n. To address this for mean regression with sparse β, they advocated the use
of non-local mixture density presenting continuous “slab” density with value 0 at null-value 0
because these priors, called non-local mixture priors here, obtain selection consistency when
the dimension p is O(n). Castillo et al. (2014) provided several conditions to ensure selection consistency even when p  n. However, none of these Bayesian methods specifically
deal with skewed and heteroscedastic response, contamination of few observations with large
errors and variable selection for median and other quantile functionals.
In this article, we accommodate skewed and heteroscedastic response distribution using transform-both-sides model (Lin et al., 2012) with sparsity inducing prior for the vector
of regression coefficients. Our key observation is that, under such models with generalized
Box-Cox transformation (Bickel and Doksum, 1981), even a local mixture prior on aftertransform regression coefficients induces non-local priors on the original regression function
for certain choices of the transformation parameter. Similar to moment and inverse moment
non-local priors in Johnson and Rossell (2010b), this method enables clear demarcation between the signal and the noise coefficients in the posterior leading to consistent posterior
selection even when p = O(n). Addition to that, our use of standard local priors on the
transformed regression coefficients facilitates straightforward posterior computation which
can be implemented in publicly available softwares. We later extend this model to accommodate cases when the observations are contaminated with a small number of observations
with very large (or small) errors. Our approaches are shown to out-perform well-known
competitors in simulation studies as well as for analyzing and interpreting a real-life medical
cost study.

2
2.1

Bayesian variable selection model
Transform-both-Sides Model

For the skewed and heteroscedastic response Yi for i = 1, . . . , n, we assume the transformboth-sides model (Lin et al., 2012)
gη (Yi ) = gη (xTi β) + ei ,

(1)

where β = (β1 , . . . , βp )0 , xi is the observed p-dimensional covariate vector, gη (y) is the monotone power transformation (Bickel and Doksum, 1981),
y|y|η−1 − 1
,
gη (y) =
η

(2)

with unknown parameter η ∈ (0, 2). This transformation in (2) is an extension of Box-Cox
power family that has a long history and success in dealing with skewed and heteroscedastic response. We assume that under an optimal η, the transformed response gη (y) has a
symmetric and unimodal distribution with mean and median gη (xTi β). Thus ei ’s are independent mean 0 errors with common symmetric density function fe and variance σ 2 . The
transformation gη (y) in (2) is monotone with derivative gη0 (y) = |y|η−1 ≥ 0. Model (1) can
be expressed as a linear model
Yi = xTi β + i ,
(3)
3

where i has a skewed heteroscedastic density with median 0 because P [i > 0] = 1/2, and
approximate variance is σ 2 |xTi β|2−2η . Hence, the median of the skewed and heteroscedastic
response Yi in (1) is xTi β. For the time being, we consider a Gaussian N(0, σ 2 ) density for fe
in (1). Later in §4, we consider other densities to accommodate a heavy tail for fe .
For the model of (1), any sparsity inducing prior for β should depend on the transformation parameter η since η has a significant effect on the range and scale of Yi (approximate
variance σ 2 |xTi β|2−2η ). Based on this argument, we specify an conditional mixture prior
gη (βj ) given η using a “local” φ(·; 0, σβ2 ) density for the “slab” when gη (βj ) is non-zero with
discrete prior probability (1 − π0 ), where φ(·; µ, v 2 ) is the Gaussian density with mean µ and
variance v 2 . This conditional mixture prior gη (βj ) given η (a local mixture prior according
to definition of Johnson and Rossell (2010b)) results in a conditional mixture prior
fβ (βj | η) = π0 δ0 (βj ) + (1 − π0 )φ(gη (βj ); 0, σβ2 )|gη0 (βj )|

(4)

for βj given η independently for j = 1, · · · , p, where π0 ∈ [0, 1] is the probability of βj being
zero, δ0 (·) is the discrete measure at 0. When η > 1, gη0 (0) = 0 and hence the resulting
unconditional marginal prior for the nonzero βj in (4) turns out to be a non-local mixture
prior of Johnson and Rossell (2010b). However, the prior of transformed gη (βj ), a mixture of
discrete measure at 0 and φ(·; 0, σβ2 ) density, is a local mixture prior. This is demonstrated
via the plots of two resulting unconditional priors of βj when η = 0.5 and η = 1.5 in Figure 1.
Our simulation study in §5 shows that the model selection and estimation procedures for our
Bayesian method perform substantially better than competing methods when η > 1 (case
with non-local unconditional prior for βj ) compared to, say, when η = 0.5 ∈ (0, 1) (case with
a local prior for β given η). Thus, heteroscedasticity and the possibly non-local property
of π(β | η) come as a bi-product of the transform-both-sides model of (1). This implicit
non-local mixture prior modeling of unconditional βj may be the reason for some desirable
asymptotic properties of our method even when p = O(n) (discussed in §3). However, our
methods’ ability to use a local prior for gη (βj ) significantly reduces computational complexity
of the associated Markov chain Monte Carlo (MCMC) algorithms, while facilitating the
desirable asymptotic property.

Figure 1: Density Plot for Different η
When the error density fe in (1) is Gaussian φ(·; µ, σ 2 ) with mean 0 and variance σ 2 ,
we can specify a hierarchical Bayesian model using the prior π(β | η) in (1) along with known
4

priors for σ 2 and η,
σ 2 ∼ IGa(a, b),

η/2 ∼ Beta(c1 , d1 )

(5)

and a hyperprior for the hyperparameters (π0 , σβ2 ). For computational simplicity, in §3, we
establish variable selection consistency of this hierarchical Bayesian model using (1) along
with the non-local prior (4) on β and the prior for σ 2 in (5).

3

Consistent variable selection for large p

In this section, we investigate selection consistency of the proposed model when the number
of covariates p is grows with sample size n, with p ≤ n and the true error density is skewed and
heteroscedastic, but follows the same specification as in (1). Unlike the Gaussian likelihood
with mixture priors of Johnson and Rossell (2010b), our Bayesian model described in (1) - (5)
does not admit a closed form expression of the marginal likelihood. We derive appropriate
bounds of the marginal likelihood to obtain the desired Bayesian consistency results. Here
we only present a brief outline of our assumptions, developments and practical implications
of our theoretical results. Supporting results and details of the proofs are deferred to §7.
For brevity of exposition, we consider a design matrix X which are nearly orthogonal in the
sense that there exists constant 0 < c1 ≤ c1 such that c1 n ≤ λ1 (X T X) ≤ λp (X T X) ≤ c2 n
where λ1 (A) ≤ λ2 (A) ≤ · · · ≤ λp (A) denote ordered eigen values of the matrix A. The
assumption ensures identifiability of the regression coefficient β and is commonly assumed
in the selection consistency literature, refer for example to Johnson and Rossell (2010b).
To estimate the posterior probability assigned to the correct model, we use Laplace
approximation (Rossell and Telesca, 2015) to obtain probabilistic bounds of marginal likelihood p(y|γ). Here γ = (γ1 , . . . , γp ) denotes the predicted model P
for which γj = I(βj 6= 0) for
j = 1, 2, · · · , p are the indicators of active coefficients and pγ = pi=1 γi denotes the number
of active variables. Denote by β ∗ the vector of true regression coefficients with γ ∗ defined
accordingly. Assuming β̃ as the posterior mode and η is the optimal power parameter, the
Laplace approximation of p(y|γ) around β̃ is
p(y|γ) ≈

(2π)pγ /2
exp{log L(β̃) + log(π(β̃|γ)))}
|H(β̃)|1/2

(6)

where the projection of the prior fβ in (4) onto the support of γ, denote by π(β | γ) and is
given as
Y
π(β | γ) =
φ(gη (βj ); 0, σβ2 ) gη0 (βj ) .
(7)
j∈γ

We note again that (7) becomes a non-local prior (Johnson and Rossell, 2010b) for η ∈ (1, 2)
and is the primary reason for controlling false positives. Based on definitions of key concepts
in (6) and (7), we state our main theorem on selection consistency of our Bayesian method
even when p is of the order O(n).
Theorem 1 When the observations are generated from (1) for a known η ∈ (1, 2) and σ > 0,
and model (1) is fitted with priors (4) and (5) for the optimal η and with any π0 ∈ (0, 1),
then for p ≤ n, the posterior probability P (γ = γ ∗ | y) → 1 almost surely as n, p → +∞ .
5

The detailed proof of Theorem 1, given in §7, is a non-trivial extension of the proof of
Theorem 1 in Johnson and Rossell (2010b) which used non-local priors to obtain variable
selection consistency. Unlike them, we use a local prior gη (β) to induce a possible nonlocal prior for β given η in (4). To the best of our knowledge, this is the first result on
Bayesian selection consistency when the response distribution is skewed and heteroscedastic
and our method of proof opens the theoretical investigation of sparse Bayesian methods for
transformable models and heterodasticity.

4

Accommodating extremely large errors

Presence of few observations with extremely large errors and their influences on final analysis
for various application areas have been emphasized by many authors including Hampel et al.
(2011). The assumption of Gaussian error density fe in (1) may not be valid due to the
presence of a small number of observations with large errors even after optimal Box-Cox
transformation. To address this, we extend the model (1) to a random location-shift model
with
gη (Yi ) = gη (xTi β) + γi + ei ,
(8)
where γi is nonzero if the ith observation has large error, and zero otherwise. We assume
the vector γ = (γ1 , . . . , γn )T to be sparse to ensure only a small probability of the response
having a large error after transformation. Similar idea of location-shift model, however, with
un-transformed response, is popular in the recent literature on robust linear models (for
example, She and Owen (2011) and McCann and Welsch (2007)). To ensure that gη (xTi β)
is the mean and median of gη (Yi ), we require the mean and median of γi + ei to be zero,
that is, we need a symmetric distribution for γi . For this purpose, we use another spike-andslab mixture prior fγ (γi ) = πγ δ0 + (1 − πγ )φ(γi ; 0, σγ2 ) independently for i = 1, . . . , n, where
0 < πγ < 1.
To induce a heavy-tailed error density after transformation, we also consider another
extension of the model (1) as
−1/2

gη (Yi ) = gη (xTi β) + Ui

with Ui ∼ H(· | ν) ,

ei

(9)

where H(· | ν) is a positive mixing distribution indexed by a parameter ν and ei ’s are again
independent N(0, σ 2 ). This class of heavy-tailed error distributions of (9) is called normal
independent (NI) family (Lange and Sinsheimer, 1993). We consider three kinds of heavy
tailed distribution, Student’s-t, slash and contaminated normal (CN) respectively, for the
marginal error density in (9) using the following specific choices of H(· | ν) (Lachos et al.,
2011): χ2ν /ν distribution with possibly non-integer ν > 2, H(u | ν) = uν for u ∈ [0, 1],
and discrete H(u | ν) with P [ρ < 1] = 1 − P [ρ = 1] = ν. For student-t error, we use the
prior for the degrees of freedom parameter ν to be a truncated exponential on the interval
(2, ∞). For ν of the slash distribution marginal error, we use a Gamma(a, b) prior with
small positive values of a and b with b  a. For contaminated normal marginal error, we
assign Beta(ν0 , ν1 ) and Beta(ρ0 , ρ1 ) priors respectively for ν and ρ. In §5, we compare the
performances of Bayesian analyses under these competing models for highly skewed and
heteroscedastic responses.
6

5

Simulation Studies

Simulation model with no outliers: We use different simulation models to compare
our Bayesian methods under model (1) with LASSO (Tibshirani, 1996) and the penalized
quantile methods (Koenker, 2005). From each simulation model, we simulated 50 replicated
datasets of sample size n = 50. For both simulation studies, the observations are sampled
from the model (1) with ei ∼ N(0, σ02 ) with σ02 = 1. The hyperparameters for priors in (5)
are set as a = 2, b = 2, c1 = 1 and d1 = 1. The tuning parameters for LASSO and penalized
quantile regression are selected via a grid search based on the 5-fold cross-validation. We
compare the estimators from different methods based on following criteria: the mean masking proportion M (fraction of undetected true βj 6= 0), the mean swamping proportion S
(fraction of wrongly selected βj = 0), and the joint detection rate JD (fraction of simulations with 0 masking). We also compare the goodness-of-fit of estimation methods using an
influence measure L/L∗ − 1, where
L=

n
X

T

2

(gη0 (yi ) − gη0 (xi β̂))

/(2σ02 )

−

n/2 log(2πσ02 )

+ (η0 − 1)

i=1

n
X

log(|yi |).

(10)

i=1

is the log-likelihood under (β̂, η0 , σ0 ) and L∗ is the same log-likelihood under (β0 , η0 , σ0 ), and
(β0 , η0 , σ0 ) are the known true parameter values (of the simulation model). The results of
our study using simulated data from TBS model (1) with different values of η are displayed
in Table 1 with p = 8.
Table 1: Results of simulation studies for using different methods of analysis:
Simulation model of (1) with η0 = 0·5, p = 8, β0 = (3, 1·5, 0, 0, 2, 0, 0, 0).
Method used
TBS-SG
Penalized Quantile
LASSO
TBSt-SG
TBSS-SG
TBSCN-SG

L/L∗ − 1
-0·02
0·04
0·04
-0·01
-0·02
-0·02

# of non-zeros
3·16
5·84
4·98
3·16
3·16
3·14

M(%)
0
0
0
0
0
0

S(%)
3·2
56·8
3·96
3·2
3·2
2·8

JD(%)
100
100
100
100
100
100

Simulation model of (1) with η0 = 1·8, p = 8, β0 = (3, 1·5, 0, 0, 2, 0, 0, 0).
Method used
TBS-SG
Penalized Quantile
LASSO
TBSt-SG
TBSS-SG
TBSCN-SG

L/L∗ − 1
-0·06
0·04
0·66
-0·05
-0·06
-0·06

# of non-zeros
3·02
5·82
4·56
3
3
3

M(%)
0
0
0
0
0
0

S(%)
0·4
56·4
3·12
0
0
0

JD(%)
100
100
100
100
100
100

M : masking proportion (fraction of undetected true βj 6= 0); S: swamping proportion S (fraction
of wrongly selected βj with true value 0); JD: joint detection rate.

In Table 1, we compare our Bayesian TBS model (1) with prior (4) for β (called TBSSG in short) to frequentist methods of penalized quantile and LASSO. From the results
in Table 1, it is evident that our TBS-SG method provides better results than competing
7

frequentist methods in terms of average number of non-zeros and swamping error rate. We
also compare TBS-SG method with other Bayesian TBS models with heavy tailed normal
independence (NI) error in (1). These competing NI models in (9) include TBSt-SG model
(in short) with t distribution for H(· | ν), TBSS-SG model (in short) with slash distribution
for H(· | ν) and TBSCN-SG model (in short) with contaminated normal distribution for
H(· | ν). All our Bayesian methods have “SG” in their end of acronym to indicate the spike
Gaussian prior of (4) for β. TBS models accommodating heavy tailed response perform
the best in competing models with ideal masking, swamping and joint outlier detection
rates. Both our methods and frequentist methods provide comparable performances based
on average L/L∗ − 1 values, although the L/L∗ − 1 values from penalized quantile estimates
using different datasets are highly variable. All methods perform desirable with respect
to masking and joint detection. Also, we found that our Bayesian methods provide better
results when true η value is η0 = 1.8, compared to η0 = 0.5.
To compare the performances for different η, we set p = 20 and the number of non-zero
coefficient to be 12. Denote by (x)k the vector formed by appending k copies of x. Consider
case i) β0 = {(2)12 , (0)8 }, case ii) β0 = {(−10)6 , (4)6 , (0)8 }, case iii) β0 = {(−10)10 , (4)2 , (0)8 },
case iv) β0 = {(−10)2 , (−4)2 , (−2)2 , (2)2 , (4)2 , (10)2 , (0)8 }, case v) β0 = {(−10)6 , (2)6 , (0)8 },
case vi) β0 = {(−10)2 , (−8)2 , (−6)2 , (−4)2 , (−2)2 , (2)2 , (0)8 }. We use only TBSCN-SG model
for analysis because these three TBS models accommodating heavy tailed response have
similar performance.
From the results in Table 2, we can clearly see that for all the cases, all the four methods
perform better when η0 = 1.8 compared to when η0 = 0·5, with respect to average number
of non-zeros, masking, swamping and joint detection rate. This can be explained by the fact
that when η0 = 1·8, we expect the posterior draws of η to be close to 1·8 which corresponds
to a non-local prior for β (see Figure 1). When the range of signals is large and when there
are many groups of small coefficients (see case (iv) and case (vi)), all of the methods do not
perform well. Considering only variable selection results (average number of non-zeros), our
TBS model clearly out performs penalized quantile method and LASSO.
Studies using simulation model with outliers and heavy-tailed distribution: Our
simulation models are similar to previous simulation model of (1) except that a few of the
observations are now have large errors even after transformation. Although the Bayesian
TBS methods with NI error in (9) do not provide the identification and estimation of these
observations, we wonder whether they ensure robust variable selection and estimation of β,
particularly in comparison to the Bayesian method using random location-shift model of (8).
For the sake of brevity of the presentation, we omit the tables for results of simulation
studies using data simulated from models (8) and (9), and only summarize the results here.
When we use the simulation model (8) with η0 = 0·5, p = 8, γ(1:2) = 8, γ3 = −8, γ(4:50) = 0,
and β0 = (3, 1·5, 0, 0, 2, 0, 0, 0), our Bayesian method with model (8) obtains 3·32 non-zero
γi ’s on average. The masking (M), swamping (S) and joint detection (JD) rates are 1·33%,
0·77% and 98%. Also, our method provides 3·22 non-zero estimated βj on average with the
masking, swamping and joint detection rates of 0·67%, 4·8% and 98% respectively.
For Simulation 4, we choose η0 = 1·8, p = 8, γ(1:2) = 8, γ3 = −8, γ(4:50) = 0 and
β0 = (3, 1·5, 0, 0, 2, 0, 0, 0). For our Bayesian method with model (8), we have 3·28 nonzero γi ’s on average with the masking, swamping and joint detection rates of 0%, 0·6% and
8

Table 2: Results of simulation studies for using different methods of analysis when the true
model is (1) with p = 20:

Case i)

Case ii)

Case iii)

Case iv)

Case v)

Case vi)

Measurement
L/L∗ − 1
# of non-zeros
M (%)
S(%)
JD(%)
L/L∗ − 1
# of non-zeros
M (%)
S(%)
JD(%)
L/L∗ − 1
# of non-zeros
M (%)
S(%)
JD(%)
L/L∗ − 1
# of non-zeros
M (%)
S(%)
JD(%)
L/L∗ − 1
# of non-zeros
M (%)
S(%)
JD(%)
L/L∗ − 1
# of non-zeros
M (%)
S(%)
JD(%)

TBS-SG
η0 =0·5
-0·05
13·6
0
20
100
-0·03
13·26
0·67
16·75
94
-0·03
12·94
0·33
12·25
96
-0·01
12·58
6·17
16·5
46
0·01
11·32
11·17
8·25
24
-0·01
12·1
6·67
11·25
44

η0 =1·8
2·72
12·02
0
0·25
100
0·03
12
0
0
100
-0·06
12
0
0
100
0·09
12
0
0
100
0·04
12
0
0
100
0·10
12
0
0
100

Penalized Quantile
η0 =0·5 η0 =1·8
0·02
0·7
16·36
15·58
0·5
0
55·25
44·75
94
100
0·05
-0·24
17·12
15·02
0·67
0
65
37·75
94
100
0·03
-0·15
16·94
14·52
0·33
0
62·25
31·5
96
100
0·05
-0·22
16·42
15·32
3·67
0
60·75
41·5
64
100
0·05
-0·21
16·12
14·96
5
0
59
37
62
100
0·06
-0·27
15·62
14·96
4·83
0
52·5
37
62
100

LASSO
η0 =0·5
0·07
14·02
0
25·75
96
0·1
14·96
0
37·5
96
0·11
14·7
0·83
35
90
0·09
14·2
5·67
36
54
0·1
13·56
7·17
30·25
34
0·12
13·52
8·33
31·5
36

η0 =1·8
-82·81
12·24
0
3
100
-1282
12·14
0
1·75
100
-2736
12·24
0
3
100
-689·6
12·08
0·1
2
92
-1072
12·04
0·83
1·75
90
-763·3
12·14
1·33
3·75
84

TBSCN-SG
η0 =0·5 η0 =1·8
-0·04
2·74
13·74
12·02
0
0
21·75
0·25
100
100
-0·02
0·09
13·24
12
8·33
0
16·75
0
92
100
-0·02
0·04
12·92
12
0·5
0
12·25
0
94
100
-0·01
0·13
12·18
12
7·67
0
13·75
0
34
100
0·01
0·09
11·16
12
12
0
7·5
0
14
100
0·01
0·13
11·76
12
8·33
0
9·5
0
36
100

M : masking proportion; S: swamping proportion S; JD: joint detection rate.

9

100%. Also, our method provides 3·22 non-zero estimated βj on average with the masking,
swamping and joint detection rates of 0%, 0% and 100% respectively.
When we simulate data from model (9), the Bayesian method leads to results similar
to the results obtained for using simulation model (8). However, only the Bayesian method
using (8) provides the identification of the observations with errors of large magnitude. In
practice, identification of such observations will facilitate further investigations regarding
their measurement accuracy, influence on inference and other exploratory diagnostics. All of
our Bayesian models provide better results than the penalized quantile regression and LASSO
with respect to average number of non-zeros, masking, swamping and joint detection.

6

Analysis of medical expenditure study

Our motivating study is the Medical Expenditure Panel Survey (Cohen, 2003; Natarajan
et al., 2008), called the MEPS study in short, where the response variable is each patient’s
‘total health care expenditures in the year 2002’. Previous analyses of of this study (Natarajan et al., 2008) suggest that the variance of the response is a function of the mean (heteroscedasticity). Often in practice, medical cost data are typically highly skewed to the right,
because a small percentage of patients may accumulate extremely high costs compared to
other patients, and the variance of total cost tends to increase as the mean increases.
In this article, we focus only on one large cluster because every cluster of MEPS study
has different sampling design. After removing only a few patients with missing observations,
we have 173 patients and 24 potential predictors including age, gender, race, disease history,
etc. The minimum cost is 0 and the maximum is $79660, with a mean $4584 and median
$1342. For the convenience of computation, we standardize the response (cost) and five
potential predictors of the patient: age in 2002, highest education degree attained, perceived
health status, body mass index (BMI), and ability to overcome illness (OVERCOME). Rest
of the potential predictors are binary variables with values 0 and 1. We analyze this study
using our proposed Bayesian models and compare the results with the penalized quantile
regression method of Koenker (2005). For Bayesian methods, we use our transform-bothsides model (1), the model of (8) with sparse large errors (TBSO-SG in short) and the
model of (9) with contaminated normal marginal error. For each method, we compute an
observed residual yi0 − xTi β̂, where yi0 is the observed un-transformed response and xTi β̂ is
the estimated median. The Q-Q plots for the residuals are in Figure 2.
From Q-Q plots in Figure 2, it is obvious that the normality assumption about untransformed response is untenable. We now compare the goodness-of-fit of the three Bayesian
methods to evaluate their abilities from handling skewness and heteroscedasticity. For this
purpose, we use the residual gη̂ (yi0 ) − gη̂ (xTi β̂) for the Bayesian TBS-SG model of (1) and
the TBSCN-SG model of (9), and the residual gη̂ (yi0 ) − gη̂ (xTi β̂) − γ̂i for Bayesian TBSO-SG
model of (8), and then display their Q-Q plots in Figure 3. It is evident from the Q-Q plots
that TBSO-SG model of (8) has the best justification to use it for Bayesian analysis, and
TBSCN-SG model of (9) also performs well except may be for some observations in both
tails.
Using our Bayesian model of (8), we find large posterior evidence of effects of OVERCOME variable with posterior mean=−0·17 and 95% credible interval (−0·21, −0·12), stroke
10

Figure 2: Q-Q plots of observed residuals obtained from 4 methods

11

Figure 3: Q-Q plots for residuals of transformed responses obtained from 3 Bayesian models

12

with posterior mean=0·92 and 95% credible interval (0·66, 1·23), and medication with posterior mean=−0·35 and 95% credible interval (−0·41, −0·24). Model(9) identifies these same
predictors of model (8) with slightly different interval estimates. Model (1) also identifies
three predictors with large posterior evidence of effects: perceived health status, stroke and
the indicator of major ethnic group. Stroke is the only variable with large posterior evidence
of effects in all three models. Even though penalized quantile regression based analysis selects a larger number of predictors compared to the number of predictors selected by our
models, the only statistically significant variable from quantile regression analysis is age (estimate of 0·08 with standard error 0·03). This may be explained by the larger estimated
standard errors of the estimates from quantile regression compared to the posterior standard
deviations of the corresponding parameters obtained via Bayesian analysis.
In order to better understand the prediction performance on observed data, we present
a scatter plot with overlaid quantile lines for each Bayesian method in Figure 4. For each
method, we display scaled xTi β̂ and scaled observed response yi0 , along with estimated 25th
percentile and 75th percentile curves using gη̂−1 {gη̂ (xTi β) + Zα∗ }, where Zα∗ is estimated αpercentile of fe (·). Figure 4 shows that the method using (8) explains the observed data
better than methods using (1) and (9). The observations with large errors identified by
analysis using (8) are marked by asterisk signs in the second plot. We find that all the observations identified by (8) are outside the estimated interquartile ranges. It shows that our
transform-both-sides model of (8) is successful in handling data with skewness, heteroscadesticity as well as very large errors in few subjects. Model of (9) is also able to handle skewness
and heteroscadesticity but is not able to identify observations with extremely large errors.
We also use posterior predictive loss approach (Gelfand and Ghosh, 1998) to evaluate the prediction accuracy underPeach Bayesian method. We compute the prediction
errors of our Bayesian methods by ni=1 E[{gη̂ (yi0 ) − gη̂ (xTi β̂)}2 |D] from model (1), and by
Pn
2
T
i=1 E[{gη̂ (yi0 ) − gη̂ (xi β̂) − γi } |D] for model (8) using MCMC approximation, where D
is the observed dataset. The average prediction error from model (8) is 15 · 03, which is
considerably better than model (1) with average prediction error 174 · 45.

7

Discussion

In this article, we propose Bayesian variable selection methods for skewed and heteroscedastic response. The methods are highly suitable for modeling, computation, analysis and
interpretation of real-life health care cost studies, where we aim to determine and estimate
effects of a sparse set of explanatory variables for health care expenditures out of a large
set of potential explanatory variables. Simulation results indicate a better performance of
our Bayesian methods compared to existing frequentist quantile regression tools. Also, our
Bayesian approaches provide flexible and robust estimations to incorporate a wide variety
of practical situations. The advantages of our Bayesian methods include their practical and
easy implementation using standard statistical software. In the appendix, we prove the consistency of variable selection even when the number of potential predictors p is comparable
to, however, smaller than n. The proofs are only provided for a special case of the covariate
matrix and when the transform parameter η is known. Proof for a more general case can be
obtained following a similar, but more tedious mathematical arguments.
13

Figure 4: Scatter plots of scaled observed responses and quantile regression functions obtained from 3 Bayesian models

14

References
Bach, F. R. (2008), Bolasso: model consistent lasso estimation through the bootstrap, in
‘Proceedings of the 25th International Conference on Machine learning’, ACM, pp. 33–40.
Bhattacharya, A., Pati, D., Pillai, N. and Dunson, D. (2014), ‘Dirichlet-laplace priors for
optimal shrinkage’, JASA .
Bickel, P. J. and Doksum, K. A. (1981), ‘An analysis of transformations revisited’, JASA
76(374), 296–311.
Bondell, H. D. and Reich, B. J. (2012), ‘Consistent high-dimensional bayesian variable selection via penalized credible regions’, JASA 107(500), 1610–1624.
Carvalho, C., Polson, N. and Scott, J. (2010), ‘The horseshoe estimator for sparse signals’,
Biometrika 97(2), 465–480.
Castillo, I., Schmidt-Hieber, J. and van der Vaart, A. W. (2014), ‘Bayesian linear regression
with sparse priors’, ArXiv e-prints .
Chen, L., Pourahmadi, M. and Maadooliat, M. (2014), ‘Regularized multivariate regression models with skew-t error distributions’, Journal of Statistical Planning and Inference
149, 125 – 139.
Cohen, S. B. (2003), ‘Design strategies and innovations in the medical expenditure panel
survey’, Medical Care 41(7), III.
Efron, B., Hastie, T., Johnstone, I. and Tibshirani, R. (2004), ‘Least angle regression’, The
Annals of statistics 32(2), 407–499.
Fan, J. and Li, R. (2001), ‘Variable selection via nonconcave penalized likelihood and its
oracle properties’, JASA 96(456), 1348–1360.
Gelfand, A. E. and Ghosh, S. K. (1998), ‘Model choice: A minimum posterior predictive loss
approach’, Biometrika 85(1), 1–11.
George, E. I. and McCulloch, R. E. (1997), ‘Approaches for Bayesian variable selection’,
Statistica Sinica 7(2), 339–373.
George, E. and McCulloch, R. (1993), ‘Variable selection via Gibbs sampling’, JASA
88(423), 881–889.
Hampel, F. R., Ronchetti, E. M., Rousseeuw, P. J. and Stahel, W. A. (2011), Robust Statistics: the Approach Based on Influence Functions, Vol. 114, Wiley.
Hogg, R., McKean, J. and Craig, A. (2013), Introduction to Mathematical Statistics, Pearson.
Johnson, V. E. and Rossell, D. (2010a), ‘On the use of non-local prior densities in bayesian
hypothesis tests’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72(2), 143–170.
15

Johnson, V. and Rossell, D. (2010b), ‘On the use of non-local prior densities in Bayesian hypothesis tests’, Journal of the Royal Statistical Society: Series B (Statistical Methodology)
72(2), 143–170.
Koenker, R. (2005), Quantile Regression, Vol. 38, Cambridge University Press.
Lachos, V. H., Bandyopadhyay, D. and Dey, D. K. (2011), ‘Linear and nonlinear mixed-effects
models for censored hiv viral loads using normal/independent distributions’, Biometrics
67(4), 1594–1604.
Lange, K. and Sinsheimer, J. S. (1993), ‘Normal/independent distributions and their applications in robust regression’, Journal of Computational and Graphical Statistics 2(2), 175–
198.
Lin, J., Sinha, D., Lipsitz, S. and Polpo, A. (2012), ‘Semiparametric bayesian survival analysis using models with log-linear median’, Biometrics 68(4), 1136–1145.
McCann, L. and Welsch, R. E. (2007), ‘Robust variable selection using least angle regression
and elemental set sampling’, Computational Statistics & Data Analysis 52(1), 249–257.
Mitchell, T. J. and Beauchamp, J. J. (1988), ‘Bayesian variable selection in linear regression’,
JASA 83(404), 1023–1032.
Natarajan, S., Lipsitz, S. R., Fitzmaurice, G., Moore, C. G. and Gonin, R. (2008), ‘Variance
estimation in complex survey sampling for generalized linear models’, Journal of the Royal
Statistical Society: Series C (Applied Statistics) 57(1), 75–87.
Park, T. and Casella, G. (2008), ‘The Bayesian lasso’, JASA 103(482), 681–686.
Rossell, D. and Telesca, D. (2015), ‘Non-local priors for high-dimensional estimation’, Journal of the American Statistical Association 0(just-accepted), 1–33.
She, Y. and Owen, A. B. (2011), ‘Outlier detection using nonconvex penalized regression’,
JASA 106(494), 626–639.
Tibshirani, R. (1996), ‘Regression shrinkage and selection via the lasso’, Journal of the Royal
Statistical Society. Series B (Methodological) pp. 267–288.
Walker, A. M. (1969), ‘On the asymptotic behaviour of posterior distributions’, Journal of
the Royal Statistical Society. Series B (Methodological) 31(1), 80–88.
Zhang, C.-H. (2010), ‘Nearly unbiased variable selection under minimax concave penalty’,
The Annals of Statistics pp. 894–942.
Zhao, P. and Yu, B. (2006), ‘On model selection consistency of lasso’, J. Mach. Learn. Res.
7, 2541–2563.
Zou, H. (2006), ‘The adaptive lasso and its oracle properties’, JASA 101(476), 1418–1429.
Zou, H. and Hastie, T. (2005), ‘Regularization and variable selection via the elastic net’,
Journal of the Royal Statistical Society, Series B 67, 301–320.
16

Proof of Theorem 1
For γ ∈ {0, 1}p and x ∈ Rp , define xγ ∈ Rpγ to be the vector (xj : γj 6= 0). The log-likelihood
corresponding to (1) has the expression
`(β) =

n
X
i=1

 n
1
log gη0 (yi ) − log(2π) − n log(σ) − 2 kgη (y) − gη (Xβ)k2 ,
2
2σ

(11)

T

T
where gη (y) = gη (y1 ), · · · , gη (yn ) , gη (Xβ) = gη (xT1 β), · · · , gη (xTn β) . The gradient of
(11) is given by

1
∇`(β) = 2 X T D̄η gη (y) − gη (Xβ) ,
(12)
σ
η−1
and D̄η = diag(d¯η ). β 7→ ∇`(β)
where d¯η = (gη0 (xT1 β), · · · , gη0 (xTn β)) with gη0 (xTi β) = xTi β
is continuous. Since the number of true signals βj∗ is finite, we assume βj∗ ∈ [ML , MU ]. Hence,
η−1
as long as β is in a neighborhood of β ∗ , ∇`(β) is bounded since xTi β
in (12) is bounded
when η ∈ (1, 2). The Hessian of (11) is defined as
H(β) =

1 T
T
¯ T M (y, X)X,
X
−
D̄
D̄
+
D̄
η
β
η
η
σ2

(13)

where d¯η is defined as the element-wise second derivative of gη (xTi β) on β with gη00 (xTi β) =
η−2
¯ = diag(d¯ ). Meanwhile, m (y, X) = g (y )−g (xT β), · · · , g (y )−
and D̄
xTi β
sgn(xTi β)(η−1)
η n
η
η
η
η 1
η
1

T
gη (xn β) and Mη (y, X) = diag(mη (y, X)). Using Laplace approximation, the Bayes factor
can be approximated as

pγ −pγ ∗
|H(βγ∗∗ |γ ∗ )|1/2
π(β̃γ |γ)
p(y|β, γ)
{`(β̃γ )−`(βγ∗∗ )}
2
=
(2π)
×
e
×
×
.
(14)
p(y|β ∗ , γ ∗ )
π(βγ∗∗ |γ ∗ )
|H(β̃γ |γ)|1/2
n
The first term in the r.h.s of (14) is Op (1) with p ≤ n. For the second term, we denote
by λ(y) = `(β̃) − `(β ∗ ) the likelihood ratio statistic. As in Rossell and Telesca (2015,
Proposition 3), it is straightforward to verify that our sampling model (1) satisfies Walker’s
conditions (A1)-(A5) and (B1)-(B4) (Walker, 1969). Hence, our MLE β̃i is consistent and
the Hessian matrix H(β̃γ |γ) in (14) converges in probability. We consider two cases below.
∗
P
→
When γ ∗ 6⊂ γ, i.e., γ misses some true active coefficients, the second term e(logL(β̃γ )−logL(βγ ∗ )) −
∗
∗
e−nKL(p(y|β ,γ ),p(y|β,γ)) where KL(p(y|β ∗ , γ ∗ ), p(y|β, γ)) is the Kullback-Leibler divergence between optimal p(y|β ∗ , γ ∗ ) and p(y|β, γ) under γ. Here the minimum KL divergence KL(p(y|β ∗ , γ ∗ ), p(y|β, γ))
is strictly positive since,
n Z
X
∗
∗
[log(p(yi |β ∗ , γ ∗ )) − log(p(yi |β, γ))]p(yi |β ∗ , γ ∗ )dyi
KL(p(y|β , γ ), p(y|β, γ)) =

i=1
n
2

1 X
T
T ∗ 2
=
E gη (yi ) − gη (xi β) − E gη (yi ) − gη (xi β )
2σ 2 i=1

=

n
2
1 X
gη (xTi β ∗ ) − gη (xTi β) > 0
2
2σ i=1

17

(15)

when β 6= β ∗ satisfies the eigenvalue conditions in §3 that indicates no linear dependency
among covariates xi . Therefore the second term is Op (e−n ) when γ ∗ 6⊂ γ.
When γ ∗ ⊂ γ, we denote the likelihood-ratio statistic by Λ(y) = logL(β̃γ ) − logL(βγ∗∗ ).
Under appropriate regularity conditions (Hogg et al., 2013), our likelihood ratio statistic Λ(y)
is asymptotically chi-square distributed. The regularity conditions relevant to the argument
are listed as (R0)-(R9) in Hogg et al. (2013) where (R0)-(R2) and (R6)-(R8) can be obtained
trivially. Conditions (R3) and (R4) are related to Fisher information and are satisfied by
Hessian matrix (13). Conditions (R5) and (R9) essentially guarantee that the remainder of
a second order Taylor expansion around β is bounded in probability. To that end, note that
∂3
`(β)
∂βj ∂βk ∂βl
n
h
 T η−3 i
 T 2η−3
η−1X
T
T
+ (η − 2) gη (yi ) − gη (xi β) |xi β|
=
xij xik xil − 2 + sgn(xi β) |xi β|
σ 2 i=1
n
h
i
η−1X
T
2η−3
T
T
η−3
≤
3|x
β|
+
|
g
(y
)
−
g
(x
β)
|
|x
β|
:= Mjkl (y; X)
(16)
x
x
x
η i
η
ij ik il
i
i
i
σ 2 i=1
where E[Mjkl (y; x)] < ∞ for all j, k, l ∈ 1, · · · , p. Therefore our model (1) satisfies all
regularity conditions implying Λ(y) ∼ χ2pγ −pγ ∗ and hence Op (1) as required. When γ has
moderate size with pγ ≤ n, we will show the first term in (14) is dominated by the second
term later.
Next consider the second term π(β̃|γ)/π(β ∗ |γ ∗ ) under a non-local prior with η ∈ (1, 2)
is known.
Qpγ
2
2
η−1
√ 1
i=1 2πσβ exp(−β̃γi /2σβ )|β̃γi |
π(β̃|γ)
= Qpγ ∗ 1
∗2
2
∗ η−1
√
π(β ∗ |γ ∗ )
j=1 2πσβ exp(−βγj∗ /2σβ )|βγj∗ |
Q
X
X
 |β̃γi |η−1
1
2
∗2
2
pγ −pγ ∗
= (√
β̃γi −
βγj∗ )/2σβ Q ∗ η−1
(17)
)
exp − (
|βγj∗ |
2πσβ
First if γ ∗ ⊂ γ, given that β̃γi = Op (n−1/2 ) and η − 1 > 0, the second term is ensentially Op (n−(η−1)(pγ −pγ ∗ )/2 ). When γ ∗ 6⊂ γ, denote s∗ = γ ∩ γ ∗ , the second term is
Op (n−(η−1)(pγ −ps∗ )/2 ) and upper bounded by Op (n−(η−1)/2 ).
To conclude the proof, we need to deal with the the third term in (14). The Hessian
matrix H(β̃) is given by (13) with each element Hij = XiT Xj Op (1). Since β˜γ converge in
p
probability to βγ ∗ , we have n−1 H(β̃|γ) →
− H(β̃|γ). Therefore by continous mapping theorem,
the third term in (14) is approximated as
|H(β ∗ |γ ∗ )|1/2
1/2

H(β̃|γ)



npγ ∗ /2 |n−1 H(β0 |γ ∗ )|
= Op (n(pγ ∗ −pγ )/2 ).
npγ /2 n−1 H(β̃|γ)

(18)

To conclude, the Bayes factor (14) is Op (n−η(pγ −pγ ∗ )/2 ) when γ ∗ ⊂ γ, and Op (n(pγ ∗ −pγ )/2 e−n )

18

when γ ∗ 6⊂ γ. Then the posterior probability p(γ = γ ∗ | y) can be lower-bounded as
p(y | γ ∗ )π(γ ∗ )
P
p(y | γ ∗ )π(γ ∗ ) + γ6=γ ∗ p(y | γ)π(γ)
h
X p(y | γ)π(γ) i−1
X p(y | γ)π(γ)
+
= 1+
p(y | γ ∗ )π(γ ∗ ) γ ∗ 6⊂γ p(y | γ ∗ )π(γ ∗ )
γ ∗ ⊂γ
h
i−1
X
X
≥ 1+
Op (n−η(pγ −pγ ∗ )/2 ) +
Op (n(pγ ∗ −pγ )/2 e−n )

p(γ = γ ∗ | y) =

γ ∗ ⊂γ

γ ∗ 6⊂γ

h
i−1
X
X
a.s.
−(pγ −pγ ∗ )/2
−n
≥ 1+
Op (n
)+
Op (e )
−−→ 1
γ ∗ ⊂γ

γ ∗ 6⊂γ

which concludes the proof.

19

