12.2 Regularization
In this century, the amount of available data continues to grow. This means we have access to more covariates for prediction, and we can also generate additional inputs to enhance the predictive power of our models. As a result, we often encounter wide datasets, where the number of inputs may exceed the number of observations. Even in modest settings, we might have hundreds of inputs, and we can use them to identify which ones contribute to accurate predictions. However, we generally avoid using all inputs at once due to the risk of overfitting. Thus, we require a class of methods for input selection or regularization.
In the standard linear regression setting,
\[ \mathbf{y} = \mathbf{i}_N \beta_0 + \mathbf{W}\boldsymbol{\beta} + \boldsymbol{\mu}, \]
where \(\mathbf{i}_N\) is an \(N\)-dimensional vector of ones, \(\mathbf{W}\) is the \(N \times K\) design matrix of inputs, and \(\boldsymbol{\mu} \sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I}_N)\), there has been extensive development of techniques aimed at regularization within the Frequentist inferential framework. These include methods such as Ridge regression (Hoerl and Kennard 1970); discrete subset selection techniques like best subset selection (Furnival and Wilson 1974a), forward selection, and backward stepwise selection (Hastie, Tibshirani, and Friedman 2009); as well as continuous subset selection approaches such as the LASSO (Tibshirani 1996), the elastic net (Zou and Hastie 2005), and OSCAR (Bondell and Reich 2008).
It is important to note, however, that Ridge regression does not perform variable selection; rather, it shrinks coefficient estimates toward zero without setting them exactly to zero.
Ridge regression and the LASSO can be viewed as special cases of a more general class of estimators known as Bridge regression (Fu 1998), which also admits a Bayesian interpretation. Consider the following penalized least squares criterion in the linear regression setting:
\[ \hat{\boldsymbol{\beta}}^{\text{Bridge}} = \arg\min_{\boldsymbol{\beta}} \left\{ \sum_{i=1}^N \left( y_i - \beta_0 - \sum_{k=1}^K \tilde{w}_{ik} \beta_k \right)^2 + \lambda \sum_{k=1}^K |\beta_k|^q \right\}, \quad q \geq 0, \]
where \(\tilde{w}_{ik}\) denotes the standardized inputs. Standardizing inputs is important in variable selection problems to avoid issues caused by differences in scale; otherwise, variables with larger magnitudes will be penalized less and disproportionately influence the regularization path.
Interpreting \(|\beta_k|^q\) as proportional to the negative log-prior density of \(\beta_k\), the penalty shapes the contours of the prior distribution on the parameters. Specifically:
- \(q = 0\) corresponds to best subset selection, where the penalty counts the number of nonzero coefficients.
- \(q = 1\) yields the LASSO, which corresponds to a Laplace (double-exponential) prior.
- \(q = 2\) yields ridge regression, which corresponds to a Gaussian prior.
In this light, best subset selection, the LASSO, and ridge regression can be viewed as maximum a posteriori (MAP) estimators under different priors centered at zero (Park and Casella 2008). However, they are not Bayes estimators in the strict sense, since Bayes estimators are typically defined as the posterior mean. While ridge regression coincides with the posterior mean under a Gaussian prior (Ishwaran and Rao 2005), the LASSO and best subset selection yield posterior modes.
This distinction is important because the Bayesian framework naturally incorporates regularization through the use of proper priors, which helps mitigate overfitting. Specifically, when proper shrinkage priors are used, the posterior balances data likelihood and prior information, thus controlling model complexity.
Furthermore, empirical Bayes methods, where the marginal likelihood is optimized, or cross-validation can be used to estimate the scale parameter of the prior covariance matrix for the regression coefficients. This scale parameter, in turn, determines the strength of regularization in ridge regression.
Note that regularization introduces bias into the parameter estimates because it constrains the model, shrinking the location parameters toward zero. However, it substantially reduces variance, as the estimates are prevented from varying excessively across samples. As a result, the mean squared error (MSE) of the estimates, which equals the sum of the squared bias and the variance, is often lower for regularization methods compared to ordinary least squares (OLS), which remains unbiased under the classical assumptions. This trade-off is particularly important when the goal is to identify causal effects, where unbiasedness may be preferred over predictive accuracy (see Chapter 13).
12.2.1 Bayesian LASSO
Given the popularity of the LASSO as a variable selection technique, we present its Bayesian formulation in this subsection (Park and Casella 2008). The Gibbs sampler for the Bayesian LASSO exploits the representation of the Laplace distribution as a scale mixture of normals. This leads to the following hierarchical representation of the model:
\[ \begin{aligned} \mathbf{y} \mid \beta_0, \boldsymbol{\beta}, \sigma^2, \mathbf{W} &\sim \mathcal{N}(\mathbf{i}_N \beta_0 + \mathbf{W} \boldsymbol{\beta}, \sigma^2 \mathbf{I}_N), \\ \boldsymbol{\beta} \mid \sigma^2, \tau_1^2, \dots, \tau_K^2 &\sim \mathcal{N}(\mathbf{0}_K, \sigma^2 \mathbf{D}_{\tau}), \\ \tau_1^2, \dots, \tau_K^2 &\sim \prod_{k=1}^K \frac{\lambda^2}{2} \exp\left\{ -\frac{\lambda^2}{2} \tau_k^2 \right\}, \\ \sigma^2 &\sim \frac{1}{\sigma^2}, \\ \beta_0 &\sim c, \end{aligned} \]
where \(\mathbf{D}_{\tau} = \operatorname{diag}(\tau_1^2, \dots, \tau_K^2)\) and \(c\) is a constant.
After integrating out \(\tau_1^2, \dots, \tau_K^2\), the conditional prior of \(\boldsymbol{\beta} \mid \sigma^2\) is:
\[ \pi(\boldsymbol{\beta} \mid \sigma^2) = \prod_{k=1}^K \frac{\lambda}{2 \sqrt{\sigma^2}} \exp\left\{ -\frac{\lambda}{\sqrt{\sigma^2}} |\beta_k| \right\}, \]
which implies that the log-prior is proportional to \(\lambda \sum_{k=1}^K |\beta_k|\), matching the penalty term in the LASSO optimization problem.
The conditional posterior distributions for the Gibbs sampler are (Park and Casella 2008):
\[ \begin{aligned} \boldsymbol{\beta} \mid \sigma^2, \tau_1^2, \dots, \tau_K^2, \tilde{\mathbf{W}}, \tilde{\mathbf{y}} &\sim \mathcal{N}(\boldsymbol{\beta}_n, \sigma^2 \mathbf{B}_n), \\ \sigma^2 \mid \boldsymbol{\beta}, \tau_1^2, \dots, \tau_K^2, \tilde{\mathbf{W}}, \tilde{\mathbf{y}} &\sim \text{Inverse-Gamma}(\alpha_n/2, \delta_n/2), \\ 1/\tau_k^2 \mid \boldsymbol{\beta}, \sigma^2 &\sim \text{Inverse-Gaussian}(\mu_{kn}, \lambda_n), \\ \beta_0 \mid \sigma^2, \tilde{\mathbf{W}}, \tilde{\mathbf{y}} &\sim \mathcal{N}(\bar{y}, \sigma^2 / N), \end{aligned} \]
where:
\[ \begin{aligned} \boldsymbol{\beta}_n &= \mathbf{B}_n \tilde{\mathbf{W}}^{\top} \tilde{\mathbf{y}}, \\ \mathbf{B}_n &= \left( \tilde{\mathbf{W}}^{\top} \tilde{\mathbf{W}} + \mathbf{D}_{\tau}^{-1} \right)^{-1}, \\ \alpha_n &= (N - 1) + K, \\ \delta_n &= (\tilde{\mathbf{y}} - \tilde{\mathbf{W}} \boldsymbol{\beta})^{\top} (\tilde{\mathbf{y}} - \tilde{\mathbf{W}} \boldsymbol{\beta}) + \boldsymbol{\beta}^{\top} \mathbf{D}_{\tau}^{-1} \boldsymbol{\beta}, \\ \mu_{kn} &= \sqrt{ \frac{ \lambda^2 \sigma^2 }{ \beta_k^2 } }, \\ \lambda_n &= \lambda^2, \end{aligned} \]
and \(\tilde{\mathbf{W}}\) is the matrix of standardized inputs, and \(\tilde{\mathbf{y}}\) is the centered response vector.
Note that the posterior distribution of \(\boldsymbol{\tau}\) depends on the data through \(\boldsymbol{\beta}\) and \(\sigma^2\), which is a typical feature of hierarchical models. In this formulation, we can interpret \(\tau_k\) as local shrinkage parameters, while \(\lambda\) acts as a global shrinkage parameter. Higher values of \(\tau_k\) and \(\lambda\) imply stronger shrinkage toward zero. Park and Casella (2008) propose two approaches for specifying the global shrinkage parameter: empirical Bayes estimation or a fully Bayesian hierarchical specification, where \(\lambda^2\) is assigned a Gamma prior.
We should acknowledge that the Bayesian LASSO is more computationally expensive than the Frequentist LASSO. However, it provides credible intervals for the parameters automatically. In contrast, obtaining standard errors in the Frequentist LASSO is more challenging, particularly for parameters estimated to be exactly zero (Kyung et al. 2010).
Example: Simulation exercise to study the Bayesian LASSO performance
We simulate the process
\[\begin{equation*}
y_i = \beta_0 + \sum_{k=1}^{10} \beta_k w_{ik} + \mu_i,
\end{equation*}\]
where \(\beta_k \sim \mathcal{U}(-3, 3)\), \(\mu_i \sim \mathcal{N}(0, 1)\), and \(w_{ik} \sim \mathcal{N}(0, 1)\), for \(i = 1, 2, \dots, 500\).
Additionally, we generate 90 extra potential inputs from a standard normal distribution, which are included in the model specification. Our goal is to assess whether the Bayesian LASSO can successfully identify the truly relevant inputs.
We use the bayesreg package in R to perform the Bayesian LASSO, using 5,000 MCMC draws and 1,000 burn-in iterations. The following code illustrates the simulation exercise and compares the posterior means with the true population values.
The summary of the fit and the plot comparing the population parameters with the posterior means show that the Bayesian LASSO is able to identify the variables that are relevant in the data generating process.
In Exercise 1, we propose programming the Gibbs sampler from scratch, assuming a hierarchical structure for the global shrinkage parameter, and comparing the results with those obtained using the monomvn package.
## Warning: package 'bayesreg' was built under R version 4.5.3## Loading required package: pgdraw## Warning: package 'pgdraw' was built under R version 4.5.3## Loading required package: doParallel## Loading required package: foreach##
## Attaching package: 'foreach'## The following objects are masked from 'package:purrr':
##
## accumulate, when## Loading required package: iterators## Loading required package: parallel# Parameters
n <- 500 # sample size
p <- 100 # number of predictors
s <- 10 # number of nonzero coefficients
# Generate design matrix
X <- matrix(rnorm(n * p), nrow = n, ncol = p)
# True beta: first s coefficients are non-zero, rest are zero
beta_true <- c(runif(s, -3, 3), rep(0, p - s))
# Generate response with some noise
sigma <- 1
y <- X %*% beta_true + rnorm(n, sd = sigma)
df <- data.frame(X,y)
### Using bayesreg ###
# Fit the model
fit <- bayesreg::bayesreg(y ~ X, data = df, model = "gaussian", prior = "lasso",
n.samples = 5000, burnin = 1000)
# Check summary
summary(fit)## ==========================================================================================
## | Bayesian Penalised Regression Estimation ver. 1.3 |
## | (c) Enes Makalic, Daniel F Schmidt. 2016-2024 |
## ==========================================================================================
## Bayesian Gaussian lasso regression Number of obs = 500
## Number of vars = 100
## MCMC Samples = 5000 std(Error) = 1.0277
## MCMC Burnin = 1000 R-squared = 0.9759
## MCMC Thinning = 5 WAIC = 762.14
##
## -------------+-----------------------------------------------------------------------------
## Parameter | mean(Coef) std(Coef) [95% Cred. Interval] tStat Rank ESS
## -------------+-----------------------------------------------------------------------------
## X1 | -0.06875 0.04821 -0.16883 0.02258 -1.426 19 * 5000
## X2 | -2.21521 0.05069 -2.30653 -2.10919 -43.697 3 ** 4761
## X3 | -2.64365 0.05053 -2.73372 -2.53560 -52.322 1 ** 5000
## X4 | -0.74195 0.04853 -0.83361 -0.64557 -15.287 9 ** 4803
## X5 | 1.81544 0.05074 1.71099 1.91046 35.779 6 ** 5000
## X6 | -2.12284 0.04932 -2.20982 -2.01797 -43.042 3 ** 3965
## X7 | -2.10257 0.05163 -2.19467 -1.99403 -40.721 5 ** 5000
## X8 | 0.88326 0.04833 0.78502 0.97606 18.276 8 ** 4729
## X9 | 1.27162 0.05363 1.16420 1.37073 23.712 7 ** 5000
## X10 | -2.56227 0.05182 -2.65710 -2.45138 -49.447 2 ** 5000
## X11 | -0.02981 0.04570 -0.12172 0.05744 -0.652 46 4863
## X12 | -0.02402 0.05045 -0.12738 0.07339 -0.476 59 4648
## X13 | -0.06499 0.04941 -0.16339 0.02908 -1.315 23 * 5000
## X14 | -0.02842 0.04492 -0.11956 0.05649 -0.633 51 4655
## X15 | -0.04610 0.04756 -0.13840 0.04846 -0.969 34 5000
## X16 | -0.08182 0.04881 -0.18144 0.00919 -1.676 15 * 5000
## X17 | -0.02080 0.04694 -0.11560 0.06777 -0.443 65 5000
## X18 | -0.01180 0.04736 -0.10387 0.08468 -0.249 78 5000
## X19 | 0.06335 0.04797 -0.02838 0.15812 1.321 22 * 4966
## X20 | 0.01945 0.04766 -0.07438 0.11465 0.408 65 5000
## X21 | 0.03020 0.04709 -0.06154 0.12178 0.641 51 4847
## X22 | -0.02106 0.04768 -0.11712 0.06890 -0.442 65 4812
## X23 | 0.03314 0.04842 -0.06058 0.13271 0.685 46 5000
## X24 | -0.00294 0.04321 -0.09083 0.08217 -0.068 91 5000
## X25 | -0.00363 0.04209 -0.08651 0.08129 -0.086 91 5000
## X26 | -0.01302 0.04673 -0.10832 0.07703 -0.279 78 5000
## X27 | -0.05896 0.04782 -0.15402 0.03083 -1.233 24 * 4929
## X28 | 0.02094 0.04556 -0.06940 0.11067 0.460 78 4885
## X29 | 0.08386 0.05103 -0.01373 0.18503 1.643 15 * 4742
## X30 | 0.04346 0.04821 -0.04401 0.14059 0.901 37 5000
## X31 | 0.02826 0.04499 -0.06108 0.11469 0.628 51 5000
## X32 | -0.08035 0.04723 -0.17570 0.00641 -1.701 14 * 4615
## X33 | -0.02013 0.04535 -0.11134 0.06457 -0.444 65 4886
## X34 | -0.01977 0.04655 -0.11164 0.07098 -0.425 65 5000
## X35 | -0.01973 0.04582 -0.11197 0.06825 -0.430 78 4921
## X36 | 0.00136 0.04488 -0.08840 0.09101 0.030 91 4835
## X37 | -0.04899 0.04620 -0.14280 0.03944 -1.060 30 4851
## X38 | -0.00764 0.04769 -0.10377 0.08659 -0.160 91 4871
## X39 | 0.01468 0.04845 -0.08177 0.11165 0.303 78 4682
## X40 | -0.02047 0.04796 -0.11553 0.07224 -0.427 65 5000
## X41 | -0.00983 0.04429 -0.09744 0.07816 -0.222 91 4964
## X42 | 0.03616 0.04612 -0.04864 0.13230 0.784 45 5000
## X43 | -0.05537 0.04737 -0.14952 0.03634 -1.169 25 * 4910
## X44 | -0.01806 0.04684 -0.11602 0.07111 -0.386 78 5000
## X45 | 0.07316 0.04880 -0.01788 0.17126 1.499 18 * 4908
## X46 | -0.03822 0.04543 -0.12872 0.04861 -0.841 41 4983
## X47 | -0.05053 0.05029 -0.15086 0.04607 -1.005 32 5000
## X48 | 0.02341 0.04756 -0.06766 0.11902 0.492 59 5000
## X49 | -0.01178 0.04839 -0.11012 0.08143 -0.243 78 5000
## X50 | -0.02063 0.04461 -0.10857 0.06700 -0.462 78 4859
## X51 | -0.05355 0.04487 -0.14254 0.03216 -1.193 27 * 4985
## X52 | -0.00236 0.04744 -0.09893 0.09203 -0.050 91 5000
## X53 | -0.04902 0.04995 -0.15113 0.04615 -0.981 32 5000
## X54 | -0.04526 0.04973 -0.14500 0.04718 -0.910 34 5000
## X55 | 0.04061 0.04756 -0.05077 0.13598 0.854 37 5000
## X56 | 0.03060 0.04590 -0.05853 0.12270 0.667 46 4939
## X57 | 0.01757 0.04578 -0.07160 0.10652 0.384 65 5000
## X58 | 0.00808 0.04505 -0.08300 0.09572 0.179 78 5000
## X59 | -0.05236 0.04604 -0.14563 0.03404 -1.137 27 * 4821
## X60 | -0.00138 0.04555 -0.08964 0.09040 -0.030 78 5000
## X61 | 0.15159 0.04737 0.05801 0.24293 3.200 10 ** 5000
## X62 | 0.03098 0.04698 -0.06146 0.12343 0.659 59 4961
## X63 | -0.08365 0.05110 -0.18625 0.01110 -1.637 17 * 4971
## X64 | -0.00508 0.04725 -0.09838 0.08895 -0.108 91 5000
## X65 | -0.03447 0.04523 -0.12860 0.05332 -0.762 41 5000
## X66 | -0.02369 0.04672 -0.11848 0.06352 -0.507 59 4976
## X67 | 0.01859 0.04568 -0.06758 0.11299 0.407 65 5000
## X68 | -0.03084 0.04899 -0.12974 0.06325 -0.629 51 5000
## X69 | -0.05080 0.04634 -0.14209 0.03896 -1.096 30 4769
## X70 | -0.06600 0.04752 -0.16192 0.02268 -1.389 21 * 5000
## X71 | -0.04150 0.04583 -0.13085 0.04616 -0.906 37 5000
## X72 | -0.00590 0.04698 -0.09916 0.08813 -0.126 91 5000
## X73 | -0.03944 0.04873 -0.13602 0.05286 -0.809 41 4675
## X74 | 0.02695 0.04484 -0.05965 0.11770 0.601 65 5000
## X75 | 0.06592 0.04798 -0.02546 0.16368 1.374 20 * 5000
## X76 | -0.01316 0.04798 -0.11170 0.07994 -0.274 78 4460
## X77 | 0.04212 0.04767 -0.04704 0.14001 0.884 37 5000
## X78 | -0.02927 0.04850 -0.12789 0.06420 -0.604 51 4268
## X79 | -0.01884 0.04596 -0.11055 0.07256 -0.410 65 4310
## X80 | -0.05864 0.04894 -0.15897 0.03057 -1.198 25 * 5000
## X81 | 0.03340 0.04747 -0.05600 0.12954 0.704 46 5000
## X82 | -0.02590 0.04541 -0.11722 0.06376 -0.570 51 5000
## X83 | -0.04793 0.04793 -0.14349 0.04350 -1.000 34 5000
## X84 | -0.09502 0.05175 -0.19893 0.00392 -1.836 11 * 4858
## X85 | -0.04022 0.04675 -0.13366 0.04678 -0.860 41 5000
## X86 | -0.05434 0.04737 -0.14929 0.03387 -1.147 27 * 5000
## X87 | -0.00652 0.04611 -0.09906 0.08676 -0.141 91 4693
## X88 | 0.02042 0.04787 -0.07747 0.11476 0.427 65 5000
## X89 | -0.08417 0.04705 -0.17721 0.00435 -1.789 13 * 4504
## X90 | 0.01814 0.04651 -0.07184 0.11185 0.390 65 5000
## X91 | -0.00288 0.04629 -0.09494 0.09083 -0.062 91 5000
## X92 | 0.02695 0.04581 -0.06194 0.11691 0.588 51 4482
## X93 | -0.10024 0.05358 -0.20372 0.00530 -1.871 11 * 5000
## X94 | -0.02224 0.04393 -0.11143 0.06421 -0.506 59 4706
## X95 | 0.00423 0.04656 -0.08699 0.10026 0.091 78 4684
## X96 | 0.01956 0.04636 -0.06928 0.11135 0.422 65 5000
## X97 | -0.03601 0.04889 -0.13368 0.05712 -0.737 51 4771
## X98 | 0.02668 0.04744 -0.06456 0.12291 0.562 59 5000
## X99 | -0.03488 0.04789 -0.13177 0.05649 -0.728 46 5000
## X100 | -0.01170 0.04445 -0.09897 0.07378 -0.263 78 4772
## (Intercept) | 0.06173 0.05131 -0.03743 0.16471 . . .
## -------------+-----------------------------------------------------------------------------# Extract posterior means of beta
beta_post_mean <- rowMeans(fit$beta)
# Compare true vs estimated
plot(beta_true, beta_post_mean, pch = 19, col = "steelblue",
xlab = "True beta", ylab = "Posterior mean of beta",
main = "Bayesian LASSO Shrinkage")
abline(0, 1, col = "red", lty = 2)
12.2.2 Stochastic search variable selection
Another well-known Bayesian strategy for variable selection in the presence of a large set of regressors (inputs) is stochastic search variable selection (SSVS) (E. I. George and McCulloch 1993; E. George and McCulloch 1997). SSVS is a particular case of the broader class of spike-and-slab priors, in which the prior distribution for the location parameters is specified as a hierarchical mixture that captures the uncertainty inherent in variable selection problems (Ishwaran and Rao 2005).
The hierarchical structure of the model is given by:
\[ \begin{aligned} \mathbf{y} \mid \beta_0, \boldsymbol{\beta}, \sigma^2, \mathbf{W} &\sim \mathcal{N}(\mathbf{i}_N \beta_0 + \mathbf{W} \boldsymbol{\beta}, \sigma^2 \mathbf{I}_N), \\ \boldsymbol{\beta} \mid \sigma^2, \boldsymbol{\gamma} &\sim \mathcal{N}(\mathbf{0}_K, \sigma^2 \mathbf{D}_\gamma \mathbf{R} \mathbf{D}_\gamma), \\ \sigma^2 &\sim \text{Inverse-Gamma}\left(\frac{v}{2}, \frac{v \lambda_\gamma}{2}\right), \\ \gamma_k &\sim \text{Bernoulli}(p_k), \end{aligned} \]
where \(p_k\) is the prior inclusion probability of regressor \(w_k\), that is, \(P(\gamma_k = 1) = 1 - P(\gamma_k = 0) = p_k\), \(\mathbf{R}\) is a correlation matrix, and \(\mathbf{D}_\gamma\) is a diagonal matrix whose \((k,k)\)-th element is defined as:
\[ (\mathbf{D}_\gamma)_{kk} = \begin{Bmatrix} v_{0k}, & \text{if } \gamma_k = 0, \\ v_{1k}, & \text{if } \gamma_k = 1 \end{Bmatrix}. \]
This formulation implies that:
\[ \beta_k \sim (1 - \gamma_k) \, \mathcal{N}(0, v_{0k}) + \gamma_k \, \mathcal{N}(0, v_{1k}), \]
where \(v_{0k}\) and \(v_{1k}\) are chosen such that \(v_{0k}\) is small and \(v_{1k}\) is large. Therefore, when the data favors \(\gamma_k = 0\), the corresponding \(\beta_k\) is shrunk toward zero, effectively excluding input \(k\) from the model. In this sense, \(\mathcal{N}(0, v_{0k})\) is a spike prior concentrated at zero, while \(\mathcal{N}(0, v_{1k})\) is a diffuse slab prior.
A critical aspect of SSVS is the choice of the hyperparameters \(v_{0k}\) and \(v_{1k}\), as they determine the amount of shrinkage applied to the regression coefficients (see E. I. George and McCulloch (1993) and E. George and McCulloch (1997) for details).
The assumption \(\gamma_k \sim \text{Bernoulli}(p_k)\) implies that the prior on the inclusion indicators is given by:
\[ \pi(\boldsymbol{\gamma}) = \prod_{k=1}^K p_k^{\gamma_k} (1 - p_k)^{1 - \gamma_k}. \]
This means that the inclusion of input \(k\) is independent of the inclusion of any other input \(j \neq k\). A common choice is the uniform prior \(\pi(\boldsymbol{\gamma}) = 2^{-K}\), which corresponds to setting \(p_k = 1/2\), giving each regressor an equal chance of being included (Ishwaran and Rao 2005).
A practical choice for the correlation matrix is to set \(\mathbf{R} \propto (\tilde{\mathbf{W}}^{\top} \tilde{\mathbf{W}})^{-1}\) (E. I. George and McCulloch 1993). Regarding the hyperparameters \(v\) and \(\lambda_\gamma\), it is helpful to interpret \(\lambda_\gamma\) as a prior estimate of \(\sigma^2\), and \(v\) as the prior sample size associated with this estimate. In the absence of prior information, E. I. George and McCulloch (1997) recommend setting \(\lambda_\gamma\) equal to the least squares estimate of the variance from the saturated model, that is, the model including all regressors, and \(v\) to a small number, for instance, 0.01.
The conditional posterior distributions for the Gibbs sampler are (E. I. George and McCulloch 1993):
\[ \begin{aligned} \boldsymbol{\beta} \mid \sigma^2, \gamma_1, \dots, \gamma_K, \tilde{\mathbf{W}}, \tilde{\mathbf{y}} &\sim N(\boldsymbol{\beta}_n, \mathbf{B}_n), \\ \sigma^2 \mid \boldsymbol{\beta}, \gamma_1, \dots, \gamma_K, \tilde{\mathbf{W}}, \tilde{\mathbf{y}} &\sim \text{Inverse-Gamma}(\alpha_n/2, \delta_n/2), \\ \gamma_k \mid \boldsymbol{\beta}, \sigma^2 &\sim \text{Bernoulli}(p_{kn}), \end{aligned} \]
where:
\[ \begin{aligned} \boldsymbol{\beta}_n &= \sigma^{-2} \mathbf{B}_n \tilde{\mathbf{W}}^{\top} \tilde{\mathbf{y}}, \\ \mathbf{B}_n &= \left(\sigma^{-2} \tilde{\mathbf{W}}^{\top} \tilde{\mathbf{W}} + \mathbf{D}_{\gamma}^{-1}\mathbf{R}^{-1}\mathbf{D}_{\gamma}^{-1} \right)^{-1}, \\ \alpha_n &= N + v, \\ \delta_n &= (\tilde{\mathbf{y}} - \tilde{\mathbf{W}} \boldsymbol{\beta})^{\top} (\tilde{\mathbf{y}} - \tilde{\mathbf{W}} \boldsymbol{\beta}) + v\lambda_{\gamma}, \\ p_{kn} &= \frac{\pi(\boldsymbol{\beta}\mid \boldsymbol{\gamma}_{-k},\gamma_k=1)\times p_k}{\pi(\boldsymbol{\beta}\mid \boldsymbol{\gamma}_{-k},\gamma_k=1)\times p_k+\pi(\boldsymbol{\beta}\mid \boldsymbol{\gamma}_{-k},\gamma_k=0)\times (1-p_k)}, \end{aligned} \]
where \(\tilde{\mathbf{W}}\) is the matrix of standardized inputs, \(\tilde{\mathbf{y}}\) is the centered response vector, \(\boldsymbol{\gamma}_{-k}\) denotes the vector composed of \(\gamma_1, \dots, \gamma_K\) excluding \(\gamma_k\), and \(\pi(\boldsymbol{\beta} \mid \boldsymbol{\gamma}_{-k}, \gamma_k = \delta)\) is the posterior density of \(\boldsymbol{\beta}\) evaluated at \(\boldsymbol{\gamma}_{-k}\) and \(\gamma_k = \delta\), where \(\delta \in \{0,1\}\).
In general, it is wise to consider the inclusion of regressors jointly due to potential correlations among them; that is, the marginal frequency of \(\gamma_k = 1\) should be interpreted with caution. SSVS is more effective at identifying a good set of potential models rather than selecting a single best model.
Example: Simulation exercise to study SSVS performance
Let’s use the simulation setting from the previous example to evaluate the performance of SSVS in uncovering the data-generating process. In particular, we use the BoomSpikeSlab package to implement this example.
The analysis is performed using 5,000 posterior draws and the default prior. However, the package allows the user to modify the default prior via the SpikeSlabPrior function.
The results show that the posterior inclusion probabilities for regressors 2 through 10 are 100%, and the model with the highest posterior probability (94%) includes all of these nine variables. However, the true data-generating process, which also includes regressor 1, receives a posterior model probability of 0%. This is because the population coefficient of this regressor is essentially zero. The plot comparing the posterior means with the true population parameters indicates good performance of SSVS. In general, Bayesian methods for variable selection perform well, and the choice of the most suitable method largely depends on the prior specification (O’Hara and Sillanpää 2009).
####### Stochastic search variable selection #######
rm(list = ls()); set.seed(10101)
library(BoomSpikeSlab)## Warning: package 'BoomSpikeSlab' was built under R version 4.5.3## Loading required package: Boom## Warning: package 'Boom' was built under R version 4.5.3##
## Attaching package: 'Boom'## The following objects are masked from 'package:brms':
##
## ddirichlet, rdirichlet## The following objects are masked from 'package:LaplacesDemon':
##
## ddirichlet, dinvgamma, dmvn, rdirichlet, rinvgamma, rmvn## The following objects are masked from 'package:MCMCpack':
##
## ddirichlet, dinvgamma, rdirichlet, rinvgamma## The following object is masked from 'package:coda':
##
## thin## The following object is masked from 'package:sirt':
##
## rmvn## The following object is masked from 'package:stats':
##
## rWishart##
## Attaching package: 'BoomSpikeSlab'## The following object is masked from 'package:stats':
##
## knotslibrary(dplyr)
library(tibble)
# Parameters
n <- 500 # sample size
k <- 100 # number of predictors
s <- 10 # number of non-zero coefficients
# Generate design matrix
X <- matrix(rnorm(n * k), nrow = n, ncol = k)
# True beta: first s coefficients are non-zero, rest are zero
beta_true <- c(runif(s, -3, 3), rep(0, k - s))
# Generate response with some noise
sigma <- 1
y <- X %*% beta_true + rnorm(n, sd = sigma)
df <- data.frame(X,y)
### Using BoomSpikeSlab ###
#Scale regressors
W <- scale(X); yh <- y - mean(y)
prior <- SpikeSlabPrior(W, yh,
expected.model.size = ncol(W)/2, # expect 50 nonzero predictors
prior.df = .01, # weaker prior than the default
prior.information.weight = .01,
diagonal.shrinkage = 0) # shrink to zero
niter <- 5000
######Estimate model########
SSBoomNew <- lm.spike(yh ~ W - 1, niter = niter, prior = prior)## =-=-=-=-= Iteration 0 Wed Jun 17 15:23:00 2026
## =-=-=-=-=
## =-=-=-=-= Iteration 500 Wed Jun 17 15:23:00 2026
## =-=-=-=-=
## =-=-=-=-= Iteration 1000 Wed Jun 17 15:23:01 2026
## =-=-=-=-=
## =-=-=-=-= Iteration 1500 Wed Jun 17 15:23:01 2026
## =-=-=-=-=
## =-=-=-=-= Iteration 2000 Wed Jun 17 15:23:01 2026
## =-=-=-=-=
## =-=-=-=-= Iteration 2500 Wed Jun 17 15:23:02 2026
## =-=-=-=-=
## =-=-=-=-= Iteration 3000 Wed Jun 17 15:23:02 2026
## =-=-=-=-=
## =-=-=-=-= Iteration 3500 Wed Jun 17 15:23:03 2026
## =-=-=-=-=
## =-=-=-=-= Iteration 4000 Wed Jun 17 15:23:03 2026
## =-=-=-=-=
## =-=-=-=-= Iteration 4500 Wed Jun 17 15:23:04 2026
## =-=-=-=-=Models <- SSBoomNew$beta != 0
PIP <- colMeans(SSBoomNew$beta != 0)
# Convert the logical matrix to a data frame and then to a tibble
df <- as.data.frame(Models); df_tbl <- as_tibble(df)
# Count identical rows
row_counts <- df_tbl %>% count(across(everything()), name = "frequency") %>% arrange(desc(frequency))
sum(row_counts[1:100,101])## [1] 3948# Ensure your vector and matrix are logical
trueModel <- c(rep(1, 10), rep(0, 90)) == 1 # convert to logical if needed
# Assume your matrix is named 'mat'
matching_rows <- apply(row_counts[,-101], 1, function(row) all(row == trueModel))
# Get indices (row numbers) where the match is TRUE
row_counts[which(matching_rows), 101]## # A tibble: 1 × 1
## frequency
## <int>
## 1 3# Coefficients
SummarySS <- summary(coda::mcmc(SSBoomNew$beta))
# Extract posterior means of beta
beta_post_mean <- SummarySS$statistics[, 1]
# Compare true vs estimated
plot(beta_true, beta_post_mean, pch = 19, col = "steelblue",
xlab = "True beta", ylab = "Posterior mean of beta",
main = "SSVS Shrinkage")
abline(0, 1, col = "red", lty = 2)
The examples and exercises presented thus far have considered scenarios in which the number of inputs is smaller than the number of observations (\(K < N\)). In Exercise 4, we challenge the Bayesian LASSO and SSVS in a setting where the number of inputs exceeds the sample size (\(K > N\)). As you will observe in that experiment, both the Bayesian LASSO and SSVS perform well. However, the Bayesian LASSO requires more time to produce results compared to SSVS in this exercise. Ročková and George (2018) propose a connection between the LASSO and spike-and-slab priors for variable selection in linear models, offering oracle properties and optimal posterior concentration even in high-dimensional settings where \(K > N\).
In addition, there are other Bayesian methods for regularization, such as the spike-and-slab approach proposed by Ishwaran and Rao (2005) and non-local priors introduced by Johnson and Rossell (2012), which can be implemented using the R packages spikeslab and mombf, respectively.