The least absolute shrinkage and selection operator (LASSO) was developed by Tibshirani (1996) as an alternative to the ordinary least squares (OLS) method with two objectives in mind. The first was to improve prediction accuracy, and the second was to improve model interpretation by determining a smaller subset of regressors that exhibit the strongest effects. This example presents a fully Bayesian interpretation and implementation of the LASSO that provides estimates for the regression parameters and their variances and provides Bayesian credible intervals for the regression parameters that can guide variable selection.
The LASSO is commonly used to estimate the parameters in the linear regression model
where is the vector of responses, is the overall mean, is the matrix of standardized regressors, and is the vector of independent and identically distributed normal errors with mean 0 and unknown variance . The LASSO estimates of Tibshirani (1996) are the solution to the minimization problem
for some , where .
Tibshirani (1996) suggested that the LASSO estimates can be interpreted as posterior mode estimates when the regression parameters have independent and identical Laplace priors. Park and Casella (2008) consider a fully Bayesian analysis by using a conditional Laplace prior specification of the form
and the noninformative scale-invariant marginal prior on . Conditioning of is important because it guarantees a unimodal full posterior. Park and Casella (2008) also note that any inverse-gamma prior for maintains conjugacy.
Exploiting the fact that the Laplace distribution can be represented as a scale mixture of normal densities with an exponential mixing density, Park and Casella (2008) propose the following hierarchical representation of the full model:
The parameter can be given an independent, flat prior. After you integrate out , the conditional prior on has the desired conditional Laplace distribution.
The Bayesian LASSO parameter can be chosen by using marginal maximum likelihood or an appropriate hyperprior. The example in the next section demonstrates the latter and considers, as suggested by Park and Casella (2008), the class of gamma priors on ,
This example from Park and Casella (2008) fits a Bayesian LASSO model to the diabetes data from Efron et al. (2004). In the original study, statisticians were asked to construct a model that predicted the response variable,
Y, a quantitative measure of disease progression one year after baseline, from 10 covariates:
GLU. It was hoped that the model would produce accurate baseline predictions of response for future patients and that the form of the model would suggest which covariates were important factors in disease progression. The following SAS statements read the data and create the SAS data set
data diabetes; input age sex bmi map tc ldl hdl tch ltg glu y; sex=ifn(sex=2,1,0); datalines; 59 2 32.1 101.00 157 93.2 38.0 4.00 4.8598 87.000 151 48 1 21.6 87.00 183 103.2 70.0 3.00 3.8918 69.000 75 72 2 30.5 93.00 156 93.6 41.0 4.00 4.6728 85.000 141 ... more lines ... 60 2 24.9 99.67 162 106.6 43.0 3.77 4.1271 95.000 132 36 1 30.0 95.00 201 125.2 42.0 4.79 5.1299 85.000 220 36 1 19.6 71.00 250 133.2 97.0 3.00 4.5951 92.000 57 ;
Before specifying the model in the MCMC procedure, you need to standardize the model’s regressors, excluding the indicator variable
Sex. You can use the STDIZE procedure as follows to perform this task:
proc stdize data=diabetes out=std_diabetes; var age bmi map tc ldl hdl tch ltg glu; run;
The following statements specify the Bayesian LASSO in PROC MCMC:
ods graphics on; ods output postintervals=intervals; proc mcmc data=std_diabetes seed=45678 nmc=50000 propcov=quanew monitor=(b0 beta1-beta10 tau1-tau10 sigma2 lasso) outpost=posterior; array D[10,10]; array beta beta1-beta10; array mu0; array data age sex bmi map tc ldl hdl tch ltg glu; begincnst; call identity(D); call zeromatrix(mu0); endcnst; beginnodata; lasso=sqrt(lambda); b=lambda/2; %macro loop; %do k = 1 %to 10; tau&k = exp(omega&k); D[&k,&k]=sigma2*tau&k; %end; %mend loop; %loop; endnodata; call mult(beta, data,xb); parms lambda ; prior lambda ~ gamma(1,scale=.1); parms omega1-omega10; prior omega: ~ expexpon(iscale=b); parms sigma2 1; prior sigma2 ~ igamma(shape = .1, iscale = .1); parms b0 0; prior b0 ~ general(0); parms beta; prior beta ~ mvn(mu0,D); model y ~ normal(b0 + xb,var=sigma2); run;
The ODS OUTPUT statement saves the posterior credible intervals in the SAS data set
Intervals. The NMC= option in the PROC MCMC statement requests 50,000 MCMC iterations,
excluding the burn-in iterations. A large sample is used because the posterior samples are highly autocorrelated. The PROPCOV= option in the PROC MCMC statement requests that the quasi-Newton method be used
in constructing the initial covariance matrix for the Metropolis-Hastings algorithm. The OUTPOST= option saves the posterior samples in the data set
The next four statements create arrays that are used in the model. The array D is the covariance matrix for the regression parameters
Beta10. The array Beta is the vector of the regression parameters
Beta10. The array Mu0 is the mean vector for the prior distribution of the regression
Beta10. The array Data is the matrix of regressors, excluding the intercept.
The BEGINCNST and ENDCNST statements define a statement block within which PROC MCMC processes the programming statements only during the setup stage of the simulation. You can use the BEGINCNST and ENDCNST statement block to initialize the matrices D and Mu0. D is initially set to an identity matrix, and Mu0 is initialized as a zero vector.
The BEGINNODATA and ENDNODATA statements define a block within which PROC MCMC processes the programming statements without stepping through the entire data set. The programming statements are executed only twice: at the first and last observations of the data set. Within this statement block, the parameters
b are defined. The macro %LOOP repopulates the matrix D. The purpose of the parameters
Omega10 and their relationship with the parameters
Tau10 are explained later.
The next statement uses the MULT CALL routine to define the matrix XB, which contains the product of the regressors and the regression parameters
Beta10. That is, it contains the linear predictor, excluding the intercept.
The following block of statements declares the model parameters and assigns prior distribution to them. The parameter
Lambda, which represents
, is specified to have a gamma distribution.
Omega10 are specified to have exponential exponential distributions. The parameters
have exponential distributions, but modeling these parameters directly can cause convergence problems. Instead, the parameters
Omega10 are modeled directly, and within the macro %LOOP the parameters
Tau10, which represent , are defined as being the exponential of
Omega10, respectively. The parameter
Sigma2, which represents , is specified to have an inverse-gamma distribution. The parameter
B0, which represent , is specified to have an improper uniform distribution. The parameter vector Beta, which represents , is specified to have a multivariate normal distribution with mean equal to 0 and variance matrix equal to D.
Finally, the MODEL statement specifies that the response variable
Y have a normal distribution.
Output 1 shows that the Monte Carlo standard errors (MCSE) of each parameter are small relative to the posterior standard deviations (SD). A small MCSE/SD ratio indicates that the Markov chain has stabilized and that the mean estimates do not vary much over time.
Output 1: Monte Carlo Standard Errors
The MCMC Procedure
|Monte Carlo Standard Errors|
Output 2 shows the “Effective Sample Sizes” table. The autocorrelation times for the parameters range from 1.59 to 62.83, and most of the efficiency rates are low. These results account for the relatively small effective sample sizes, given a nominal sample size of 50,000.
Output 2: Effective Sample Sizes
|Effective Sample Sizes|
The following SAS statements use the OUTPOST data set
Posterior and the ODS OUTPUT data set
Intervals to generate a table of the Bayesian
LASSO parameter estimates, which are the modes of the posterior samples for
Beta10, and their
respective 95% HPD intervals:
proc means data=posterior mode; var b0 beta1-beta10; output out=parameters(drop=_TYPE_ _FREQ_) mode(b0 beta1-beta10)=b0 beta1-beta10; run; proc transpose data=parameters out=parameters; run;
data parameters; length parameter $ 6; set parameters(rename=(col1=mode _NAME_=Parameter)); label Parameter=; index=_N_; run;
proc sort data=parameters out=parameters; by parameter; run; proc sort data=intervals out=intervals; by parameter; run;
data parameters(where=(~missing(mode))); merge parameters intervals; by parameter; label parameter="Parameter" mode="Mode"; run;
proc sort data=parameters out=parameters; by index; run;
proc print data=parameters noobs label; var parameter mode hpdlower hpdupper; run;
Output 3 shows that the HPD intervals for the parameters
Beta10 all contain 0. Unlike what happens in the frequentist version of the LASSO, regression
parameters are not set to 0, so the inclusion of 0 in the HPD interval is the only indication that a variable is a candidate for exclusion from the model. Based on this criterion, the variables
GLU are the leading
candidates for exclusion from the model.
Output 3: Bayesian LASSO Parameter Estimates and 95% HPD Intervals