# Robust Bayesian Regularized Estimation Based ont Regression Model

Authors: Zean Li; Weihua Zhao
Journal: Journal of Probability and Statistics
 (2015)

---

## Abstract

Thet distribution is a useful extension of the normal distribution, which can be used for statistical modeling of data sets with heavy tails, and provides robust estimation. In this paper, in view of the advantages of Bayesian analysis, we propose a new robust coefficient estimation and variable selection method based on Bayesian adaptive Lasso t regression. A Gibbs sampler is developed based on the Bayesian hierarchical model framework, where we treat the t distribution as a mixture of normal and gamma distributions and put different penalization parameters for different regression coefficients. We also consider the Bayesian t regression with adaptive group Lasso and obtain the Gibbs sampler from the posterior distributions. Both simulation studies and real data example show that our method performs well compared with other existing methods when the error distribution has heavy tails and/or outliers.


---

## Body

## 1. Introduction

Following the pioneer work of Tibshirani [1], the Lasso (least absolute shrinkage and selection operator) has generated much interest in statistical literatures [2–5]. For the Gaussian linear regression,(1)yi=xiTβ+ϵi,i=1,…,n,where xi’s are either fixed or random covariates, ϵ1,…,ϵn, i.i.d.~N(0,σ2). The Lasso estimates are viewed as L1 regularized least squares estimates, which are defined as(2)minβ⁡y-xβTy-xβ+λ∑j=1pβj,for λ≥0 and with n×1 vectors y=(y1,…,yn)T, n×p matrix x=(x1,…,xn)T, and p×1 vector β. The key advantage of Lasso lies in its ability to do simultaneous parameter estimation and variable selection and the entire path of Lasso estimates for all values of λ can be efficiently computed through a modification of the LARS (least angle regression) algorithm of Efron et al. [5].The Gaussian assumption is not crucial in model (1), but it is useful to make connections to the likelihood framework for regularized t regression in Section 2. The Lasso estimator in (2) is equivalent to minimizing the penalized negative log-likelihood l(β;y) as a function of the regression coefficients β and using the L1-penalty β=∑j=1p|βj|: equivalence means here that we obtain the same estimator for a potentially different tuning parameter, but the Lasso estimator in (2) does not provide an estimate of the variance parameter σ2. In fact, the variance parameter σ2 plays an important role especially when the error in the regression model has high variance.On the other hand, the Lasso estimate in (2) based on least square loss may suffer from poor performance when the error distribution has heavy tail which is longer than normal distribution or the error has large variance or contains outliers in the linear regression model, which motivate many researches to consider regularized estimation based on robust method. The most existing robust regularized estimation methods mainly replace the least square loss function in (2) by some robust loss functions, such as Huber loss, L1 loss, and quantile loss function. Li and Zhu [6] considered quantile regression with the Lasso penalty and developed its piecewise linear solution path. Wang et al. [7] investigated the least absolute deviance (LAD) estimate with adaptive Lasso penalty (LAD-Lasso) and proved its oracle property. Wu and Liu [8] further discussed the oracle properties of the SCAD (smoothly clipped absolute deviation) and adaptive Lasso regularized quantile regression. Zou and Yuan [9] proposed a new efficient estimation and variable selection procedure which is called composite quantile regression. Chen et al. [10] studied the robust Lasso in the presence of noise with large variance. They considered two editions of robust Lasso, which use the convex combined loss and Huber loss function criterion instead of least square criterion.Besides the above existing robust regularized method, which, in essence, replace theL2 loss by other loss criteria, we can also directly consider robust estimation from the standpoint of error distribution in model (1). It is well known that when the error in the linear regression follows a normal distribution, the estimations will be sensitive to outliers, which motivate us to consider using robust error distribution in the linear regression model in (1). The t distribution as an alternative to the normal distribution has frequently been suggested in the literature as a robust extension of the traditional normal models. The t distribution provides a convenient description for regression analysis when the residual term has density with heavy tails and excess kurtosis.Since the Lasso estimate in (2) can be interpreted as a Bayesian posterior mode estimate when the regression parameters have independent Laplace priors, motivated by this connection, Park and Casella [11] discussed the Bayesian Lasso estimation for linear regression based on Gibbs sampler method. The Bayesian Lasso has several advantages; for example, it can conveniently provide interval estimates, that is, Bayesian credible intervals, that can guide variable selection, the structure of the hierarchical model provides methods for selecting the Lasso parameter, and sometimes its performances outperform non-Bayesian method. Therefore, several authors were subsequently devoted to the variable selection problem based on Bayesian method, such as Park and Casella [11], Hans [12], Li et al. [13], and Kyung et al. [14].The main goal of this paper is to extend the Bayesian regularized method fort regression model (4), which is a robust edition Bayesian regularized estimation method, and it is expected that the performances will be better than other regularized methods when the error distribution follows a symmetry distribution and there are outliers in regression model.The rest of the paper is organized as follows. In Section2, we introduce the t regression with Lasso penalty. In Section 3, we build the framework of Bayesian t regression model with adaptive Lasso penalty and then show the corresponding Bayesian hierarchical model and obtain the Gibbs sample estimation procedure. In Sections 4 and 5, we carry out some simulation studies and real data analysis to illustrate the performance of our new method. Some conclusion and remarks are contained in Section 6. We also discuss the Bayesian t regression with group Lasso penalty and its Gibbs sampler method in the Appendices.


## 2.t Regression Model and L1 Norm Regularization

### 2.1.t Regression Model

Thet distribution as an alternative to the normal distribution has frequently been suggested in the literatures; for example, Lange et al. [15] introduced t distributed errors regression models as a robust extension of the traditional normal models; Liu and Rubin [16] discussed the estimations of t distribution by EM algorithm and its extensions; and Lin et al. [17] studied heteroscedasticity diagnostics problem by score test for t linear regression models. And it has often been used in many real applications, especially in finance.The probability density function of at distributed variable y with location parameter μ, dispersion parameter σ, and degree of freedom v, denoted as t(μ,σ2,v), is(3)fy,μ,σ;v=Γ1+v/2σπv1/2Γv/21+v-1y-μσ2-1+v/2, where Γ(·) denotes the gamma function, y,μ∈R, and σ>0.The mean and variance ofy are μ(v>1) and vσ2/(v-2)(v>2), respectively. When v=∞ and 1, the t distribution reduces to the normal and Cauchy distributions, respectively. The univariate linear Student-t regression model can be expressed as(4)yi=xiTβ+ɛi,ɛi~t0,σ2,v,i.i.d.,where the xi’s are known p-vector of covariates and β=(β1,…,βp)T are the unknown p-vectors of parameters.


### 2.2.t Regression with Lasso Penalty

As discussed in Section1, the L2 loss in the Lasso model is not robust to heavy-tailed error distribution and/or outliers. This indicates that the Lasso is not an ideal goodness-of-fit measure criterion in the presence of noise with large variance. We consider following regularized t regression, which is defined as(5)θ^=argmaxβ∈Rp,σ2>0⁡∑i=1nliβ,σ2;yi,xi-λ∑j=1pβj,where θ=(βT,ϕ)T and li(β,σ2;yi,xi) is the log-likelihood function of yi:(6)liβ,σ2;yi,xi=log⁡Γ1+v/2πv1/2Γv/2-12log⁡σ2-1+v2log⁡1+v-1yi-μi2σ2.Comparing (5) with (2), we take the parameter σ2 into the optimization of the penalized maximum likelihood estimator. Though we are penalizing only the parameter β, the variance parameter estimate σ^2 is influenced indirectly by the amount of shrinkage λ in (5).There is main drawback of the estimator in (5) [18]. It is not equivalent under scaling of the response [19]. More precisely, consider the transformation(7)y¯=by,β¯=bβ,σ¯=bσb>0,which leaves model (2) invariant, but this is not the case for the estimator in (5). We address this drawback by using the penalty term λβ1/σ, leading to the following estimator:(8)θ^=argmaxβ∈Rp,σ2>0⁡∑i=1nliβ,σ2;yi,xi-λ∑j=1pβjσ.Noting the form of the regularized estimation in (8), the Lasso estimates can be interpreted as posterior mode estimates when the regression parameters have independent Laplace (i.e., double-exponential) priors:(9)fβj∣σ2,λ=λ2σ2exp⁡-λβjσ2.Park and Casella [11] argued that conditioning on σ2 is important because it guarantees a unimodal full posterior.


## 2.1.t Regression Model

Thet distribution as an alternative to the normal distribution has frequently been suggested in the literatures; for example, Lange et al. [15] introduced t distributed errors regression models as a robust extension of the traditional normal models; Liu and Rubin [16] discussed the estimations of t distribution by EM algorithm and its extensions; and Lin et al. [17] studied heteroscedasticity diagnostics problem by score test for t linear regression models. And it has often been used in many real applications, especially in finance.The probability density function of at distributed variable y with location parameter μ, dispersion parameter σ, and degree of freedom v, denoted as t(μ,σ2,v), is(3)fy,μ,σ;v=Γ1+v/2σπv1/2Γv/21+v-1y-μσ2-1+v/2, where Γ(·) denotes the gamma function, y,μ∈R, and σ>0.The mean and variance ofy are μ(v>1) and vσ2/(v-2)(v>2), respectively. When v=∞ and 1, the t distribution reduces to the normal and Cauchy distributions, respectively. The univariate linear Student-t regression model can be expressed as(4)yi=xiTβ+ɛi,ɛi~t0,σ2,v,i.i.d.,where the xi’s are known p-vector of covariates and β=(β1,…,βp)T are the unknown p-vectors of parameters.


## 2.2.t Regression with Lasso Penalty

As discussed in Section1, the L2 loss in the Lasso model is not robust to heavy-tailed error distribution and/or outliers. This indicates that the Lasso is not an ideal goodness-of-fit measure criterion in the presence of noise with large variance. We consider following regularized t regression, which is defined as(5)θ^=argmaxβ∈Rp,σ2>0⁡∑i=1nliβ,σ2;yi,xi-λ∑j=1pβj,where θ=(βT,ϕ)T and li(β,σ2;yi,xi) is the log-likelihood function of yi:(6)liβ,σ2;yi,xi=log⁡Γ1+v/2πv1/2Γv/2-12log⁡σ2-1+v2log⁡1+v-1yi-μi2σ2.Comparing (5) with (2), we take the parameter σ2 into the optimization of the penalized maximum likelihood estimator. Though we are penalizing only the parameter β, the variance parameter estimate σ^2 is influenced indirectly by the amount of shrinkage λ in (5).There is main drawback of the estimator in (5) [18]. It is not equivalent under scaling of the response [19]. More precisely, consider the transformation(7)y¯=by,β¯=bβ,σ¯=bσb>0,which leaves model (2) invariant, but this is not the case for the estimator in (5). We address this drawback by using the penalty term λβ1/σ, leading to the following estimator:(8)θ^=argmaxβ∈Rp,σ2>0⁡∑i=1nliβ,σ2;yi,xi-λ∑j=1pβjσ.Noting the form of the regularized estimation in (8), the Lasso estimates can be interpreted as posterior mode estimates when the regression parameters have independent Laplace (i.e., double-exponential) priors:(9)fβj∣σ2,λ=λ2σ2exp⁡-λβjσ2.Park and Casella [11] argued that conditioning on σ2 is important because it guarantees a unimodal full posterior.


## 3. Bayesian Formulation of Robustt Adaptive Lasso

### 3.1. Bayesian Adaptive Lassot Regression

It is well known that the ordinary Lasso has a conflict between the consistency for model selection and optimal prediction. As a solution to achieve both estimation and model selection consistency, the adaptive Lasso [20] was proposed and its model selection consistency and asymptotic normality under certain rate of shrinkage parameter λn were proved. To extend the idea of the adaptive Lasso [20] to our proposed robust t Lasso regression, we define the robust t adaptive Lasso as(10)θ^=argmaxβ∈Rp,ϕ>0⁡∑i=1nliβ,ϕ;yi,xi-∑j=1pϕ1/2λjβj,where λj’s allow unequal penalties for the coefficients and we define the precision parameter ϕ=σ-2.Motivated by the connection with Bayesian estimation and (9), we consider a fully Bayesian analysis using a conditional Laplace prior on βj:(11)fβj∣ϕ,λj=ϕλj2exp⁡-λjϕβj.For any a≤0, by the following equality [21]:(12)a2e-at=∫0∞12πsexp⁡-t22sa22exp⁡-a2s2ds,the Laplace prior (11) on β can be written as(13)fβ∣ϕ,λ1,…,λp=∏j=1p∫0∞12πsjexp⁡-βj22sjϕλj22exp⁡-ϕλj2sj2dsj. On the other hand, from Lange et al. [15] and Liu and Rubin [16], t(μ,ϕ,v) can be treated as a mixture of normal and gamma distributions; that is, if u~Γ(v/2,v/2) and z∣u~N(μ,(ϕu)-1), then z~t(μ,ϕ,v). If we further put gamma priors on the parameter λj2 and use the improper prior density on ϕ, then we have the following Bayesian hierarchical model:(14)yi=xiTβ+ϕui-1/2zi,fϕ∝1ϕ,fu∝∏i=1nuiv/2-1exp⁡-v2ui,fz∝∏i=1n12πexp⁡-12zi2,fβ,s∣λ12,…,λp2∝∏j=1p12πsjexp⁡-βj22sjϕλj22exp⁡-ϕλj2sj2,fλ12,…,λp2∣a,b∝∏j=1pbaΓaλj2a-1exp⁡-bλj2,fa,b∝1b,where u=(u1,…,un), z=(z1,…,zn), and s=(s1,…,sp).Remark 1.
The amount of shrinkage parametersλj2’s depend on the value of hyperparameters a and b. Because larger a and smaller b lead to bigger penalization, it is important to treat a and b as unknown parameters to avoid giving specific values which affect the estimates of the regression coefficients.Remark 2.
Similar to Park and Casella [11], there is also another existing method to obtain the estimation of adaptive Lasso parameters beside sampling from the posterior distribution of λj2’s. We can use the empirical Bayes estimation method to update λj2’s without setting the prior of λj2’s in the Bayesian hierarchical model. Our little experience shows that the two methods perform well, but the empirical Bayes estimation method can save some computing time; for more details see Park and Casella [11] and Leng et al. [22].


### 3.2. Gibbs Sampler for All Parameters

Letβ-i be the parameter vector β excluding the component βi, s-j the vector s excluding the component sj, and u-i the variable excluding the component ui. Based on the Bayesian hierarchical model (14), the posterior distribution of all parameters can be given by(15)fβ,ϕ,u,s,λ12,…,λp2,a,b∣y,x∝fy∣β,ϕ,u,xfϕ∏i=1nfui∏j=1pfβj,sj∣λj2fλj2∣a,bfa,b∝ϕn/2-1exp⁡-ϕ2∑i=1nuiyi-xiTβ2∏i=1nuiv-1/2exp⁡-v2ui∏j=1p12πsjexp⁡-βj22sjϕλj22exp⁡-ϕλj2sj2baΓaλj2a-1exp⁡-bλj2×b-1.Based on the expression (15), we can obtain the posterior full conditional distribution for all parameters, which can yield a tractable and efficient Gibbs sampler that works as follows.The full conditional distributionfui∣x,y,β,ϕ,u-i,s,λ12,…,λp2,a,b is(16)fui∣x,y,β,ϕ,u-i,s,λ12,…,λp2,a,b∝uiv+1/2-1exp⁡-ϕ2yi-xiTβ2+v2ui. Thus, the full conditional distribution of ui is a gamma distribution.The full conditional distributionfsj∣x,y,β,ϕ,s-j,u,λ12,…,λp2,a,b is(17)fsj∣x,y,β,ϕ,s-j,u,λ12,…,λp2,a,b∝1sjexp⁡-12ϕλj2sj+βj2sj-1. Thus, the full conditional distribution of sj is a generalized inverse Gaussian distribution. Now we consider the full conditional distribution of βj which is given by(18)fβj∣x,y,β-j,ϕ,s,u,λ12,…,λp2,a,b∝exp⁡-12ϕ∑i=1nuixij2+1sjβj2-2ϕ∑i=1nuixijy~ijβj, where y~ij=yi-∑k=1,k≠jxikβk. Let σ~j-2=ϕ∑i=1nuixij2+1/sj and μ~j=σ~2ϕ∑i=1nuixijy~ij. Then the full conditional distribution of βj is just the normal distribution N(μ~j,σ~j2). The full conditional distribution of ϕ is(19)fϕ∣x,y,β,s,u,λ12,…,λp2,a,b∝ϕn/2+p-1exp⁡-ϕ∑i=1nuiyi-xiTβ2+∑j=1pλj2sj2. So the full conditional distribution of ϕ is a gamma distribution.The full conditional distribution ofλj2 is(20)fλj2∣x,y,β,ϕ,s,u,λk,k=1,…,p,k≠j2,a,b∝λj2aexp⁡-ϕsj2+bλj2. That is, the full conditional distribution of λj2 is a gamma distribution. At last, the full conditional distributions of a and b are(21)fa∣x,y,β,ϕ,s,u,λ12,…,λp2,b∝baΓap∏j=1pλj2a,fb∣x,y,β,ϕ,s,u,λ12,…,λp2,a∝bap-1exp⁡-b∑j=1pλj2. Obviously, the full conditional distribution of b is just a gamma distribution. However, the full conditional posterior distribution of c does not have a closed form. Fortunately, we note that the full conditional posterior distribution of c is a log-concave function, so the adaptive rejection sampling algorithm [23] can be used to sample from this distribution. The adaptive rejection sampling algorithm is an efficient sampling method from any univariate probability density function which is log-concave.Remark 3.
The Lasso was originally designed for variable selection for the posterior mode ofβj’s can be exactly zeros due to the nature of the laplace prior in (11). However, the posterior draws of βj’s cannot be exact zeros since βj’s are drawn from the continuous posterior distribution. A post hoc thresholding rule may overcome this difficulty. Alternatively, kyung et al. [24] recommended using the credible interval on the posterior mean to guide variable selection.Remark 4.
For variable selection, sometimes, several explanatory variables may be represented by a group of derived input variables. In this case the selection of important variables corresponds to the group of variables. The group Lasso penalty [13, 25] takes the group structure into account and can do variable selection at the group level. For the Bayesian t regression model with adaptive group Lasso penalty, we give the details in the Appendices.


## 3.1. Bayesian Adaptive Lassot Regression

It is well known that the ordinary Lasso has a conflict between the consistency for model selection and optimal prediction. As a solution to achieve both estimation and model selection consistency, the adaptive Lasso [20] was proposed and its model selection consistency and asymptotic normality under certain rate of shrinkage parameter λn were proved. To extend the idea of the adaptive Lasso [20] to our proposed robust t Lasso regression, we define the robust t adaptive Lasso as(10)θ^=argmaxβ∈Rp,ϕ>0⁡∑i=1nliβ,ϕ;yi,xi-∑j=1pϕ1/2λjβj,where λj’s allow unequal penalties for the coefficients and we define the precision parameter ϕ=σ-2.Motivated by the connection with Bayesian estimation and (9), we consider a fully Bayesian analysis using a conditional Laplace prior on βj:(11)fβj∣ϕ,λj=ϕλj2exp⁡-λjϕβj.For any a≤0, by the following equality [21]:(12)a2e-at=∫0∞12πsexp⁡-t22sa22exp⁡-a2s2ds,the Laplace prior (11) on β can be written as(13)fβ∣ϕ,λ1,…,λp=∏j=1p∫0∞12πsjexp⁡-βj22sjϕλj22exp⁡-ϕλj2sj2dsj. On the other hand, from Lange et al. [15] and Liu and Rubin [16], t(μ,ϕ,v) can be treated as a mixture of normal and gamma distributions; that is, if u~Γ(v/2,v/2) and z∣u~N(μ,(ϕu)-1), then z~t(μ,ϕ,v). If we further put gamma priors on the parameter λj2 and use the improper prior density on ϕ, then we have the following Bayesian hierarchical model:(14)yi=xiTβ+ϕui-1/2zi,fϕ∝1ϕ,fu∝∏i=1nuiv/2-1exp⁡-v2ui,fz∝∏i=1n12πexp⁡-12zi2,fβ,s∣λ12,…,λp2∝∏j=1p12πsjexp⁡-βj22sjϕλj22exp⁡-ϕλj2sj2,fλ12,…,λp2∣a,b∝∏j=1pbaΓaλj2a-1exp⁡-bλj2,fa,b∝1b,where u=(u1,…,un), z=(z1,…,zn), and s=(s1,…,sp).Remark 1.
The amount of shrinkage parametersλj2’s depend on the value of hyperparameters a and b. Because larger a and smaller b lead to bigger penalization, it is important to treat a and b as unknown parameters to avoid giving specific values which affect the estimates of the regression coefficients.Remark 2.
Similar to Park and Casella [11], there is also another existing method to obtain the estimation of adaptive Lasso parameters beside sampling from the posterior distribution of λj2’s. We can use the empirical Bayes estimation method to update λj2’s without setting the prior of λj2’s in the Bayesian hierarchical model. Our little experience shows that the two methods perform well, but the empirical Bayes estimation method can save some computing time; for more details see Park and Casella [11] and Leng et al. [22].


## 3.2. Gibbs Sampler for All Parameters

Letβ-i be the parameter vector β excluding the component βi, s-j the vector s excluding the component sj, and u-i the variable excluding the component ui. Based on the Bayesian hierarchical model (14), the posterior distribution of all parameters can be given by(15)fβ,ϕ,u,s,λ12,…,λp2,a,b∣y,x∝fy∣β,ϕ,u,xfϕ∏i=1nfui∏j=1pfβj,sj∣λj2fλj2∣a,bfa,b∝ϕn/2-1exp⁡-ϕ2∑i=1nuiyi-xiTβ2∏i=1nuiv-1/2exp⁡-v2ui∏j=1p12πsjexp⁡-βj22sjϕλj22exp⁡-ϕλj2sj2baΓaλj2a-1exp⁡-bλj2×b-1.Based on the expression (15), we can obtain the posterior full conditional distribution for all parameters, which can yield a tractable and efficient Gibbs sampler that works as follows.The full conditional distributionfui∣x,y,β,ϕ,u-i,s,λ12,…,λp2,a,b is(16)fui∣x,y,β,ϕ,u-i,s,λ12,…,λp2,a,b∝uiv+1/2-1exp⁡-ϕ2yi-xiTβ2+v2ui. Thus, the full conditional distribution of ui is a gamma distribution.The full conditional distributionfsj∣x,y,β,ϕ,s-j,u,λ12,…,λp2,a,b is(17)fsj∣x,y,β,ϕ,s-j,u,λ12,…,λp2,a,b∝1sjexp⁡-12ϕλj2sj+βj2sj-1. Thus, the full conditional distribution of sj is a generalized inverse Gaussian distribution. Now we consider the full conditional distribution of βj which is given by(18)fβj∣x,y,β-j,ϕ,s,u,λ12,…,λp2,a,b∝exp⁡-12ϕ∑i=1nuixij2+1sjβj2-2ϕ∑i=1nuixijy~ijβj, where y~ij=yi-∑k=1,k≠jxikβk. Let σ~j-2=ϕ∑i=1nuixij2+1/sj and μ~j=σ~2ϕ∑i=1nuixijy~ij. Then the full conditional distribution of βj is just the normal distribution N(μ~j,σ~j2). The full conditional distribution of ϕ is(19)fϕ∣x,y,β,s,u,λ12,…,λp2,a,b∝ϕn/2+p-1exp⁡-ϕ∑i=1nuiyi-xiTβ2+∑j=1pλj2sj2. So the full conditional distribution of ϕ is a gamma distribution.The full conditional distribution ofλj2 is(20)fλj2∣x,y,β,ϕ,s,u,λk,k=1,…,p,k≠j2,a,b∝λj2aexp⁡-ϕsj2+bλj2. That is, the full conditional distribution of λj2 is a gamma distribution. At last, the full conditional distributions of a and b are(21)fa∣x,y,β,ϕ,s,u,λ12,…,λp2,b∝baΓap∏j=1pλj2a,fb∣x,y,β,ϕ,s,u,λ12,…,λp2,a∝bap-1exp⁡-b∑j=1pλj2. Obviously, the full conditional distribution of b is just a gamma distribution. However, the full conditional posterior distribution of c does not have a closed form. Fortunately, we note that the full conditional posterior distribution of c is a log-concave function, so the adaptive rejection sampling algorithm [23] can be used to sample from this distribution. The adaptive rejection sampling algorithm is an efficient sampling method from any univariate probability density function which is log-concave.Remark 3.
The Lasso was originally designed for variable selection for the posterior mode ofβj’s can be exactly zeros due to the nature of the laplace prior in (11). However, the posterior draws of βj’s cannot be exact zeros since βj’s are drawn from the continuous posterior distribution. A post hoc thresholding rule may overcome this difficulty. Alternatively, kyung et al. [24] recommended using the credible interval on the posterior mean to guide variable selection.Remark 4.
For variable selection, sometimes, several explanatory variables may be represented by a group of derived input variables. In this case the selection of important variables corresponds to the group of variables. The group Lasso penalty [13, 25] takes the group structure into account and can do variable selection at the group level. For the Bayesian t regression model with adaptive group Lasso penalty, we give the details in the Appendices.


## 4. Simulation Studies

The main goal of this section is to assess the performance of our method (BAT.L) through some simulations with comparison to several Bayesian and non-Bayesian methods. The Bayesian methods include the Bayesian Lasso (BLS.L) [11] and Bayesian adaptive Lasso (BALS.L) [22]. The non-Bayesian methods include Lasso (LASSO) [1] and adaptive Lasso (ALASSO) [20]. The data in the simulation studies are generated by(22)yi=xiTβ+ϵi,i=1,…,n,xi~N(0,Σ), where the covariance matrix Σ is an AR(1) matrix; that is, Σjk=ρ|j-k| with correlation ρ=0.5, and the error distributions of ϵi considered here are the following six situations: The first choice is a standard normal distribution,N(0,1).The second choice is a normal distribution with variance 9,N(0,32).The third choice is at distribution with degree of freedom v=3, t(3).The fourth choice is Laplace distribution,Lp(0,1).The last two choices are the mixture of two normal distributions:0.8N(0,1)+0.2N(0,62) and the mixture of two Laplace distributions: 0.9Lp(0,1)+0.1Lp(0,5).For each simulation study and each choice of the error distribution, we run 50 simulations. In each simulation, we generate a training set with 20 observations, a validation set with 20 observations, and a testing set with 200 observations. The validation set is used to choose the penalty parameters in LASSO and ALASSO. After the penalty parameters are selected, we combine the training set and validation set together to estimateβ. Since BAT.L, BLS.L, and BALS.L do not need the validation set since they estimate the penalty automatically, so we combine the training set and validation set together for estimation. The testing set is used to evaluate the performance for these methods. Within each simulation study, we consider three different settings for β.Simulation 1: β=(3,1.5,0,0,2,0,0,0)T.Simulation 2: β=(5,0,0,0,0,0,0,0)T.Simulation 3: if j=5,10,15,20, βj=5; else, βj=0, j≤20.Simulations 1, 2, and 3 are corresponding to the sparse case, very sparse case, and sparse recovery problem in the predictors, respectively.Note that we intentionally choose error distributions that are different from thet distribution and Laplace distribution to see how the Bayesian t regression adaptive Lasso estimation depends on the error assumption. Our simulation results show that, in terms of parameter estimation accuracy, the Bayesian t adaptive Lasso method still performs well even when this error distribution assumption is violated.Unless otherwise specified, the freedom degreev=4 is used in all simulations for BAT.L. Since the true model is known, we can compute the median of mean absolute deviations (MMAD), that is, median (1/200∑i=1200|xiTβ^-xiTβtrue|), where the median of mean absolute deviations is calculated based on 50 simulations. We present the results for simulation in Table 1. The parentheses in these tables are the standard deviations of the MMAD obtained by 500 bootstrap resamplings of the 50 mean absolute deviations.Table 1
Simulation results of MMAD for different distribution errors†.
Method
N
(
0,1
)
N
(
0,3
)
t
(
3
)
Lp
(
0,1
)
Normal mixture
Laplace mixture
Simulation 1
BAT.L
0.2321 (0.1211)
1.0274 (0.4014)
0.3012 (0.1287)
0.2721 (0.2201)
0.3734 (0.4105)
0.4168 (0.3710)
BLS.L
0.3613 (0.1068)
0.9694 (0.3339)
0.4980 (0.2126)
0.4958 (0.1556)
0.8623 (0.3160)
0.8498 (0.5587)
BALS.L
0.2147 (0.1145)
0.9433 (0.4071)
0.3517 (0.2010)
0.3349 (0.1770)
0.7209 (0.3799)
0.6407 (0.6589)
LASSO
0.3221 (0.1311)
0.9755 (0.4875)
0.4396 (0.2771)
0.4306 (0.1635)
0.8073 (0.3775)
0.7993 (0.5934)
ALASSO
0.2219 (0.1165)
0.9609 (0.4882)
0.3504 (0.2628)
0.3265 (0.1962)
0.8352 (0.4046)
0.6840 (0.6998)
Simulation 2
BAT.L
0.1362 (0.0658)
0.3932 (0.2607)
0.2188 (0.1658)
0.1745 (0.0906)
0.2587 (0.2985)
0.2330 (0.1060)
BLS.L
0.3081 (0.0877)
0.6778 (0.2811)
0.4881 (0.2249)
0.4709 (0.1606)
0.7045 (0.2552)
0.7356 (0.2784)
BALS.L
0.1308 (0.1100)
0.4518 (0.3194)
0.2479 (0.2349)
0.2015 (0.1759)
0.3700 (0.2961)
0.4355 (0.2568)
LASSO
0.1910 (0.1089)
0.6835 (0.3459)
0.3506 (0.2245)
0.3112 (0.1935)
0.5540 (0.3064)
0.4977 (0.3626)
ALASSO
0.0932 (0.1219)
0.4152 (0.3220)
0.2187 (0.2509)
0.1508 (0.2113)
0.2429 (0.3322)
0.3974 (0.3141)
Simulation 3
BAT.L
0.3253 (0.1049)
0.8987(0.4659)
0.3108 (0.1538)
0.3610 (0.3654)
0.3915 (0.1939)
0.4196 (0.0112)
BLS.L
0.6694 (0.1443)
1.7382 (0.3742)
0.8894 (0.3500)
0.9042 (0.2057)
1.6136 (0.6524)
1.3718 (0.0775)
BALS.L
0.3183 (0.1374)
0.9895 (0.3686)
0.5150 (0.2234)
0.4676 (0.1962)
0.9364 (0.5522)
0.6973 (0.0660)
LASSO
0.4469 (0.1449)
1.3959 (0.4807)
0.7231 (0.3199)
0.6462 (0.2530)
1.2599 (0.6710)
1.1223 (0.1315)
ALASSO
0.2803 (0.1463)
0.9267 (0.4087)
0.4601 (0.2306)
0.4513 (0.2290)
0.9204 (0.8101)
0.6814 (0.0987)
†In the parentheses are standard deviations of the MMADs obtained 500 bootstrap resampling. The bold numbers correspond to the smallest MMAD in each category.Simulations show that, in terms of the MMAD, our new regularized method performs better than the other four methods in general, especially for the nonnormal error distributions.


## 5. Land Rent Data

In this section, we demonstrate the performance of BAT.L together with other methods by a real data set.Weisberg [26] reported a data set about land rent data. The data was collected by Douglas Tiffany to study the variation in rent paid in 1977 for agricultural land planted with alfalfa. The variables are average rent per acre planted to alfalfa (y), average rent paid for all tillable land (x1), density of dairy cows (number per square mile) (x2), and proportion of farmland used as pasture (x3) and x4=1 if liming is required to grow alfalfa; x4=0, otherwise. The unit of analysis is a county in Minnesota; the 67 counties with appreciable rented farmland are included. The response variable alfalfa, which has mean 42.1661 and standard deviation 22.5866, is a high protein crop that is suitable feed for dairy cows. It is thought that rent for land planted with alfalfa relative to rent for other agricultural purposes would be higher in areas with a high density of dairy cows and rents would be lower in counties where liming is required, since that would mean additional expense.We first center all the variables so that the intercept is not considered and the predictors are standardized and then use above five regularized methods to fit the whole data set. For the three Bayesian methods, we use the Bayesian credible interval to guide variable selection. We found that the five different methods all selectx1 and x2 as important variables. Furthermore, in order to evaluate the performance of BAT.L together with other methods, we randomly select 40 samples as train data set and the remaining 27 samples as test data set; the mean square errors for test data set are reported in the Table 2 for five different regularized methods.Table 2
MSE of five different regularized methods for land rent data.
Method
BAT.L
LASSO
ALASSO
BLS.L
BALS.L
MSE
94.0764
114.1111
104.3413
109.9957
104.8580As we can see from Table2, the performance of BAT.L outperforms other variable selection procedures for this data set.


## 6. Conclusion Remarks

In this paper, we studied the regularizedt regression model based on Bayesian method. This new method extended the Bayesian Lasso [11] through replacing the least square loss function by the log-likelihood function of t distribution and the Lasso penalty by adaptive Lasso penalty. Bayesian hierarchical models are developed and Gibbs samplers are derived for our new method. Both the simulation studies and the real data analysis show that BAT.L performs better than other existing Bayesian and non-Bayesian methods when the error distribution has heavy tails and/or outliers. We also discussed the Bayesian t regression model with adaptive group Lasso penalty.There is room to improve our methods. In this paper, we treat the dispersion parameterσ in t regression model as constant. In fact, the dispersion parameter σ may be different for different observations, which need us to consider the problem of regularized t regression with varying dispersion model; that is,(23)yi=xiTβ+ɛilog⁡σi=ziTγɛi~t0,σi2,vi=1,…,n, where zi’s are covariates, which constitute, in general, although not necessary, a subset of xi’s. For the varying dispersion t regression model, it is important to consider how to simultaneously do variable selection for both regression coefficients β and γ based on Bayesian method or non-Bayesian method. Furthermore, it also be interesting to consider the problem of variable selection for varying dispersion t regression models with high dimensional covariates. Research in this aspect is ongoing.



---

*Source: 989412-2015-09-20.xml*