Bayesian Approaches to Shrinkage and Sparse Estimation

Bayesian Approaches to Shrinkage and Sparse Estimation large A guide for applied econometricians large

\setcounter{page}{1}

\onehalfspacing

Introduction

In all areas of human knowledge, datasets are increasing in both size and complexity, creating the need for richer models. This trend is also true for economic data, where high-dimensional and nonlinear/noparametric inference is the norm in several fields of applied econometric work. The purpose of this survey is to introduce the reader to Bayesian inference using shrinkage and variable selection priors. In particular we intend to demonstrate that the benefits of a Bayesian approach to high-dimensional estimation are manifold. Bayesian inference allows for a more accurate quantification of uncertainty. Parameters are treated as random variables that have their own probability density (or mass) functions. The use of a prior distribution provides a natural ground for enhancing possibly weak information in the likelihood.\footnote{Note that our interest here is in “wide” data (e.g.\ a linear regression model with more predictors than observations) where unrestricted estimation based only on the likelihood is either unreliable or impossible. In cases with “tall” data (many observations) the Bayesian posterior will tend to concentrate towards a point mass, i.e.\ uncertainty is small.} Our first aim is to explore in this review classes of priors that can recover popular penalized regression estimators, such as the lasso of Tibshirani1996. Next, we want to demonstrate how the Bayesian paradigm becomes a natural framework for combining prior forms in order to capture more complicated patterns of shrinkage and/or sparsity in the data. For example, RockovaGeorge2018 extend the lasso with ideas from the Bayesian variable selection literature in order to obtain a “spike and slab lasso” estimator that is empirically superior to shrinkage or variable selection alone, and has desirable theoretical guarantees. Finally, we aim to illustrate that the Bayesian framework is ideal for applied economists who want to use shrinkage or sparsity in more complex or unconventional settings. Economists might be interested in combining data-rigorous statistical variable selection with economic restrictions on certain parameters\footnote{For example, instead of the typical statistical shrinkage towards zero that indicates whether an effect is important or not, economists might want to shrink a parameter towards a calibrated value or a sign restriction provided by the solution of an economic model.}, or use a shrinkage estimator in a model with breaks, stochastic volatility, missing data or other complexities. Penalized and constrained maximum likelihood frameworks can deal with such cases, but computation is non-trivial because it relies on optimizing complex functions. We demonstrate emphatically in this survey paper that Bayesian computation provides numerous tools and algorithms for shrinkage and sparsity that can be incorporated in very complex statistical models with the same ease they are used in univariate linear regression settings.

Even though the notions of sparsity and shrinkage estimation are ubiquitous since the explosion of Big Data in all fields of science (e.g. we doubt there are many economists these days who haven't heard about the lasso), we want to clarify these terms before proceeding with our formal definitions. Sparsity refers to finding parameter estimates that have more zeros than not (where zeros in estimation means absence of some effect or relationship). Shrinkage means estimation where many parameter elements are suppressed towards zero, but they are not necessarily zero. While many readers might be familiar with these concepts, interpretation from a Bayesian point of view is slightly different from frequentist approaches. Sparsity is not identical for the simple reason that parameters in the Bayesian paradigm are (continuous, in many cases) random variables. Similarly, shrinkage estimation is embedded in Bayesian inference since any non-diffusing (non-flat) prior will tend to bias the likelihood; the frequentist statistician can only achieve shrinkage if they specify the estimation problem using an explicit penalized likelihood approach.

We explain these differences, and many more concepts, in this detailed review. We build our discussion gradually by introducing in this section basic components of Bayesian decision theory and estimation, and the principles of Bayesian model determination using the marginal likelihood. In Section 2 we introduce the concept of hierarchical priors and present the basic properties of a large class of hierarchical representations of Bayesian sparsity and shrinkage estimators. In Section 3 we focus on computation using hierarchical priors, and strategies for making inference in high-dimension computationally feasible. Section 4 demonstrates how the hierarchical priors and computational tools discussed in the previous sections, can be readily applied to a wide class of models that are important in economics and finance, as well as other fields of science. Section 5 concludes this review.

Throughout this review we make the assumption that the reader has a broad understanding of the concept of a prior distribution. If this is not the case, novice readers are advised to begin reading about the basics of Bayesian inference in (ref) and then move to (ref). More experienced readers, can move directly to (ref), skipping the material in this section.

Bayesian decision theory and estimation

In order to motivate shrinkage and sparsity, we first introduce the concept of loss-based estimation using a Bayesian decision theoretic approach. Detailed introductions can be found in Shrinkage2018 and Robert2007. Assume we have data $X \in \mathcal{X}$ where $\mathcal{X}$ (the sample space) is a measurable set of ${\rm I\!R}^{n}$, and parameters $\theta \in \Theta$ where $\Theta$ (the parameter space) is a measurable set of ${\rm I\!R}^{p}$. We define two probability density functions (p.d.f.) that are measurable on $\mathcal{X}$ and $\Theta$: a the likelihood function $p\left(X \vert \theta \right)$, and a prior function $\pi(\theta)$. Denote with $\widehat \theta(X)$ an estimator of $\theta$, that is, a measurable function of data $X$ that maps from ${\rm I\!R}^{n}$ to ${\rm I\!R}^{p}$.

Under these definitions we can now specify what is the loss and risk associated with the estimator $\widehat \theta(X)$. First, we can define loss functions of the form $L\left( \widehat \theta(X), \theta \right) = \rho\left( \widehat \theta(X), \theta \right)$ where $ \rho(\bullet)$ can be a symmetric loss function (the quadratic being the most popular) or any asymmetric loss function that measures how close $\widehat \theta(X)$ is to the true $\theta$. The Bayes risk associated with “decision” $\widehat \theta$ is defined as Shrinkage2018

The quantity $\mathcal{R} (\theta, \hat{\theta}) = E_{\theta} \left( L \left(\widehat \theta(X), \theta \right) \right)$ is the frequentist risk of $\widehat \theta$, which is defined as the expected value of the loss function over the data realization for a fixed $\theta$. In contrast, the Bayes risk in (ref) is the average of frequentist risk $\mathcal{R}$ with respect to the prior distribution $\pi(\theta)$. Frequentist decision theory aims at making the expected loss $\mathcal{R}(\theta, \hat{\theta})$ small, while Bayesian decision theory aims at finding the minimum of $r\left(\pi, \widehat \theta \right)$. In particular, the quantity

is the Bayes risk of the prior distribution $\pi$. Given a prior $\pi$, an associated Bayes estimator $\hat{\theta}_\pi$ is a minimizer in the sense that $r(\pi, \hat{\theta}_\pi) =r(\pi) $.

We can now define the concepts of minimaxity and admissibility. A decision rule (estimator) is admissible with respect to the loss function $L$ if and only if no other rule dominates it. That is, iff $r\left(\pi, \widetilde{\theta} \right) < r\left(\pi, \widehat \theta \right)$ then $\widetilde{\theta}$ is admissible. An estimator is $\hat{\theta}_0$ is minimax for a given loss function $L$ if

that is, it is the minimizer of the worst-case frequentist risk. For a given prior $\pi$, define an associated Bayes estimator $\hat{\theta}_\pi$. If $\sup_{\theta}\mathcal{R}(\theta, \hat{\theta}_\pi) = r(\pi,\hat{\theta}_\pi)$, then $\hat{\theta}_\pi$ can be shown to be minimax. In this case, the prior $\pi$ is least favorable in the sense that $r(\pi', \hat{\theta}_\pi) \leq r(\pi, \hat{\theta}_\pi)$ for all other priors $\pi'$. That is, $\hat{\theta}_\pi$ is the best with respect to the least favorable prior distribution $\pi(\theta)$. Minimaxity is a desirable feature for comparing estimators but, of course, it can still become a misleading measure of comparison; see a counterexample and further discussion in Robert2007. Finally, note that if a minimax estimator is a unique (Bayes) estimator, then this is also admissible.

Why is it important to think in terms of optimality of an estimator with respect to a loss function? To answer this question, consider the expected value of the squared error loss of a scalar, point estimator $\widehat \theta = \widehat \theta(X)$, which is also known as the mean squared error:

The first term in the last equation above is the variance of $\widehat \theta$, and the second term is the square of its bias. The least squares estimator, which in many simple linear settings coincides with the maximum likelihood estimator, has zero bias (unbiased) and is the “best” meaning that it has narrowest sampling distribution (minimum variance) among all unbiased estimators. Despite these two desirable properties, it is not necessarily the case that OLS will always have the lowest mean squared error. Indeed, in high-dimensional cases with fat data ($p$ large relative to $n$) the sample variance of the OLS will tend to become very large. In cases with more parameters than observations ($p>n$), the OLS estimator has infinite solutions and infinite variance. In such cases, there exist biased estimators that achieve much lower variance compared to the unbiased estimator, to the extend that this reduction in variance compensates for any increase in the square of the bias (making the total MSE of the biased estimator lower). Specifically in the case of out-of-sample prediction the MSE of our modeled variable will be larger if the estimation MSE in (ref) is high, showing that evaluating estimation loss might be more important than looking only at (minimum variance) unbiasedness.

A well-known illustration of this concept, that changed dramatically the way statisticians think about estimators, is the example of the James-Stein estimator. Assume our likelihood is $X \sim N_{p} \left(\theta, \underline{\sigma}^{2} I_{p} \right)$ where $\theta \in {\rm I\!R}^{p}$ is the unknown parameter and $\underline{\sigma}^{2}$ is assumed to be known. Stein1956 proved that the maximum likelihood estimator $\widehat \theta^{mle} = X$ is the minimum risk equivariant estimator under various loss functions, it is minimax, and it is admissible for $p=1,2$. However, for $p \geq 3$ the maximum likelihood estimator is inadmissible under a square loss function, and the James-Stein estimator

has lower risk than the MLE, that is, $\mathcal{R}(\widehat \theta^{JS}) < \mathcal{R}(\widehat \theta^{mle})$. EfronMorris1973 showed that the James-Stein estimator is a special case of an empirical Bayes estimator of $\theta$, that is, an estimator that places a Gaussian prior on $\theta$ and sets its prior variance to be a certain function of the data $X$. Stein's estimator minimizes the total quadratic risk of $\theta$, but there may be elements $\widehat{\theta}_{i}^{JS}$, $i \in [1,p]$, which have higher risk than the MLE. For that reason, EfronMorris1973 also propose a limited translation empirical Bayes estimator, which offers a compromise between Stein's estimator and the MLE.

Bayesian estimators are by default biased towards the prior expectation, which is a result of doing inference by using the information in both the likelihood and prior functions. Similarly, penalized likelihood estimators, such as the popular lasso of Tibshirani1996, constrain the likelihood function with a penalty that intends to introduce a similar bias. The purpose of this subsection is to introduce an alternative view to traditional econometric inference with small parameter space, where unbiasedness is the holy grail. In high-dimensional settings some estimation bias may be desirable, especially when the purpose is prediction in which case richly parameterized specifications are not welcome. In many instances, in-sample parameter estimation accuracy (instead of out-of-sample prediction) is of primary importance, for example, when the quantity of interest is an elasticity or a causal effect that can inform policy decisions. We show later in this survey that even in such cases Bayesian and frequentist penalized regression estimators can be desirable.

Principles of Bayesian Model Choice: A regression perspective

According to BDA2013 the process of Bayesian data analysis involves three steps

In the first part of this review, we use a simple linear regression setting as the basis for developing shrinkage and sparsity priors (step 1), for discussing posterior computation (step 2) and assessing model fit (step 3). By doing so we aim to offer the same level playing field for presenting various hierarchical prior formulations. The final section presents several extensions of shrinkage and sparsity priors in more complex settings, such as factor models, time-varying parameter regression, and cofounder selection in treatment effect estimation.

The regression model we build upon has the form

where $n$ is the number of observations, $y_{i}$ is a scalar dependent variable, $\bm X_{i}$ is a $1 \times p$ vector of covariates (or regressors or predictors) that can possibly include an intercept, dummies, exogenous variables or other effects (e.g.\ trend in a time-series setting), $\bm \beta$ is a $p \times 1$ vector of regression coefficients, and $\varepsilon_{i} \sim N(0,\sigma^2)$ is a Gaussian disturbance term with zero mean and scalar variance parameter $\sigma^2$. Within this setting our interest lies in obtaining “good” estimates of $\bm \beta$ and $\sigma^{2}$, specifically in settings with many covariates (“large $p$, small $n$” regression).

The linear regression formulation implies a certain Gaussian likelihood function $\mathcal{L}(\bm \beta,\sigma^{2} \vert \bm y, \bm X)$ that is proportional to the sampling density $p(\bm y \vert \bm \beta,\sigma^2)$. These two quantities are not identical because the likelihood is not a true density function.\footnote{The likelihood is a product of densities that lacks a normalizing constant. } The Bayesian needs to specify a joint prior distribution of the parameters, in the form $p(\bm \beta,\sigma^2)$. Bayes Theorem postulates that

but for the purpose of parameter estimation, in particular, it is easier to ignore $p(\bm y)$ since it is a normalizing constant (i.e.\ not a function of the parameters of interest $\bm \beta$, $\sigma^2$) and work instead with the formula

A default prior setting in Bayesian inference is the natural conjugate prior which is defined as

where $(\bm D, v_0, s_0)$ are prior hyperparameters chosen by the researcher. Due to the fact that the likelihood has a similar structure to this prior, it is trivial to prove (see the accompanying Technical Document) that the posterior is of the form {

}where $\bm V = \left( \bm X^{\prime} \bm X + \bm D^{-1} \right)^{-1}$, $v = v_0 + n+p$, $s^{2} = s_{0}^{2} +(\bm y - \bm X \bm \beta)'(\bm y - \bm X \bm \beta)+ \bm \beta' \bm D^{-1} \bm \beta$

$\bm X = \left[\bm X_{1}^{\prime},..., \bm X_{n}^{\prime}\right]^{\prime}$ and $\bm y = \left(y_{1},...,y_{n} \right)^{\prime}$.

Goodness of fit measures: Marginal likelihood and information criteria

While (ref) is required for the derivation of parameter posterior distributions, the quantity $p(\bm y)$ in (ref) is of paramount importance for Bayesian model determination. This is the prior predictive likelihood, more commonly known as the marginal likelihood, that is, the evidence in data $\bm y$ after we integrate out the effect of all possible values that the “random variables” $\bm \beta,\sigma^2$ can admit through their prior distribution. This can be proven via solving for $p(\bm y)$ in (ref): {

}where $\int_{-\infty}^{\infty} \int_{0}^{\infty} p(\bm \beta,\sigma^2 \vert y) d\bm \beta d\sigma^2= 1$ because this is a proper density. The marginal likelihood is the expected value of the likelihood where the expectation is taken with respect to the prior. Put differently, it is the prior mean of the likelihood function. An important characteristic of the marginal likelihood is that the integral in (ref) can only be calculated when the prior is a proper density, that is, if $p(\bm \beta,\sigma^2)$ integrates to one. The benchmark Uniform (Jeffrey's) prior on $\bm \beta$ and $\log(\sigma^2)$ is a key example where this condition fails and the marginal likelihood does not exist.

Assume we want to predict a new (future) observation $y_{n+1}$ given $\bm X_{n+1}$ using the prediction (out-of-sample) model $p(y_{n+1} \vert \bm \beta,\sigma^2, \bm y)$ which, in turn, is based on the in-sample estimated model $p(\bm y \vert \bm \beta,\sigma^2)$. We can then define the posterior predictive likelihood

which is the distribution of the out-of-sample data point marginalized over the posterior distribution of the model parameters.

Both quantities -- prior and posterior predictive distributions -- are fundamental for model assessment in Bayesian inference. In the benchmark case of the linear regression with the natural conjugate prior, the marginal likelihood can be derived analytically and is of the form

where $v_0,s_0,\bm D$ are parameters of the prior distribution (chosen by the researcher), and $v,s,\bm V$ are parameters of the posterior distribution whose values are provided in (ref) and $\bm \mu^* = \bm V (\bm X' \bm y)$.

The predictive likelihood is also available analytically and it is of the form

where we define the $p$-dimensional t-density with location $\bm \mu$, scale matrix $\bm \Sigma$, and degrees of freedom $d$ as

The marginal likelihood is rarely available analytically, and in most cases the integral in (ref) has to be approximated using Monte Carlo or numerical methods.\footnote{Two early examples are GelfandDey1994 and Chib1995; see also ChibJeliazkov2001 for a review.} In cases of either a complex model or a complex prior structure, or both, evaluating the marginal likelihood can become challenging, if not impossible. In such cases it might be easier to calculate the posterior predictive likelihood in (ref) using a procedure called leave one out cross-validation (LOO-CV). This would involve fitting the model in training data and then using a hold-out sample to evaluate the posterior predictive likelihood. Notice that if MCMC samples from the parameter posterior are available, evaluation of (ref) is straightforward using Monte Carlo integration.\footnote{Recognizing the numerical and computational shortcomings of model choice based on marginal likelihoods, there are several early studies that propose model choice criteria that are based on variants of the posterior predictive distribution, see Davison1986, GelfandGhosh1998, Gelmanetal1996, LaudIbrahim1995, IbrahimLaud1994 and SanMartiniSpezzaferri1984.}

When marginal or posterior predictive likelihoods are difficult to obtain, a (computationally) straightforward alternative strategy is to rely on information criteria. For example, the Bayesian information criterion (BIC), is a first-order approximation to the marginal likelihood. Performing a Taylor expansion around the posterior mode\footnote{The posterior mode is chosen such that the first derivative of the posterior is zero, which simplifies terms when taking the Taylor expansion; see Raftery1995 for a detailed proof.} $(\widetilde{\bm \beta},\widetilde{\sigma}^{2})$ for the logarithm of the term $p \left( \bm y \vert \bm \beta, \sigma^{2} \right) p \left(\bm \beta,\sigma^{2} \right)$ in (ref), we can write the log-marginal likelihood as

where $J_{n}\left(\widetilde{\bm \beta},\widetilde{\sigma}^{2} \right)$ is the expected Fisher information matrix of $p\left( \bm y \vert \bm \beta, \sigma^{2} \right) p \left(\bm \beta,\sigma^{2} \right)$ evaluated at the posterior mode $(\widetilde{\bm \beta},\widetilde{\sigma}^{2})$. In large samples, the posterior mode coincides with the MLE $(\widehat{\bm \beta},\widehat{\sigma}^{2})$. Considering this approximation and removing from (ref) any terms of order $O\left(1\right)$ or less, we obtain

The approximation above provides the basis for defining the Bayesian information criterion

where $\mathcal{L} \left( \widehat{\bm \beta}, \widehat{\sigma}^2 \vert \bm y, \bm X \right)$ is the likelihood function evaluated at the MLE.

The BIC is only a crude approximation to the marginal likelihood and it is based on a point estimate. An alternative popular criterion is the deviance information criterion (DIC) proposed by Spiegelhalteretal2002 which is of the form

The first term is the expectation of the data density with respect to the posterior\footnote{For that reason, the DIC is related to the posterior predictive likelihood, i.e.\ the integral in (ref), rather than the marginal likelihood.} which can be evaluated numerically from the MCMC output by taking the mean of $p(\bm \beta, \sigma^2 \vert \bm y)$ over all MCMC samples of the parameters. The second term is the value of the data density evaluated at the posterior mode $(\widetilde{\bm \beta},\widetilde{\sigma}^{2})$. For more information on the DIC see also ChanGrant2016, Spiegelhalteretal2014 and vanderLinde2005.

ChenChen2008 propose a modification to the Bayesian information criterion for high-dimensional spaces, which they call the extended Bayesian information criterion (EBIC). In the context of a proportional hazards model, VolinskyRaftery2000 propose a modification of the BIC penalty term that is consistent with a conjugate unit-information prior under this model. FosterGeorge1994 propose the risk inflation criterion (RIC) while FosterGeorge2000 present empirical Bayes selection criteria. Watanabe2010,Watanabe2013 derives the widely applicable information criterion (WAIC), also known as the Watanabe-Akaike information criterion since this criterion can be considered to be a Bayesian variant of the popular Akaike information criterion. Gelmanetal2014 and Vehtarietal2017 perform informative comparisons of the properties of BIC, DIC, WAIC and LOO-CV in a Bayesian context.

Testing hypotheses: Bayes factors

Consider now the case of two competing models, model one (denoted as $M_{1}$) and model two (denoted as $M_{2}$). For example, a key scenario that fits this setting, is that of testing hypotheses of the form $H_0: \beta_{j} = 0$ vs $H_1: \beta_{j} \neq 0$, for some $j=1,...,p$. Evidence in favor of either $H_{0}$ or $H_{1}$, corresponds to how good is the fit of two corresponding nested regression models ($M_{1}$ is unrestricted, and $M_{2}$ has the restriction $\beta_{j} = 0$ imposed). In this setting it is convenient to condition parameter posteriors and marginal likelihoods for each model on the random variable $M_{i}$, $i=1,2$, that indexes each of the two models. For example, $p(\bm \beta, \sigma^{2} \vert \bm y, M_{1})$ and $p(\bm y \vert M_{1})$ denote the parameter posterior and marginal likelihood, respectively, of regression model $1$. Consequently, the quantity

is the Bayes Factor between models $1$ and $2$. The quantity

is the posterior odds between models $1$ and $2$. It is defined as the product of the Bayes factor and the prior odds. If we assign equal model probabilities a-priori, then $p\left(M_1 \right)=p\left(M_2 \right)=\frac{1}{2}$ and the Bayes factor is identical to the posterior odds ratio. The Bayes factor above is a primary tool for assessing evidence in favor of a statistical model versus a competing model.

KassRaftery1995 provide a rule-of-thumb on how to interpret the statistical evidence against model $2$ based on ranges of values of $BF_{12}$: for values higher than three the evidence is substantial, for values higher than 10 it is strong, and for values higher than 100 it is decisive. Given that marginal likelihoods are not available with improper priors (even if the posterior is proper), there has been plenty of interest in calculating Bayes factors when such priors are used. Aitkin1991 proposes to calculate Bayes factors based on integrating the likelihood with the posterior -- this is equivalent to replacing $p(\bm \beta,\sigma^2)$ with $p(\bm \beta,\sigma^2 \vert \bm y)$ in (ref). This formulation allows to calculate “posterior” Bayes factors regardless of the prior structure of each model, and at the same time it avoids Lindley's paradox Aitkin1991. BergerPericchi1996,BergerPericchi1998 suggest the use of the intrinsic Bayes factor. Their suggestion involves splitting the data into $n$ subsets, such that one can obtain the marginal likelihood of the $i^{th}$ subset conditional on all other subsets. Subsequently, either the arithmetic or geometric average of the Bayes factors estimated in all $n$ subsets of the data can be used as the final estimate.

For nested model comparisons, VerdinelliWasserman1995 show that Bayes factors can be calculated using the Savage-Dickey density ratio (SDDR) approach. Consider two regression models as in (ref) but for notational simplicity set $p=1$, that is, only a single covariate is available. The first model, $M_{1}$, is an unrestricted model while model $M_{2}$ imposes the restriction $\beta = \beta^{\star}$ for some scalar value $\beta^{\star}$ (the previous example of testing of $H_0:\beta = 0$ vs $H_1:\beta \neq 0$ fits this setting). In this case the Bayes factor can be written as

that is, SSDR is the ratio of the marginal posterior and prior of $\beta$ under model $M_{2}$, evaluated at the point $\beta = \beta^{\star}$. In general it will be easy to evaluate these two distributions, especially when the Gibbs sampler is used for approximating the posterior distribution. This is because evaluation of the numerator using Monte Carlo integration would be fairly straightforward. Additionally, in the case of an independent prior of the form $p(\beta,\sigma^{2}) = p(\beta)p(\sigma^{2})$ the denominator above becomes $\int_{0}^{\infty} p\left( \beta^{\star},\sigma^{2} \vert M_{2} \right)d\sigma^2 = p\left( \beta^{\star} \vert M_{2} \right) \int_{0}^{\infty} p\left( \sigma^{2} \vert M_{2} \right)d\sigma^2 = p\left( \beta^{\star} \vert M_{2} \right)$, i.e.\ we only need to evaluate the (Gaussian) prior of $\beta$ at the point $\beta^{\star}$.

There are of course numerous other ways of obtaining approximations to the Bayes factors that do not explicitly involve calculating ratios of marginal likelihoods. GoutisRobert1998 propose an alternative procedure for testing nested models based on the Kullback-Leibler divergence. The idea is to compute the projection of the unrestricted model to the restricted parameter space, and use the corresponding minimum distance to judge whether or not the restricted model is appropriate. The same way we used the BIC to obtain a first-order approximation to the marginal likelihood, we can also use the BIC to obtain approximations to Bayes factors -- this approach is illustrated in Raftery1995. Notable early studies on the topic of Bayes factors include KassWasserman1995, DeSantisSpezzaferri1997, OHagan1995, BergerPericchi2001, BergerMortera1999, LewisRaftery1997, Raftery1996 and DiCiccioetal1997. A systematic review of methods for calculating Bayes factors can be found in KadaneLazar2004.

Finally, it is worth noting that in the case of nested hypothesis testing we can derive an optimal Bayesian point estimate by minimizing expected loss averaged over the two hypotheses, using posterior model probabilities as weights. That is, considering again the simple case with $p=1$ and ignoring the variance parameter $\sigma^2$ for simplicity, we aim to find point estimate $\widehat{\beta}$ such that the joint expected loss under the two models/hypotheses

achieves a minimum. Under a quadratic loss function $L \left( \beta,\widehat{\beta}\right)$, the posterior means are optimal meaning that the optimal estimator is

This estimator can be considered a Bayesian pre-test estimator, hence the acronym BPE in the equation above; see judgeetal1985 for a detailed discussion. In the next section we will generalize this result to the case of $K$ models, in order to motivate model choice in the presence of many models.

Model choice with many models: Bayesian model averaging

Model choice can have many forms, but the benchmark scenario that will motivate later in this paper to focus on shrinkage and sparse estimation, is that of model determination among many nested models. In particular, consider the problem of deciding which of $p$ variables in the covariate matrix $\bm X$ should be in the “optimal” regression model. Each covariate can have two outcomes, either it is included in a model or it is excluded, meaning that the model space in the presence of $p$ covariates is $2^{p}$. We denote the model set as $\mathcal{M}=\left\{ M_r: r=1,\ldots, 2^p \right\}$. The covariates that pertain to model $M_{r}$ are denoted in this subsection as $\bm X_{r}$ and their associated coefficients as $\bm \beta_{r}$. That is, $\bm X_{r}$ is a matrix that is constructed using only a subset of the columns in $\bm X$. Therefore, we denote regression model $M_r$ as\footnote{For simplicity we do not explicitly allow for an intercept. If an intercept is present in all competing models, then it is important to remove the sample mean from all covariates $\bm X$ (and, as a result, in all subsets $\bm X_{r}$) in order to ensure that the estimated intercept has exactly the same interpretation in all models. With demeaned covariates and the use of a flat prior, the intercept term becomes identical to the sample mean of $\bm y$ in all $2^{p}$ competing models.}

where $\bm X_{r}$ is $n \times p_{r}$ and $\bm \beta_{r}$ is $p_{r} \times 1$ with $p_{r} \in \{ 1,...,p\}$. Now with $2^p$ models, even for small $p$, pairwise model comparison based on Bayes factors is impractical and alternative computational methods are needed. Most importantly, in the presence of many models the researcher might not want to give the same weight to each and every model. For example, she might want to give more weight on parsimonious models or models that include a certain predictor suggested by some theory or common sense. For that reason we define prior model probabilities $p(M_r)$ with $\sum_{r=1}^{2^p} p(M_r) = 1$. Based on Bayes theorem, prior model probabilities combined with marginal likelihoods $p(\bm y \vert M_r)$ give posterior model probabilities

Bayesian model selection (BMS) corresponds to selecting the best model, that is, the model $M_r$ with the highest $p(M_r \vert \bm y)$. Bayesian model averaging (BMA) involves averaging over many models using $p(M_r \vert \bm y)$ as weights. That is, for a quantity of interest $\Delta$ (e.g.\ an out-of-sample observation $y_{n+1}$ of $\bm y$) BMA is constructed as the following weighted average

For small model spaces, typically when $p<30$ posterior model probabilities can be calculated analytically such that we can enumerate and estimate all $2^{p}$ available models. For $p>30$ it is impossible to enumerate and estimate all models in a deterministic way. In such cases, one can rely on Markov chain Monte Carlo algorithms which are able to “visit” in each iteration, in a stochastic way, the most probable models. Hoetingetal1999 and Fragosoetal2018 provide two systematic reviews on the topic.

While model selection and model averaging with an arbitrary number of models are straightforward extensions of the case with only two models, prior elicitation in multi-parameter and multi-model settings is anything but straightforward. In order to explain the intuition behind why this is the case, consider the natural conjugate prior defined previously, which in the case of model $M_r$ can be written as

Prior elicitation involves choice of $\bm D_{r}, v_0, s_0$. The hyperparameters $v_0, s_0$ are scalar in all regression models can be simply set to a small value close to zero, implying a Jeffrey's (diffuse) prior on $\sigma^2$. However, $\bm D_{r}$ is a matrix that changes size based on the number of predictors in model $M_{r}$. Assume for simplicity we define $\bm D_{r} = \tau \bm I_{p_{r}}$, with $\bm I_{p_{r}}$ the $p_{r} \times p_{r}$ identity matrix. In this case, prior elicitation breaks down to choosing a single hyperparameter $\tau$. We can't use the diffuse choice $\tau \rightarrow \infty$ because the marginal likelihood in (ref) will become infinite, hence, $\tau$ should be finite in the multi-model case. However, using the same finite value of $\tau$ in all models, doesn't mean that the effect of this prior is identical (that is, “objective”) for each model. Consider for instance two models, one with two predictors $\bm X_{2} = (\bm x_{1}, \bm x_{2})$ and a restricted model with only the first predictor $\bm X_{1} = \bm x_{1}$. The posterior variance is $\bm V_{r} = \sigma^{2} \left( \bm X_{r}^{\prime}\bm X_{r} + \left(\tau \bm I_{p_{r}}\right)^{-1} \right)^{-1}$ for each model $r=1,2$, so that the impact of $\tau$ on the common predictor in the two models will be identical only if $\bm x_{1}$ is not correlated with $\bm x_{2}$ and $\bm X_{2}^{\prime}\bm X_{2}$ becomes diagonal. If this is not the case, the correlation between the two predictors will imply that the effect of $\tau$ on the regression coefficient of $\bm x_{1}$ will not be the same in the two models. This issue complicates prior elicitation further when considering $p \gg 2$ correlated covariates, that also potentially have different units of measurement.\footnote{The scaling issue in $\bm X$ can be dealt with by standardizing the data, that is, dividing each column with its sample standard deviation. High correlation in columns of $\bm X$ can also be dealt with by orthogonalizing this matrix. While standardization is easy to apply and is recommended in all model averaging and variable selection algorithms, orthogonalization of the columns of $\bm X$ is only feasible when $n>p$. Therefore this latter procedure is not available in the high-dimensional case ($p>n$), which is exactly where there is higher chance of encountering many correlated predictors!}.

For that reason, many researchers have proposed empirical Bayes priors, in the spirit of the empirical Bayes formulation of Stein's estimation rule; see equation (ref) and discussion of EfronMorris1973. Empirical Bayes procedures allow to choose prior hyperparameters as a function of the data observations, sometimes also chosen to optimize some criterion (e.g.\ maximum marginal likelihood). A default prior for multi-model settings is the g-prior due to Zellner1986. The $g$-prior for model $M_{r}$ takes the form

where $\sigma^{2} \left(\bm X_{r}^{\prime} \bm X_{r} \right)^{-1}$ is essentially the covariance matrix associated with the OLS estimator $\widehat{\bm \beta}_{r}$ and $g$ a scalar tuning parameter. Under this prior, the posterior variance of $\bm \beta$ conditional on $\sigma^2$ becomes $\bm V_{r} = \frac{1}{1 + g} \times \sigma^{2}\left( \bm X_{r}^{\prime}\bm X_{r}\right)^{-1}$, such that the posterior variance is uniformly affected by selection of $g$. Consequently, the posterior mean/mode is

When $g \rightarrow 0$ the posterior mean tends to the OLS estimate of model $M_r$ ($\widehat{\bm \beta}_{r}$) while when $g \rightarrow \infty$ the posterior contracts towards zero. While the effect of the prior now depends in a straightforward, transparent way\footnote{We avoid using the term “objective”, first, because as GelmanHenning2017 argue it is counterproductive to do so and, second, because the $g$-prior is not in any way an objective prior.} on a single hyperparameter, choice of this hypeparameter is very important for determining marginal likelihoods and model probabilities.

Fernandezetal2001a,Fernandezetal2001b propose default values of $g$ in the context of Bayesian model averaging and Eicheretal2011 expand this discussion by considering further values of $g$. A benchmark suggestion of Fernandezetal2001b is to set $g \equiv g_{r} = p_{r}/n$, that is, a value of $g$ that is the ratio of the number of coefficients in each model $r$ over the total number of observations. Wide models with many covariates models will have larger $g$, thus, tending to shrink their posterior towards zero more aggressively. Put differently, the prior variance is getting smaller meaning that the information in the prior increases relative to the information in the likelihood. This is a basic principle of shrinkage and variable selection estimators: when $p$ is large and especially when $p>n$, the information in the likelihood is not sufficient to estimate all $p$ coefficients and the prior becomes increasingly important for determining posterior outcomes. That is, for both Bayesian and non-Bayesian approaches, the concepts of shrinkage and sparsity amount to the prior expectation that increasingly many coefficients a priori will be zero or close to zero.

Of course, there are more rigorous ways of selecting $g$. A key contribution is that of Liangetal2008 who put hyper-priors on $g$, treating it as a random variable. Such hierarchical approaches are the topic of close examination of the next section, so we won't expand on it here. Krishnaetal2009 extend the $g$-prior into an adaptive powered correlation prior of the form

where $\lambda \in \mathbb{R}$ controls the prior's response to collinearity in predictors. $\lambda=-1$ gives the original prior proposed by Arnold Zellner, while $\lambda=0$ gives the ridge regression prior.

While the $g$-prior addresses the issue of setting a prior on different regression models that might be nested and have correlated covariates, another important issue is how to define a prior on model space. For both conceptual and computational reasons Bayesians prefer to index all possible $2^p$ models using dummy variables $\bm \gamma = (\gamma_{1},...,\gamma_p)^{\prime}$. When $\gamma_j=0$ a covariate is excluded from a model and when $\gamma_{j}=1$ it is included. Therefore, the model with no predictors is indexed as $\bm \gamma = (0,...,0)^{\prime}$ and the model with all predictors is indexed as $\bm \gamma = (1,...,1)^{\prime}$. All intermediate models are indexed by vectors $\bm \gamma$ that are sequences of zeros and ones. Instead of placing priors on the model space, we can now explicitly consider priors on $\bm \gamma$, and the binomial distribution is a good candidate for a parameter that takes $0/1$ values. The binomial prior can become both uniform but also more informative when this is desirable (e.g.\ in high-dimensional spaces, where our prior is that only a small number of predictors will be important).

This setting that combines the $g$-prior on regression coefficients with a binomial prior on model space, is the major workhorse model for implementing Bayesian variable selection. While its theoretical underpinnings are well-understood (see Hoetingetal1999 for a thorough description), it provides the ground for some of the most interesting Bayesian work on computation in high-dimensional settings.\footnote{See for example, BottoloRichardson2010, Clydeetal2011, Dellaportasetal2002, Hansetal2007, JiSchmidler2013, MadiganYork1995, NottKohn2005 and Peltolaetal2012b.} At the same time this setting possesses implicitly the benefits of a hierarchical prior approach. Therefore, we use this brief discussion of BMA as a stepping stone for introducing in the next the concept of full-Bayes/hierarchical Bayes priors that result in shrinkage and sparse estimators.

Hierarchical (full Bayes) priors

When interest lies in models with many parameters, simple priors such as the benchmark natural conjugate prior presented in the previous section, are inadequate for learning interesting features about our parameters and for quantifying uncertainty. In statistics, the concept of hierarchical or multi-level modeling refers to the process of enhancing a simpler model with a richer specification that allows for learning interesting features of a multi-parameter vector, such as groupings or sparsity and shrinkage towards zero, where the latter being the main focus of this review. The Bayesian interpretation of hierarchical modeling involves specifying prior distributions for the prior hyperparameters of regression coefficients, especially when $p$ is large. A simple hierarchical specification for the regression coefficients $\bm \beta$,\footnote{Ignore estimation uncertainty of $\sigma^2$ for the moment, e.g.\ assume it is known and fixed.} takes the form

where $F(a,b)$ denotes some distribution function with hyperparameters $(a,b)$. Due to the fact that choice of $\tau^{2}$ is so crucial for the posterior outcome of $\beta_{j}$, the idea behind this hierarchical specification is to treat the hyperparameter $\tau^{2}$ as a random variable and learn about it from the data, via Bayes Theorem. For that reason, a prior such as the one in equations (ref) - (ref) is many times referred to as a full-Bayes prior, as it allows for full quantification of uncertainty around parameters of interest. While the example above pertains to linear regressions with Gaussian likelihood and prior distributions, Section 4 demonstrates that the concept of hierarchical priors is much more powerful and can be applied to numerous multivariate, non-Gaussian, nonlinear or other settings. Additionally, adaptive hierarchies can be defined in which $\beta_{j}$ depends on hyperparameters specific to this $j$-th element ($\tau^{2}_{j}$) that have their individual hyperprior distributions. Finally, if needed, further layers of the hierarchy can be defined: for instance, if choice of the hyperperameter $a$ of $\tau^{2}$ is not straightforward, we can define another level for the prior distribution of $a$, or we could introduce two variance parameters for $\beta_j$ in (ref).\footnote{For example, a powerful class of hierarchical priors called global-local shrinkage priors PolsonScott2010 provides an excellent benchmark for specifying appropriate hierarchical priors. Such priors are of the form

where $\tau^{2}$ is a global shrinkage parameter (applying the same shrinkage to the whole parameter vector $\bm \beta$) and $\lambda_{j}$ is a local shrinkage parameter (applying shrinkage only to $\beta_j$). As we see next, such priors will typically have at least three hierarchical layers, but in practical situations they tend to have many more (e.g.\ by putting priors on some or all of the hyperparameters $a, b, c, d$). }

An important feature of the hierarchical prior in equations (ref) - (ref) is that, while the conditional prior $p \left(\beta_{j} \vert \tau^{2} \right)$ is Gaussian, unconditionally the prior for $\beta_{j}$ is non-Gaussian. Indeed the marginal prior for $\beta_j$ becomes

that is, a scale mixture of normals representation that allows to approximate very complex prior shapes for $\beta_{j}$.\footnote{It is trivial to show that if $\tau^{2}$ is not a fixed parameter, then unconditionally the prior for $\beta_{j}$ always has excess kurtosis higher than zero, thus, being a leptokurtic distribution with tails thicker than the normal distribution.} Mixtures have the benefit of allowing for classification and grouping of parameters. In the case of identifying sparsity and shrinkage, we can think of the mixture prior as grouping parameters into “important” and “non-important”. Therefore, it is this implied mixture representation of hierarchical modeling with prior distributions that allows to extract interesting features in a multi-parameter setting. Finally, the posterior mode of $\bm \beta$ under a hierarchical prior has a penalized likelihood representation. For the linear regression model, penalized likelihood problems admit the following regularized least squares form

where the first term gives the solution to the usual least squares problem and the term $g\left( \bm \beta, \lambda \right)$ defines the penalty as a function of the regression parameters $\bm \beta$ and a scalar (or possibly vector) tuning parameter $\lambda$. Numerous penalized estimators, such as ridge (Tikhonov regularization), lasso, and elastic net fall under the general form in (ref), and Bayesian modal estimators under suitable hierarchical priors can fully recover all of them.

In order to understand the ability of hierarchical priors to classify parameters as important and non-important (or non-penalized and penalized), we plot in (ref) a normal prior with fixed variance vs three cases of a normal prior with variance parameter distributed as $\chi^{2}$ with one degree of freedom, exponential with rate parameter $\lambda=0.5$, and binomial with one trial and probability $\pi = 0.9$ (that is, a Bernoulli distribution). The simple normal prior provides more probability at the origin (zero) relative to its tails, however, it is fairly flat (diffuse) in a small area around zero. What the three mixture priors are introducing, is a more pronounced peak at zero such that when a parameter is in the region of zero it can be shrunk at a faster rate. At the same time, all three mixture distributions have fat tails, providing positive probability to parameter values that are far from zero. That is, these shapes allow for a clearer separation and classification of a parameter as being zero or non-zero. The extreme case of the Bernoulli prior on $\tau^{2}$ (bottom right panel of (ref)) creates a distribution that looks normal but also has a point mass at zero with high probability. Therefore, all three examples of hierarchical priors provide sharper inference in favor or against the groups of interest (important and non-important parameters).

Computation with hierarchical priors is reviewed in detail in the next section. For now it suffices to note that because of the conditional structure of hierarchical priors, conditional posteriors are typically easy to derive even if the joint parameter posterior is intractable. Sampling from these conditional posteriors using Markov chain Monte Carlo (the Gibbs sampler, in particular) is equivalent to taking samples from the intractable joint posterior. Additionally, several approximate methodologies such as variational Bayes and maximum a-posterior (MAP) estimation rely on similar conditional distributions. Therefore, in our discussion in this section we present various hierarchical priors, explain their properties and focus on deriving conditional posteriors. In the next section we discuss in more detail how to use these conditional posteriors to estimate the desired parameters.\footnote{Additional derivations and computational details can be found in the accompanying Technical Document.}

Diffusing hierarchical prior

A natural choice for the variance parameter $\tau^{2}$ in the hierarchical model of equations (ref) - (ref) is a prior distribution that is diffuse. Similar to Jeffrey's prior for the regression variance $\sigma^{2}$, the choice $\tau^{2} \sim U(0,\infty) $ equivalently $p(\tau^{2}) \propto \tau^{-2}$ can be thought as a default prior choice that reflects our lack of information about sparsity patterns in the data. We might want to also allow for each $\beta_j$ to be determined adaptively, in which case a Jeffrey's prior on hyperparameters $\tau_j^{2}$, $j=1,...,p$ can be defined. Therefore, the full hierarchical prior specification for the regression model is of the form

where $\bm D_{\tau} = diag(\tau_{1}^{2},...,\tau_{p}^{2})$. While a Jeffrey's prior on $\tau_j^2$ is a first natural attempt towards hierarchical prior modeling, as Lindley1983 notes, “a prior for $\tau^{2}$ that behaves like $\tau^{-2}$ will cause trouble" meaning it will lead to an improper posterior. gelman2006 examines this issue in more detail and explains why a $Uniform(-\infty,\infty)$ prior on $\log\left( \tau^{2}\right)$ would also not work. However, as KahnRaftery1992 and gelman2006 note, under certain conditions, Jeffrey's prior on $\tau^{2}$ yields a limiting proper posterior density. Note that the same improper density can be obtained from the prior $\tau^{2} \sim Inv-Gamma(\epsilon,\epsilon)$ for $\epsilon \rightarrow 0$ (see also (ref) below). gelman2006 argues that the $Inv-Gamma(\epsilon,\epsilon)$ prior does not have any proper limiting posterior distribution, such that inference becomes sensitive to the choice of $\epsilon$ -- simply setting $\epsilon$ to any “small” value is not a reliable solution.

Figueiredo2003 and BaeMallick2004 are examples of empirical studies that rely on shrinkage using a uniform hyperprior distribution. Tipping2001 specifies an inverse gamma prior on $\tau^{2}$ (and calls the resulting hierarchical structure a sparse Bayesian learning prior) and adopts the limiting case $\epsilon = 10^{-4}$ as the default hyperparmeter choice. Diffusing priors should not be the first choice in empirical settings especially in high-dimensional and ultra-high-dimensional settings. There are numerous other hyperprior distributions that are interpretable and have better theoretical guarantees gelman2006.

Student-t shrinkage

While we just argued that it is not desirable to use the inverse gamma distribution as a way of imposing a diffusing prior on $\tau^{2}$, informative inverse gamma priors provide flexible parametric shrinkage. Following the specification of the normal-inverse gamma prior in ArmaganZaretzki2010, we write this prior using the following form

where $\bm D_{\tau} = diag(\tau_{1}^{2},...,\tau_{p}^{2})$. This is a scale mixture of normals representation of the fat-tailed and leptokurtic Student-t distribution. Similar to our arguments in (ref) the excess kurtosis of the Student-t results in shrinkage towards zero at a faster rate than the simple normal distribution. At the same time the fatter tails accommodate values of $\tau^2$ that can be far from zero. In (ref) we illustrate the shape of the marginal distribution of $\beta_j$ for various values of the parameters $\rho,\xi$.

Similar to an inverse gamma prior for the variance parameter $\sigma^2$, the conjugacy of this distribution allows for numerous methods of inference using this prior. For example, Tipping2001 uses type-II maximum likelihood methods Berger1985, but (as we discuss in the following section) variational Bayes and other approximate algorithms are also trivial to derive. ArmaganZaretzki2010 show that conditional posteriors are of the form

where $\bm V = \left( \sigma^{-2} \bm X^{\prime} \bm X + \bm D_{\tau}^{-1} \right)^{-1}$ and $\Psi = (\bm y - \bm X \bm \beta)^{\prime}(\bm y - \bm X \bm \beta)$. The Gibbs sampler can be used to sample sequentially from these conditional posteriors, as these samples are guaranteed to be samples from the desired joint parameter posterior.

For the conditional posterior (ref) of the prior precisions, $\frac{1}{\tau_{j}^{2}} $, we have

where the proportional sign is with respect to $\frac{1}{\tau_{j}^{2}}$'s. It can be seen that the conditional posterior of $\frac{1}{\tau_{j}^{2}}$'s is independent across $j$ and that it has the form in (ref).

Normal-gamma priors

CaronDoucet2008 proposed the normal-gamma family of hierarchical priors, and GriffinBrown2010,GriffinBrown2017 established further results and their excellent properties. This prior takes the following hierarchical form

where again $\bm D_{\tau} = diag(\tau_{1}^{2},...,\tau_{p}^{2})$. The pdf of $\tau_j$ is

such that the marginal pdf of $\beta_j$ is

where $\mathcal{K}_v$ is the modified Bessel function of the second kind, and the tails of this distribution decrease in $\vert \beta_j \vert^{(\lambda - 1)}\exp\left( \gamma \vert \beta_j \vert \right)$.

Due to the connection of the gamma distribution with a wide array of other distributions (e.g.\ inverse gamma, inverse Gaussian, $\chi^{2}$, etc) choice of the hyperparameters $\lambda$ and $\gamma^{2}$ can result in various shapes for the unconditional distribution of $\bm \beta$ that have different shrinkage properties. This prior becomes diffusing when $\lambda,\gamma^2 \rightarrow 0$, however, this choice falls under the same critique of gelman2006 for the diffusing inverse gamma prior. This is due to the fact that when $\lambda < 1/2$ the normal-gamma prior places infinite mass in the vicinity of zero, that is, $\lim_{\beta_j \rightarrow 0} p(\beta_j) = \infty$.

LASSO prior and extensions

The least absolute shrinkage and selection operator (lasso) of Tibshirani1996 has been established as a key workhorse of scientists in all fields working with high-dimensional settings. The estimator takes the form

where $t$ is a prespecified free parameter that determines the degree of regularization. The Lagrangian form of this program is

where $\vert\vert x \vert\vert_{1} = \sum \vert x_{i} \vert $ is the $\mathcal{l}_{1}$ norm and $\vert\vert x \vert\vert_{2} = \sqrt{\sum x_{i} ^{2} }$ is the $\mathcal{l}_{2}$ norm. $\lambda$ is a tuning parameter related to $t$, controlling for how strongly shrinkage is exercised. As $\lambda \rightarrow 0$ the penalty term vanishes and the lasso becomes indistinguishable from the least squares problem. This optimization formula is related to basis pursuit denoising, which is the preferred term for the lasso among researchers in computer science and signal processing.

Tibshirani1996 first noted that the lasso estimate can be derived as a Bayes posterior mode under the following Laplace prior distribution

However, as castillo2015 note the full posterior distribution under a Laplace prior does not contract at the same rate as its mode, making uncertainty quantification using the Bayesian lasso unreliable. The intuition behind this is that the $\lambda$ coefficient above needs to be large enough to penalize coefficients $\beta_j$ to zero, but not too large such that nonzero coefficients can be modeled. This issue is addressed by modifications such as the adaptive lasso (AlhamzawiAli2018; see end of this section) and the spike and slab lasso (RockovaGeorge2018; see section on spike and slab priors) and is related to the motivating arguments of JohnsonRossell2010 for proposing the non-local priors (see relevant section below).

The first application of the lasso prior stems from computing science and is due to Girolami2001. While the joint parameter posterior under a Laplace prior is not of standard form, Girolami2001 used variational Bayes inference (which at the time was not popular in mainstream statistics) to approximate the posterior mean and variance. Figueiredo2003 used the fact that the Laplace prior admits a hierarchical representation in the form of a normal-exponential (double exponential) mixture. The hierarchical representation of this prior is of the form

where the exponential distribution has the functional form $p(\tau^{2} \vert \lambda^2) = \left( \frac{\lambda^{2}}{2} \right) \exp\left( \frac{\lambda^{2}}{2} \tau_{j}^{2} \right)$. The marginal distribution for $\beta$ conditional on $\lambda^2$ is of the form

which is the desired Laplace distribution for $\beta_j$. Figueiredo2003 derived an EM algorithm for obtaining the posterior mode (MAP estimator).

A formal Bayesian treatment of the Bayesian lasso using MCMC can be found in ParkCasella2008. These authors choose to specify the Bayesian lasso as a normal-exponential mixture but conditional on the regression variance $\sigma^2$. This is because a hierarchical prior on $\beta_{j}$ that is independent of $\sigma^{2}$ results in a multimodal posterior for $\beta_{j}$. The ParkCasella2008 Laplace prior takes the form

where $\bm D_{\tau} = diag(\tau_{1}^{2},...,\tau_{p}^{2})$. Conditional posteriors under this hierarchical representation are trivial to derive and more details can be found in the accompanying Technical Document.

The approach in ParkCasella2008 is probably the most widely used but it is not the only one available. Hans2009 specified the lasso in terms of the normal orthant distribution. Let $\mathcal{Z} = \lbrace -1, 1 \rbrace^{p}$ represent the set of all $2^p$ possible vectors of length $p$ whose elements are $\pm 1$. For any realization $z \in \mathcal{Z}$ define the orthant $\mathcal{O}_{z} \subset {\rm I\!R}^p$. If $\bm \beta \in \mathcal{O}_{z}$, then $\beta_{j} \geq 0$ if $z=1$ and $\beta_{j} <0$ if $z=-1$. Then $\bm \beta$ follows the normal-orthant distribution with mean $m$ and covariance $S$, which is of the form

The Hans2009 prior takes the form

and, using the definition of the normal orthant distribution, conditional posteriors are of the form

where:

The conditional posterior of $\sigma^2$ is not of a standard form and, therefore, cannot be sampled directly. Hans2009 suggests a simple accept/reject step within the Gibbs sampler that allows to obtain approximate samples from the posterior of $\sigma^2$. Finally, MallickYi2014 propose a third hierarchical representation of the Laplace prior, this time as a mixture of Uniform distributions (see our Technical Document for details of this algorithm).

There are numerous extensions to the basic lasso that come in various forms. For example, the elastic net combines the benefits of ridge regression ($\mathcal{l}_{2}$ penalization) and the lasso ($\mathcal{l}_{1}$ penalization) by solving the problem

where now $\lambda_{1}$ and $\lambda_{2}$ are tuning parameters. The Bayesian prior that provides the solution to the elastic net estimation problem is of the form

LiLin2010 start from this prior and derive a mixture approximation and a Gibbs sampler that has the minor disadvantage that requires an accept-reject step for obtaining samples from the conditional posterior of $\sigma^{2}$ (similar to the sampler of Hans2009 for the lasso). The formulation of the elastic net prior in Kyungetal2010 is slightly different to the one above, but they manage to derive a slightly different mixture representation and a slightly more straightforward Gibbs sampler.

Other popular extensions to the lasso include the group lasso that allows for group shrinkage; the fused lasso that allows for spatial or temporal relationships between neighbouring parameters; and the adaptive lasso that fixes some variable selection consistency issues with the regular lasso. All these extensions have straightforward hierarchical forms, and we refer the reader to discussions in Kyungetal2010, GriffinBrown2011, Lengetal2014 and AlhamzawiAli2018, among several other studies. Our Technical document provides details of posterior inference using the elastic net, group lasso, fused lasso and adaptive lasso.

Generalized double Pareto shrinkage

Armaganetal2013a propose the following generalized double Pareto (GDP) prior on $\bm \beta$

This distribution can be represented using the familiar, from the Bayesian lasso, normal-exponential-gamma mixture. The only difference is that, while the Exponential component has the same rate parameter for all $j=1,...,p$, in the representation of the GDP mixture this parameter is adaptive. The generalized double Pareto distribution has a spike at zero with Student’s t-like heavy tails.

The generalized double Pareto prior takes the form

where $\bm D_{\tau} = diag(\tau_{1}^{2},...,\tau_{p}^{2})$.

The conditional posteriors are of the form

where $\bm V = \left( \bm X^{\prime} \bm X + \bm D_{\tau}^{-1} \right)^{-1}$, $\bm D_{\tau}^{-1} = diag(\tau_{1}^{-2},...,\tau_{p}^{-2})$ and $\Psi = (\bm y - \bm X \bm \beta)^{\prime}(\bm y - \bm X \bm \beta)$. Paletal2017 show, both theoretically and numerically, that the above “three-block” Gibbs sampler is less efficient than a modified two-block Gibbs sampler they propose.

Dirichlet-Laplace

The Dirichlet-Laplace prior was introduced in Bhattacharyaetal2015, and ZhangBondell2018 studied its posterior consistency as well as consistency in variable selection in the context of a linear regression model. The Dirichlet-Laplace hierarchical prior, which is a generalization of the Laplace prior, takes the form

where $\bm D_{\lambda,\tau,\psi} = diag(\lambda^{2} \tau_{1}^{2}\psi_{1}^{2},...,\lambda^{2} \tau_{p}^{2}\psi_{p}^{2})$.

The conditional posteriors are of the form

where $a^*=(n+p)/2$, $b^*=( \Psi +\bm \beta^{\prime} \bm D_{\tau,\lambda,\psi}^{-1} \bm \beta )/2 $, $c^*=\sqrt{\lambda^{2} \psi_{j}^{2} \sigma^2 / \beta_{j}^{2}}$, $\bm V = \left( \bm X^{\prime} \bm X + \bm D_{\tau,\lambda,\psi}^{-1} \right)^{-1}$, and $\Psi = (\bm y - \bm X \bm \beta)^{\prime}(\bm y - \bm X \bm \beta)$. $IG$ is the two-parameter inverse Gaussian distribution, and $GIG$ is the three-parameter generalized inverse Gaussian.

Horseshoe prior

The horseshoe prior was first introduced by Carvalhoetal2010 and it since its inception has been the most popular and influential hierarchical prior in Bayesian inference. The survey paper by Bhadraetal2020 provides a thorough review of the applications of this prior in numerous inference problems in statistics and machine learning, including nonlinear models and neural networks. The Horseshoe is a prime representative of the class of global-local shrinkage priors (see Footnote (ref)) and it can be represented as a scale mixture of normals with half-Cauchy mixing distributions. That is, the prior has the following formulation

where $\bm \Lambda = diag(\lambda_{1}^{2},...,\lambda_{p}^{2})$, and $C^{+}(0,\alpha)$ is the half-Cauchy distribution on the positive reals with scale parameter $\alpha$. That is, $\lambda_{j}$ has conditional prior density

Under this hierarchical specification, the marginal prior for each $\beta_{j}$ is unbounded at the origin and has tails that decay polynomially.

There are numerous theoretical results established for this prior, most notably DattaGhosh2013 and vanderPasetal2014, and the reader is referred to Bhadraetal2020 for a more detailed discussion. There are also various computational approaches to the Horseshoe (see the accompanying Technical Document for details), but the most straightforward is the one proposed by MakalicSchmidt2016. These authors note that the half-Cauchy distribution can be written as a mixture of inverse-gamma distributions. In particular, if

then $x \sim C^{+}(0,\alpha)$. Therefore, the MakalicSchmidt2016 prior takes the form

where $\bm \Lambda = diag(\lambda_{1}^{2},...,\lambda_{p}^{2})$.

The conditional posteriors are of the form

where $\bm V = \left( \bm X^{\prime} \bm X + \bm D_{\tau,\lambda}^{-1} \right)^{-1}$, $\bm D_{\tau,\lambda}= diag(\tau^{2}\lambda_{1}^{2},...,\tau^{2}\lambda_{p}^{2}) = \tau^{2} \bm \Lambda$ and $\Psi = (\bm y - \bm X \bm \beta)^{\prime}(\bm y - \bm X \bm \beta)$.

Generalized Beta mixtures of Gaussians

Armaganetal2011 motivate the use of a three-parameter beta (TPB) distribution for the prior variance parameter, as a flexible class of shrinkage priors. The TPB distribution takes the form

for $0<x<1$, $a,b,\varphi>0$. The TPB normal scale mixture representation for the distribution of random variable $\beta_{j}$ is given by

Proposition 1 in Armaganetal2011 shows that this distribution can either be written as normal-inverted beta mixture, or a normal-gamma-gamma mixture. The second choice gives a very straightforward Gibbs sampler scheme, and it can be seen as a special case of the normal-gamma class of priors GriffinBrown2017.

The Generalized Beta mixtures of Gaussians prior takes the form

where $\bm D_{\tau} = diag(\tau_{1}^{2},...,\tau_{p}^{2})$. Note that setting $a=b=1/2$ we can obtain the horseshoe prior of Carvalhoetal2010. For other choices we can recover popular cases of shrinkage priors.

The conditional posteriors are of the form

where $\bm V = \left( \sigma^{-2} \bm X^{\prime} \bm X + \bm D_{\tau}^{-1} \right)^{-1}$, $\bm D_{\tau} = diag(\tau_{1}^{2},...,\tau_{p}^{2})$ and $\Psi = (\bm y - \bm X \bm \beta)^{\prime}(\bm y - \bm X \bm \beta)$.

The TPB normal mixture includes as special cases Strawderman-Berger and horseshoe priors.

Non-local priors

Non-local priors have been proposed by JohnsonRossell2010 in the context of hypothesis testing of the form $H_{0}: \beta_j = 0$ vs $H_{1}: \beta_j \neq 0$. From a frequentist perspective, such testing procedures are used in order to find out how likely it would be for a set of observations to occur under the null hypothesis. However, in a Bayesian setting the data are assumed to be observed once, and parameters are continuous random variables. Traditional (local) priors put significant probability in both the null and alternative hypotheses, thus, making it harder for the (continuous) posterior distribution to detect-non zero coefficients asymptotically. Non-local densities place zero probability at zero, and this feature allows such priors to separate more clearly between the null and alternative hypotheses. That is, such priors do not place any prior probability under the null.\footnote{As JohnsonRossell2010 note:

}

Any distribution that “decreases to 0 near the boundaries between disjoint null and alternative parameter spaces might be considered” JohnsonRossell2010 to be a non-local prior density. Within the context of a linear regression setting similar to the one defined in (ref), JohnsonRossell2012 propose two specific classes of priors. The first class of prior densities for $\bm \beta$ consists of product moment (pMOM) densities, which are defined as

(ref) plots the pMOM density for $\tau^{2}=\sigma^{2}=1$ and for three values of $r$ ($r=1,2,3$). This graph clearly shows the shapes that this prior can achieve, especially with regards to the rate at which this prior decreases in the region of zero. The second class of prior densities consists of the product inverse moment (piMOM) densities, which are defined as

In both of these two priors, $\tau^{2}$ is a scale parameter that determines dispersion of the prior around zero. Therefore, this parameter determines the size of the regression coefficients that will be shrunk to zero, and it is of prime importance. JohnsonRossell2012 and Shinetal2018 treat $\tau^2$ to be fixed and show that high-dimensional model selection consistency is achieved under the pMOM prior, as long as $\tau^{2}$ is of a larger order than $\log p$ and it increases subexponentially in $n$. However, fixing this parameter might not be desirable in most applied high-dimensional problems\footnote{For example, JohnsonRossell2012 note that if the covariate matrix $\bm X$ is not standardized, then it would be important to define an adaptive shrinkage parameter $\tau^2_j$ for each $j=1,...,p$. In such a case, choice of each individual $\tau^2_j$ for large $p$ becomes inconvenient, if not infeasible.}, and a hierarchical approach might be desirable. Caoetal2020 propose a hyperprior density for $\tau^{2}$ of the form

The hierarchical pMOM (or “hyper-pMOM”) prior they propose achieves strong model selection consistency when $p$ increases at a polynomial rate with $n$. Unfortunately, neither the pMOM, hyper-pMOM or piMOM priors allows for a closed-form computation of joint, marginal or conditional posteriors. Therefore, Caoetal2020 rely on Laplace approximations.

Spike and slab priors

Similar to non-local priors, spike and slab priors allow for variable selection and testing of the hypotheses $H_{0}: \beta_j = 0$ vs $H_{1}: \beta_j \neq 0$. Unlike non-local priors, spike and slab prior densities place significant probability into both hypotheses. In a regression context, the spike and slab prior MitchellBeauchamp1988 takes the form

for each $j=1,...,p$, where $\delta_{0}(\beta_j)$ is the Dirac delta function placing point mass at zero and $\gamma_j$ are 0/1 (dummy) variables indicating whether column $j$ of $\bm X$ is included in the regression or not. The mechanism with which this prior classifies predictors as important or not, is simple: when $\gamma_j=1$ the prior for $\beta_j$ is $N(0,\tau^2)$, that is, estimation is not restricted by the prior for reasonably large values of $\tau^{2}$; when $\gamma_j=0$ the prior becomes a point mass function concentrated at zero and it dominates the likelihood such that the posterior is also concentrates its mass at zero. The concept of variable selection is fully determined by the indicator random variables $\gamma_j$'s. Samples from the posterior of each $\gamma_j$ will be sequences of zeros and ones, and the posterior mean denotes the posterior inclusion probability of each predictor in the best model. For example, if we sample MCMC 10,000 draws and find that 2,000 times $\gamma_j=1$, then the posterior mean is simply $2000/10000=0.2$ which translates into $20\%$ posterior inclusion probability of predictor $j$. BarbieriBerger2004 show that the median probability model, that is, the model where only variables with probabilities larger than 0.5 are selected/retained, is optimal for prediction. OHaraSilanpaa2009 suggest that as a variable selection mechanism such variable selection priors should work well up to cases where $p$ is 10-15 times larger than $n$, but of course this proportion is only a rule of thumb that is heavily determined by the informativeness of the data and modeling choices.

The spike and slab prior belongs to the general class of hierarchical full-Bayes priors introduced earlier in this section, since it can be written in the form

If, in addition, we introduce a hyperprior distribution on $\tau^2$ (e.g.\ inverse-gamma, see IshwaranRao2003), then the spike and slab prior is not only a hierarchical prior, but also belongs to the class of local-global shrinkage priors with global shrinkage parameter $\tau^2$ and local shrinkage parameters $\gamma_j$. In signal processing and similar fields, the spike and slab is known as a “normal-Bernoulli” or “Gaussian-Bernoulli” prior.

A third parametric formulation of this particular spike and slab prior is due to KuoMallick1998. In their formulation the regression model with variable selection prior is written as

where $\beta_{j}$ is the coefficient on predictor $j$ and $\gamma_j$ is a 0/1 variable indicating whether predictor $j$ is included in the model. This formulation is equivalent to the previous two, but it implies that the vector of indicators $\bm \gamma$ enters only via the likelihood and not through the (hierarchical) prior for $\bm \beta$. In the KuoMallick1998 formulation each $\beta_{j}$ will simply have a typical Gaussian prior with variance $\tau^2$. Notice that when $\gamma_j=1$, $\beta_j$ will be sampled from its posterior, but when $\gamma_j=0$, $\beta_j$ is not identified. In this case what happens is -- as is the case with any unidentified parameter in a Bayesian setting (e.g.\ mutlicollinearity) -- that $\beta_j$ is sampled from its prior. This lack of identification of $\beta_j$ is not a problem, as what we care about is the joint effect $\gamma_j \times \beta_j$ and the fact that predictor $j$ simply has to be removed whenever $\gamma_j=0$. This detail means that in variable selection a-la KuoMallick1998 the posterior of $ \beta_j$ with $\gamma_j=0$ will be equal to its normal prior, while the posterior of the same parameter under the spike and slab prior of equation (ref) is a point mass at zero. Other than this (possibly minor) difference, Bayesian variable selection using all three forms presented above is conceptually and empirically comparable.

The class of spike and slab priors and its theoretical properties have been studied extensively in the literature; see JohnstoneSilverman2004, IshwaranRao2005b, Jiang2006, Bogdanetal2011 and castillo2012. From an applied scientist's point of view, the spike and slab prior is very versatile and can take numerous useful forms.\footnote{For example, KoopKorobilis2016 specify a spike and slab prior that is able to search for homogeneities in panel data. That is, the spike and slab prior is modified in order to test the hypothesis of the form $H_{0}: \beta_i = \beta_j$ vs $H_{1}: \beta_i \neq \beta_j$.} We next briefly review possible formulations of the spike and slab prior, and their implications for modeling coefficients and selecting variables in a linear regression. We finish this section with a discussion of some key computational aspects of this class of priors.

Tuning of parameters in the spike and slab prior

In the formulation in (ref) one only has to choose the variance parameter $\tau^2$. This cannot be zero because the slab will become identical to the spike component, and it cannot become infinity because it would also be impossible to separate the spike from the slab component (remember from the previous section that Bayes factors with diffuse priors do not exist). Therefore, $\tau^2$ has to be quite different from zero and not too large (e.g.\ $\tau^2=4$ is a reasonable choice). Of course one can use any of the hyperprior distributions already explored in the previous sections, e.g.\ the choice $\tau^2 \sim exponential(\lambda^{2}/2)$ will convert the slab into a Laplace prior. However, one should be careful not to overshrink the slab (e.g.\ by setting $\lambda$ too large in the Laplace prior) because then the spike and slab will be indistinguishable and posterior inclusion probabilities will be meaningless.

A computationally more efficient formulation of the spike and slab prior (at least within an MCMC setting) is the one proposed by GeorgeMcCulloch1993,GeorgeMcCulloch1997, where both the spike and slab distributions are continuous

where $\tau_0^{2}$ is a “small” variance parameter (corresponding to the spike) and $\tau_1^{2}$ is a “large” variance parameter (corresponding to the slab). In the limit, when $\tau_{0}^{2} = 0$, the spike becomes the Dirac delta at zero, but for any other values of $\tau_{0}^{2}$ close but different from zero the spike distribution is unable to shrink $\beta_j$ exactly to zero. That is, this version of the spike and slab is appropriate for testing $H_{0}: \beta_j \approx 0$ vs $H_{1}: \beta_j \neq 0$, that is, it provides a soft thresholding rule. Chipmanetal2001 provide the threshold value above (below) which a regression coefficient is classified as belonging to the slab (spike) component and is not shrunk (shrunk) to zero:

Therefore, elicitation of $\tau_{0}^{2},\tau_{1}^{2}$ becomes very important for variable selection in the GeorgeMcCulloch1993 prior. NarisettyHe2014 show that fixing these two variance hyperparameters may result in variable selection inconsistency, and propose values that are functions of $n$ and $p$ that ensure good performance of the prior when the data dimensions increase. IshwaranRao2005b set $\tau_{0}^2 = \tau^{2}$ and $\tau_{1}^2 = c \tau^{2}$ where $c>>1$ and $\tau^{2} \sim Inv-Gamma$, although $\tau^2$ could also follow any of the hierarchical distributions defined previously, e.g.\ Horseshoe or Laplace. SylviaHelga2010 go one step further by motivating a mix-and-match strategy where $\tau_{0}^{2}$ has a Laplace prior, while $\tau_{1}^{2}$ has an inverse-gamma prior. More recently, RockovaGeorge2018 showed that, under mild conditions, a spike and slab lasso prior produces posterior distributions that concentrate asymptotically around the true regression coefficients at nearly the minimax rate. In their formulation both the spike and the slab are based on Laplace distributions (represented as normal-exponential mixtures), with the spike distribution shrunk more aggressively than the slab distribution.

An important feature of variable selection priors is the prior on $\gamma_j$. As in (ref) this is typically Bernoulli with prior probability $\pi_0$, or equivalently a binomial prior for the full vector $\bm \gamma = \left( \gamma_1,...,\gamma_p \right)^{\prime}$. Unfortunately, the choice $\pi_0=0.5$ in a binomial prior is not uniform as it implies a prior expectation that half of the $p$ predictors will be included in the final model. Therefore, in high-dimensional settings it is customary to set this parameter to a value closer to zero, e.g.\ $\pi_0=0.1$. If desired, a prior can be placed on this parameter and a conjugate choice is the beta distribution, that is, $\pi_0 \sim Beta(1,\alpha_0)$. The choice $\alpha_0=1$ makes this prior uniform, but in high-dimensional cases it will be preferable to set $\alpha_0$ to become proportional to the number of predictors $p$. Note that in the presence of a beta hyperprior on $\pi_0$, it is not necessary to use indicator variables $\gamma_j$. For example, following Dunsonetal2008 we can specify a spike and slab of the form\footnote{See also Korobilis2013a,Korobilis2013b,Korobilis2016 for related priors applied to econometric contexts such as dynamic regressions and vector autoregresions.}

that provides a smoother mixture of the two components. (We can, of course, specify an equivalent formulation for the GeorgeMcCulloch1993 continuous spike and slab formulation.) Finally, Carvalhoetal2008 turn this latter formulation into a sparsity inducing variable selection prior by replacing (ref) with

that is, a spike and slab prior for $\pi_0$. Finally, YuanLin2005 propose a prior for $ \bm \gamma$ that accounts for correlation in predictors, such that if two predictors are highly correlated only one is included in the selected model. In their formulation they multiply the standard binomial prior for $\bm \gamma$ with the determinant of the Gram matrix of predictors, that is, $\vert \bm X^{\prime} \bm X \vert$. High-correlated predictors have small $\vert \bm X^{\prime} \bm X \vert$ and are discouraged from being selected. Such enhancements of the base spike and slab prior are important for variable selection, because marginal inclusion probabilities may be poor under high correlation. In particular, highly correlated predictors may be jointly selected often but each predictor only a small number of times.

Computation with spike and slab priors

Computation with spike and slab priors is as straightforward as is the case with most other hierarchical priors. Conditional on $\gamma_j$ being either zero or one, the prior for $\beta_j$ is either a point mass at zero or normal MitchellBeauchamp1988 or it is one of two normal components GeorgeMcCulloch1993. Therefore, conditional on $\gamma_j$, results for the normal linear model can be used. The same holds in the case where the components of the spike and slab are non-normal, rather they are Student-t, Laplace etc: as long a hierarchical prior structure is used and the prior can be written in conditionally normal form, derivation of conditional posteriors is straightforward.

Regarding posterior computation of $\gamma_j$'s this usually has to be done element-by-element, that is, we need to derive $\gamma_j$ conditional on $\bm \gamma_{-j}$ (the set $\bm \gamma$ with the $j$-th element removed).\footnote{For that reason, when the Gibbs sampler is used to sample from the conditional posterior of $\gamma_j$ given $\bm \gamma_{-j}$, it is advisable in each Gibbs iteration to sample in random order $j$ to avoid high autocorrelation of samples.} However, in the case of the MitchellBeauchamp1988 prior of (ref), the conditional posterior $p(\gamma_j \vert \bm \gamma_{-j}, \bm \beta, \sigma^{2}, \bm y)$ cannot be used to obtain samples from the posterior of $\gamma_j$. Intuitively, this is because when we sample $\gamma_j = 0$ then the prior for $\beta_j$ is the Dirac delta function that puts infinite mass at zero. Therefore, in the next iteration $p(\gamma_j \vert \bm \gamma_{-j}, \bm \beta, \sigma^{2}, data)$ will give $\gamma_j=0$ with probability one, meaning that the sampler will get stuck in a loop where the only possible outcome is $\beta_j=\gamma_j=0$. This is not an issue in the continuous spike and slab prior of GeorgeMcCulloch1993, since the spike is a continuous normal distribution and allows samples of $\beta_j$ to be slightly different from zero.

To see this, let's derive $p(\gamma_j \vert \bm \gamma_{-j}, \bm \beta, \sigma^{2}, \bm y)$ in the case of the spike and slab prior of equations (ref) - (ref), which we rewrite for convenience

For simplicity, we do not introduce prior distributions on $\tau^{2}$ and $\pi_0$, so we assume these are fixed and chosen by the researcher. Using Bayes theorem, the posterior of $\gamma_j=0$ is

In this decomposition, the first term is provided by the likelihood where we set the $j$-th element of $\bm \beta$ equal to zero (since $\gamma_j=0$), regardless of what the sampled value $\beta_j$ is in the previous iteration of the Gibbs sampler. This is a normal distribution with mean $\bm X \bm \beta^{\star}$, where $\bm \beta^{\star}$ is equal to $\beta$ with the $j$-th element equal to zero, and variance $\sigma^{2}$. The second term is the prior for $\beta_j$ under the restriction $\gamma_j=0$, that is, the Dirac delta density $\delta_{0}(\beta_j)$. The last term is given simply by the Bernoulli prior for $\gamma_j$ and it is equal to $(1-\pi_0)$. Therefore, this posterior is:

Using similar arguments, we have that

Therefore, the conditional posterior of $\gamma_j$ is {\tiny

}Notice how the Dirac delta $\delta_{0}(\beta_j)$ enters the denominator term. If in the sampling process it happens to sample $\gamma_j=0$, then $\beta_j=0$ and any subsequent $\gamma_j$'s will also be zero for ever. This is because once a $\beta_j=0$ is observed, $\delta_{0}(\beta_j)$ becomes infinite and the ratio in the Bernoulli posterior is zero.

The solution to this problem is integration. That is, we need to remove dependence to $\beta_j$, and instead of the posterior $p(\gamma_j \vert \bm \gamma_{-j}, \bm \beta, \sigma^{2}, \bm y)$ we compute $p(\gamma_j \vert \bm \gamma_{-j}, \bm \beta_{-j}, \sigma^{2}, \bm y)$, that is, we integrate out $\beta_j$ and condition only on $\bm \beta_{-j}$. Intuitively, because $\gamma_j$ depends only to $\beta_j$ through the spike and slab prior (i.e.\ it is independent to $\bm \beta_{-j}$), the ratio in the Bernoulli posterior of (ref) will only involve the densities $p( \bm y \vert \gamma_j=0,\bm \gamma_{-j},\bm \beta, \sigma^{2})$, $p( \bm y \vert \gamma_j=1,\bm \gamma_{-j},\bm \beta, \sigma^{2})$ and $p(\gamma_j = 0)$, $p(\gamma_j = 1)$. The accompanying Technical document provides details of conditional posteriors under various forms of spike and slab prior distributions, including cases with more complex hiearchical layers such as the spike and slab lasso of RockovaGeorge2018.

Monte Carlo study: Specification of spike and slab priors for variable selection

Consider a GeorgeMcCulloch1993,GeorgeMcCulloch1997 type spike and slab prior

The conditional posteriors of $\bm \beta , \sigma^2, \bm \gamma$, and $ \pi_0$ are

where $\phi(\cdot \vert m, v)$ is the normal density with mean $m$ and variance $v$ and $\bm D$ is a diagonal matrix with diagonal elements $\{ (1-\gamma_j)\tau_{0j}^2 + \gamma_j \tau_{1j}^2 \}_{j=1}^p$.

SSVS-Lasso

Suppose we employ a Laplace density for the slab component

and consider three different ways of defining priors for the spike component that are commonly used in practice, which we define as SSVS-Lasso 1-3.

In SSVS-Lasso-1, $\tau^2_{0j}$ is fixed i.e.\ $\tau^2_{0j}=c_1$ for some small $c_1>0$ and in SSVS-Lasso-2, it is proportional to the prior variance for the slab component i.e.\ $\tau^2_{0j} =c_2 \tau^2_{1j}$ for some small $c_2>0$. In both SSVS-Lasso-1 and 2, with the prior $\lambda_1^{2} \sim Gamma(r_1,\delta_1) $, the prior variance for the slab is updated according to

In SSVS-Lasso-3, we place two separate Laplace densities on the components i.e.\

with $\lambda_0 \gg \lambda_1$ so that the density for $N(0,\sigma^2 \tau^2_{0j})$ is the “spike” and $N(0,\sigma^2 \tau^2_{1j})$ is the “slab”. This is similar to the spike-and-slab Lasso in RockovaGeorge2014 and BaiRockovaGeorge2021 \footnote{They propose an EM algorithm for estimation.}. The prior variances are updated according to

Which specifications are appropriate in applications? In order to investigate this question, we consider simulation by generating data from the regression model $\bm y = \bm X \bm \beta + \bm \epsilon$ with $\bm \epsilon \sim N_n(\bm 0,\sigma^2 \bm I_n)$. We let $n=100$ and $\sigma^2=3$. We construct the true vector of slope parameters $\bm \beta =c \tilde{\bm \beta}$ by assigning values $\{1.5, -1.5, 2, -2, 2.5, -2.5 \}$ to the first 6 elements of $ \tilde{\bm \beta}$ and setting others to zero. We choose a constant $c>0$ to achieve a desired level of signal-to-noise ratio \footnote{ In a general linear regression $\bm y = \bm X \bm \beta+ \bm \epsilon$, the signal-to-noise ratio (SNR) is defined as $SNR=\frac{|| \bm \Sigma_X^{1/2} \bm \beta ||^2 }{ \sigma^2 }$ where $\sigma^2$ is the error variance and $\bm \Sigma_X$ is a $p \times p$ covariance matrix of $\bm X$. $||\bm \Sigma_X^{1/2} \bm \beta||^2 = \bm \beta' \bm \Sigma_X\bm \beta$ measures the overall signal strength. A related quantity is $R^2_{pop}$, the population value of $R^2$, defined as $\frac{SNR}{1+SNR}$. }. The data matrix $\bm X$ is generated from the multivariate normal distribution with mean zero and covariance matrix being an identity matrix. The covariates are standardized for estimation. We examine different values of the number of covariates $p\in \{50, 100,300 \}$ and the signal-to-noise ratio $R^2_{pop}\in \{ 0.4,0.8 \}$.

The analysis was repeated 100 times with new covariates and responses generated each time. For each, the metrics recorded were: the bias and MSE of the first 6 elements of the coefficients vectors, the number of false negatives (FN), the number of false positives (FP), and the number of true positives (TP). Posterior means were used as point estimates of the slope coefficients and the error variance. We utilize the post-processing approach of LiPati2017 in order to categorize the covariates into signals and noises.

We compare performance of SSVS-Lasso 1-3 under different data generating processes. We fix $c=d=1$ so that the prior on the inclusion probability $\theta$ is uniform. We also let $a=b=0.1$. For SSVS-Lasso 1 and 2, the hyperparameters for the prior on $\lambda^2_1$ are fixed as $r_1=1$ and $\delta_1=1$. We let $c_1=10^{-4}$ for SSVS-Lasso-1, $c_2=10^{-4}$ for SSVS-Lasso-2, and $\lambda_0=20, \lambda_1=1$ for SSVS-Lasso-3. We also show results of NarisettyHe2014 which is a two component mixture of normals prior with fixed prior variances and KuoMallick1998 which can be seen as a spike-and-slab prior with the spike component being a point mass at zero. (ref) summarizes results.

Under a relatively strong signal i.e.\ $R^2_{pop}=0.8$, generally speaking, SSVS-Lasso-3 and Narisetty-He outperform others in all measures, and this tendency becomes more apparent in high-dimensional case i.e.\ $p=300$. Kuo-Mallick tends to have larger FPs than others. When the signal-to-ratio is lower, i.e.\ $R^2_{pop}=0.4$, SSVS-Lasso-3 outperforms others in terms of bias/MSE of the signals and shows reasonable performance in other metrics.

\FloatBarrier

Bayesian Computation with hierarchical priors

We have established that hierarchical priors obey conditional structures that make derivation of conditional posteriors a straightforward business. As a consequence, the Gibbs sampler is the primary computational tool for variable selection and shrinkage problems using hierarchical priors. However, exactly because such full Bayes shrinkage estimators are mostly needed in high and ultra-high dimensions, the Gibbs sampler and related Monte Carlo-based methods become computational costly. In such cases there are numerous other strategies that allow for faster computation. These strategies include approximate methods for computing marginal posterior distributions, or iterative, non-sampling methods that approximate the posterior mode or mean. Many of these algorithms originate in computing science, where data dimensions have always been larger than traditional economic data sets. Currently, Bayesian computation in high-dimensional spaces -- especially in the presence of hierarchical priors -- is the topic of an expanding research agenda in mainstream statistics as well as in the field of machine learning. In this section, we summarize this vibrant research, focusing on both MCMC and fast approximate algorithms.

Brute-force/analytical algorithms

Analytical algorithms for hierarchical priors, in general, do not exist apart from a few special cases that can be fairly restrictive. In the context of estimating a normal mean $\theta$ (see our discussion of EfronMorris1973 in (ref)), KahnRaftery1992 put uniform hyperpriors on the mean and variance hyperparameters of a normal prior distribution of $\theta$. In order to obtain the posteriors of these hyperparameters they need to integrate $\theta$, something they are able to do numerically since in their case only univariate integrals are involved on the support $[0,1]$. In the context of a regression with spike and slab prior, Clyde1999 shows that if the design is orthogonal (that is if the Gram matrix $\bm X^{\prime} \bm X \propto I$) and the regression variance $\sigma^2$ is known, variable selection indicators $\bm \gamma$ can be obtained without resorting to Monte Carlo methods (either Gibbs sampler or Monte Carlo). PapaspiliopoulosRossell2017 also derive an efficient non-sampling algorithm for Bayesian model averaging in regressions with a block diagonal design.\footnote{Examples of modeling scenarios with block-diagonal matrix $\bm X^{\prime} \bm X$ include time-series regressions with time-varying parameters, and vector autoregressions written in “seemingly unrelated regressions” form; see (ref) for more details.} Their methods involve calculation of model probabilities and parameter estimates using one-dimensional numerical integration.

An interesting case that allows for approximately analytical posterior inference using hierarchical priors is provided in vandenBoometal2015a and vandenBoometal2015b. These authors use rotation matrices to partition the regression model into a component explained by predictor $X_{j}$ and all remaining predictors $\bm X_{(-j)}$. In particular, they split the regression into two components

These authors use a non-shrinking natural conjugate prior in the first part of the rotated regression in order to obtain analytically an estimate of $\bm \beta_{(-j)}$ and $\sigma^2$. Then conditional on these estimates, they introduce in the second part a hierarchical shrinkage prior on $\beta_j$ and derive analytically its posterior, since the regression variance is known using its estimate from the first partition. This procedure requires to repeat the rotation and partition of the regression model for each predictor $j$, $j=1,...,p$. The outcome is an analytical derivation of the posterior of each element $\beta_j$ of $\bm \beta$ under a hierarchical prior that would otherwise require posterior simulation. This algorithm is of course approximate because it requires to obtain the posterior of $\beta_j$ by integrating out the influence of each $\bm \beta_{(-j)}$ using a natural conjugate prior, rather than the same hierarchical prior used on $\beta_j$. KorobilisPettenuzzo2019 extend this idea to various several hierarchical priors, including normal-Jeffrey's, spike and slab, and normal-gamma.

Gibbs sampler

We have already established that complex distributions (e.g.\ Student-t, Laplace, normal-half Cauchy) can be written in a conditionally conjugate form by using hierarchical representations. Depending on whether we also condition on the regression variance parameter, or not, we obtain the following two hierarchical prior formulations

where $\bm D_{\tau} = diag(\tau_{1}^2,...,\tau_{p}^2)$ and depending on the structure of the distribution $p(\bm \tau^2)$ (which itself can be a hierarchical mixture of several distributions) we obtained the various interesting cases we explored so far.

Because of the conditional structure of the prior, posterior conditionals are easy to derive. For example, consider the case of the independent prior, then the joint posterior is of the form

We can derive the conditional posterior of $\bm \beta$ as

because $p(\bm \tau^{2})$ and $p(\sigma^{2})$ are constants when conditioning on $\bm \tau^{2} , \sigma^{2}$ (because they do not involve the random variable $\bm \beta$). The prior distribution $p(\bm \beta \vert \tau^{2})$ is normal and, due to the modeling assumptions, $p \left(\bm y \vert \bm \beta, \sigma^{2} \right)$ is also normal. As a result, the conditional posterior for $p \left( \bm \beta \vert \bm \tau^{2} , \sigma^{2}, \bm y \right)$ is identical to the conditional posterior under the non-hierarchical version of the same prior BDA2013. Similarly, the conditional posterior for $\sigma^{2}$ becomes

which is also identical to the case of the non-hierarchical independent normal-inverse gamma prior. Finally, the conditional posterior for $\tau$ becomes

The data density $p \left(\bm y \vert \bm \beta, \sigma^{2} \right)$ does not contain information about $\bm \tau$ so it becomes a constant. Instead the “model” for $\bm \tau^2$ is provided by the density $p(\bm \beta \vert \bm \tau^{2})$ where $\bm \beta$ are observed “data” (fixed to their sampled values), since in this conditional posterior the only random variable is $\bm \tau^2$.

It becomes apparent that because of the conditional structure of hierarchical priors, for the vast majority of hierarchical priors we have common formulas for the conditional posteriors of $\bm \beta$ and $\sigma^{2}$, while the formulation of the conditional posterior of $\bm \tau^{2}$ will depend on how complicate its prior is. Under the natural conjugate prior the conditional posteriors are of the form

where $\bm V = \left( \bm X^{\prime} \bm X + \bm D_{\tau}^{-1} \right)^{-1}$ and $\bm \Psi=\left( \bm y - \bm X \bm \beta\right)'\left( \bm y - \bm X \bm \beta\right)$. Under the independent prior the posteriors are of the form

where $\bm V = \left( \bm X^{\prime} \bm X/\sigma^{2} + \bm D_{\tau}^{-1} \right)^{-1}$ and $\bm \Psi=\left( \bm y - \bm X \bm \beta\right)'\left( \bm y - \bm X \bm \beta\right)$. Notice how $\sigma^2$ affects the posterior of $\bm \tau^{2}$ in the conjugate prior case, while $\sigma^{2}$ doesn't show up in the posterior of $\bm \tau^{2}$. See the Technical Document for derivations.

A Gibbs sampler will cycle through equations (ref) or equations (ref), obtaining a sample of parameters conditional on all others, until a large enough sample from the posterior of each parameter is available. Results in PalKhare2014 and KhareHobert2013 establish, for a large class of hyper-prior distributions $p(\bm \tau^{2})$, that the above Gibbs sampler is ergodic and has the joint posterior $p\left(\bm \beta,\bm \tau^{2}, \sigma^{2} \vert \bm y \right)$ as its stationary density. Despite its ergodicity, the basic Gibbs sampler for models with hierarchical priors may suffer from slow mixing and convergence to the desired posterior. The conditional structure of a hierarchical prior implies a long chain of dependence of $\bm \beta$ on $\bm \tau^2$ and its hyper-priors. For example, the representation of the Horseshoe prior suggested by MakalicSchmidt2016 as a hierarchical mixture of a normal distribution and four inverse gamma hyper-prior distributions (see (ref)) is one example where slow mixing might become a serious issue. For that reason, in the case of the Horseshoe in particular, several authors propose to use more efficient slice sampling schemes, some of which we explore in detail in the accompanying Technical Document. Another disadvantage of the Gibbs sampler is the fact that sampling becomes cumbersome as the dimension $p$ of covariates increases. In the next we explore various approaches for speeding up MCMC and for dealing with convergence issues, particularly in the case where $p$ is large or even $p>>n$.\footnote{All the approaches we explore propose novel ways of sampling from the parameter posterior under the Gibbs sampling scheme. However, given a specific algorithm, the ability of the programming language to handle large matrices is also important. This is illustrated, for example, in Matusevichetal2016, where in the context of Bayesian variable selection they combine array database management systems (DBMS) and R processing capabilities, allowing R to handle large matrices without running out of RAM.}

Fast sampling from Normal posteriors

In high-dimensional settings with $p$ large, the most cumbersome step in a Gibbs sampler for the linear regression model with hierarchical priors is sampling from the $p$-variate normal conditional posterior distribution of $\bm \beta$. This step involves an inversion of precision matrix $\bm Q = \bm V^{-1} = \left( \bm X^{\prime} \bm X + \bm D_{\tau}^{-1} \right)$ (in the case of the natural conjugate prior) or $\bm Q = \bm V^{-1} = \left( \bm X^{\prime} \bm X/ \sigma^{2} + \bm D_{\tau}^{-1} \right)$ (in the case of the independent prior) in order to obtain the posterior covariance matrix $\bm V$.\footnote{Note that in the case of the natural conjugate prior without hierarchical structure, the matrix $\bm D_{\tau}$ is known (calibrated by the researcher), as is the data information $\bm X^{\prime} \bm X$. In this case, one could calculate and invert $\bm Q$ once, outside the loop of the Gibbs sampler. However, the presence of a hierarchical prior for $\bm \tau^{2}$ means that the matrix $\bm D_{\tau}$ changes values in each Gibbs iteration. Therefore, $\bm V^{-1}$ should be computed and inverted in each iteration (regardless of whether we used the independent prior on $\bm \beta$ or not).} Next, the Cholesky decomposition of $\bm V$ is needed in order to sample from the desired normal distribution. While the inversion step can be sped up (e.g.\ by using Woodbury's identity), standard built-in algorithms (in various programming languages) for obtaining the Cholesky decomposition of a $p \times p$ matrix have a worst asymptotic complexity (measured in flops) of $\mathcal{O}(p^{3})$. Therefore, simple sampling from a normal posterior is deemed to become computationally cumbersome, if not infeasible, as $p$ increases.

Rue2001 provides a precision-based sampler in order to obtain samples from a normal distribution efficiently when the precision matrix $\bm Q = \left( \bm X^{\prime} \bm X + \bm D_{\tau}^{-1} \right)$ is known. Due to the fact that the Bayesian conditional posterior of $\bm \beta$ (ignoring $\sigma^{2}$, e.g.\ assume it is fixed to value 1) is of the form $\bm \beta \vert \bullet \sim N( \bm V \bm X^{\prime} \bm y, \bm V)$, the procedure proposed by Rue2001 takes the following form

It is trivial to show that $E(\bm \beta) =\bm \mu= (\bm L^{\prime -1}\bm L^{-1}) \bm X^{\prime} \bm y= ( \bm L\bm L^{\prime} )^{-1}\bm X^{\prime} \bm y = \bm V \bm X^{\prime} \bm y$ and $cov(\bm \beta) = cov(\bm \mu + \bm u) = \bm L^{\prime -1} cov(\bm Z)\bm L^{-1} = \bm L^{\prime -1} \bm L^{-1} = \bm V $, which means that the above procedure provides valid samples from the desired normal distribution. The main feature of this algorithm is that it requires to invert the Cholesky factor of $\bm Q$, instead of inverting $\bm Q$ itself to obtain $\bm V$. While the worst case complexity of this algorithm is also $\mathcal{O}(p^3)$, it provides high efficiency gains in certain classes of models, e.g.\ when the Gram matrix $\bm X^{\prime} \bm X$ is block-diagonal (assuming the prior variance $\bm D_{\tau}$ is diagonal or at most block-diagonal of similar structure).

More recently, Bhattacharyaetal2016 proposed an efficient algorithm that makes full use of Woodbury matrix inversion lemma in the context of generating normal variates from a distribution of the form $\bm \beta \vert \bullet \sim N( \bm V \bm X^{\prime} \bm y, \bm V)$ (again for simplicity, ignore $\sigma^{2}$). Their algorithm takes the following form

It is also easy to show that the sample of $\bm \beta$ comes from the desired normal distribution with mean $\bm V \bm X^{\prime} \bm y$ and variance $\bm V$. Note that this algorithm requires inversion of the $n \times n$ matrix $( \bm X \bm D_{\tau} \bm X^{\prime} + \bm I_{n} )^{-1}$, while sampling directly from the normal posterior requires inversion of the $p \times p$ matrix $\bm Q$. Therefore, step 3 in this algorithm only becomes efficient for $p>>n$. The main efficiency gains in this algorithm stem from the fact that it requires to sample from two normal distributions with diagonal covariance matrices (a $p$-variate distribution with covariance $\bm D_{\tau}$ and an $n$-variate distribution with identity covariance). Generating uncorrelated normal draws is much more efficient than sampling directly from the $p$-variate normal posterior of $\bm \beta$ using the full covariance matrix $\bm V$. In particular, the worst-case complexity (asymptotic upper bound) of this algorithm is $\mathcal{O}(n^{2}p)$, that is, it is only linear in $p$. Therefore, for $n>p$ this algorithm will perform worse than the algorithm of Rue2001 since the term $n^{2}$ will dominate, but this algorithm shines in the $p > n$ case where it can offer some dramatic improvements in computation times.

Note that for very large $p$, computation of $\bm X \bm D_{\tau} \bm X^{\prime}$ in step 3.\ above will become cumbersome. In such ultra high-dimensional cases, Johndrowetal2020 provide an approximate version of Bhattacharyaetal2016. This involves removing “irrelevant” columns of $\bm X$ such that the above product can be computed using a significantly smaller number of algorithmic operations.

Scalable Gibbs

The standard form of the Gibbs sampler for the linear regression model with hierarchical prior contains three blocks as in (ref). There is one block for each set of parameters, namely $\bm \beta$, $\bm \tau^{2}$ and $\sigma^{2}$. In the context of hierarchical priors, Rajaratnam2019 propose to sample $(\bm \beta, \sigma^{2})$ in one block. Their proposed Gibbs sampler is more efficient, as reducing the number of blocks to sample from, also reduces correlation among draws from parameter posteriors. When working with the natural conjugate form of a hierarchical prior, the scalable Gibbs algorithm requires only to sample from $p((\bm \beta, \sigma^{2}) \vert \bm \tau^{2}, \bm y)$ and $p(\bm \tau^{2} \vert (\bm \beta, \sigma^{2}), \bm y)$. The first joint distribution for $(\bm \beta, \sigma^{2}) $ can be approximated by first sampling from $p(\sigma^{2} \vert \bm \tau^{2}, \bm y)$ (notice the lack of dependence on $\bm \beta$) and subsequently from $p( \bm \beta \vert \sigma^{2}, \bm \tau^{2} , \bm y)$. The scalable Gibbs algorithm has the following form

where $\bm V$ has the same definition as in (ref). The proof why the posterior for $\sigma^{2}$, after integrating out $\bm \beta$, has the form shown above, can be found in the Appendix of Rajaratnam2019; see also Paletal2017.

Skinny Gibbs

In the context of the spike and slab prior with continuous spike and slab distributions, Narisettyetal2018 propose an efficient sampling scheme that also separates the posterior into a mixture distribution with independent components. We remind that the continuous spike and slab prior GeorgeMcCulloch1993,NarisettyHe2014 can be written in matrix form as

where $\bm \Gamma = diag(\bm \gamma)$, or more compactly

where $\bm D_{\gamma} = diag\left( (1-\gamma_{1})^{2} \tau_{0}^{2} + \gamma_{1}^{2} \tau_{1}^{2},..., (1-\gamma_{p})^{2} \tau_{0}^{2} + \gamma_{p}^{2} \tau_{1}^{2} \right) $. The conditional posterior under this prior is of the form

where $\bm V = \left( \bm X^{\prime} \bm X/\sigma^2 + \bm D_{\gamma} ^{-1} \right)^{-1}$ and $\vert \bullet$ denotes conditioning on other parameters in the model as well as the data. Computing $\bm V$ requires an inversion as well as obtaining the Cholesky decomposition in order to sample from the $p$-dimensional normal posterior. As we already saw, when $p$ is large these operations can be extremely cumbersome.

The skinny Gibbs algorithm of Narisettyetal2018 solves this issue by splitting the conditional posterior into two independent components, an active (A) and an inactive (I) and sample $\bm \beta$ as the union of the following conditionals

where $\bm \beta_{A}$ is the $p_{A}$-dimensional vector of elements of $\bm \beta$ corresponding to $\gamma_{j}=1$, and $\bm \beta_{I}$ are the remaining $p_{I} = p - p_{A}$ elements that correspond to $\gamma_{j}=0$. In the first posterior the covariance matrix is $\bm V_{A} = \left( \bm X^{\prime} \bm X/\sigma^2 + \frac{1}{\tau_{1}^{2}}\bm I \right)^{-1} $, while in the second posterior $\bm V_{I} = \left(n + \frac{1}{\tau_{0}^{2}} \right)^{-1} \bm I_{p_I} $. In sparse settings we would expect to find that $p_{I} >> p_{A}$, meaning that the bulk of the elements of $\bm \beta$ would be sampled as restricted elements $\bm \beta_{I}$. This means that we can sample very efficiently $p_{I}$ coefficients from a normal posterior with diagonal convariance matrix, and sample the remaining $p_{A}$ elements from a normal posterior with a full covariance matrix.

Of course, in practical situations it is expected that the two sets of coefficients, $\bm \beta_{A}$ and $\bm \beta_{I}$, will be correlated with each other. As a result, sampling the full vector $\bm \beta$ from two independent conditional posteriors could leave us with a significant approximation error. For that reason, Narisettyetal2018 add a “compensation term” in the conditional posterior of $\bm \gamma$ that accounts for the approximation involved in sampling the coefficients $\bm \beta$. This term ensures that the skinny Gibbs converges to a stationary distribution, while keeping computational complexity minimal. Narisettyetal2018 show that the skinny Gibbs posterior possesses strong selection consistency property.

Orthogonal Data Augmentation

In the conext of variable selection using spike and slab priors, GhoshClyde2011 note that when the design matrix $\bm X$ is orthogonal, a stochastic search using MCMC can become very efficient when $p$ is large. Similarly for penalized likelihood problems (whether Bayesian or not), orthogonal designs can be very efficient and result in consistent variable selection regardless of how large $p$ is relative to $n$. The proposal of GhoshClyde2011 is to augment the $n \times p$ correlated design matrix $\bm X$ with an $n_{a} \times p$ matrix $\bm X_{a}$, such that the $(n+n_{a}) \times p $ “complete” design matrix

has orthogonal columns, that is, $\bm X_{c}^{\prime} \bm X_{c} = \bm X^{\prime} \bm X + \bm X_{a}^{\prime} \bm X_{a} = \bm W$, where $\bm W = diag(w_{1},...,w_{p})$ is a diagonal matrix with $w_{j}>0$. Furthermore, we have the restriction that the augmented data matrix $\bm X_{a}$ has real entries, and $\bm X_{a}^{\prime} \bm X_{a}$ must be a positive semidefinite symmetric matrix. GhoshClyde2011 select a diagonal matrix $\bm W$ which then implies the value of $\bm X_{a}$ from the orthogonality condition $\bm A \equiv \bm X_{a}^{\prime} \bm X_{a} = \bm W - \bm X^{\prime} \bm X $. $\bm X_{a}$ can be obtained as the symmetric matrix square root of $\bm A$, thus, ensuring that its entries are real. GhoshClyde2011 discuss in detail ways of choosing $\bm W$; for example, since in most variable selection settings the columns of $\bm X$ are typically standardized to have unit norm or variance, one can set $w_{1}=...=w_{p} =w$, such that choice of $\bm W$ collapses into a choice of a scalar $w$.

Once $\bm X_{a}$ has been specified, we can estimate the augmented orthogonal regression model

where $y_{c} = \left[ \bm y^{\prime} , \bm y_{a}^{\prime} \right]^{\prime}$ with $\bm y_{a}$ latent data which we need to sample and $\bm \varepsilon_{c} \sim N_{(n+n_a)}(\bm 0, \sigma^{2} \bm I)$. In the most general case, GhoshClyde2011 consider a hierarchical variable selection prior, which combines a spike at zero with a component that is Student-t, obtained via normal-inverse-gamma mixture (see (ref))

for all $j=1,...,p$. For $\alpha=1$ this prior becomes a heavy-tailed Cauchy distribution of the form $\beta_j \sim C(0,\sigma^{2} \gamma_{j})$. GhoshClyde2011 propose a Gibbs sampling scheme the iterates over the following conditional distributions

All the conditional distributions have standard forms, and details can be found in GhoshClyde2011.

Approximate computation with hierarchical priors

Approximate inference methods typically involve optimization algorithms for approximating posterior moments (typically the mean or the mode, and the variance) instead of sampling from the full posterior distribution. Then one can proceed their analysis using only these moments, treating the Bayesian estimator similar to a frequentist point estimator. This approach to Bayesian inference has been popularized in computing science, e.g.\ in estimation of high-dimensional Bayes networks, where large datasets is the norm and MCMC inference is extremely costly. An obvious critique of approximate Bayesian inference of this sort, is that we can't fully take into account parameter uncertainty by characterizing the full parameter posterior distribution. However, it is in high-dimensional and the so-called ultra high-dimensional models Shinetal2018, that shrinkage and sparsity via a hierarchical prior is necessary. Since analytical results for most classes of hierarhical priors are not available, and Monte Carlo sampling is costly in very high dimensions, it is not surprising that approximate methods have become very popular. Additionally, at the conceptual level, we saw that only Bayesian posterior medians/modes correspond to penalized likelihood estimators, but there are not always good theoretical guarantees for the tails of the posterior.\footnote{See for example our discussion of the results of castillo2015 in the Bayesian lasso prior.} Finally, the concept of sparsity is indeed more interpretable in a setting with point estimates of coefficients, that is, it is more straightforward to test and interpret $H_0: \beta_j = 0$ when $\beta_j$ is approximated with a point estimate rather than when we have thousands of samples from the full posterior of $\beta_j$.

All these reasons lead us to review some of the most popular approximate methods for posterior inference. While many of these methods have been popularized and used extensively in computing science often without theoretical justifications, investigation of their theoretical/asymptotic properties is currently a topic of vivid research in mainstream statistics.

Variational Bayes

Variational Bayes is probably the most prominent of algorithms, at least when it comes to high-dimensional inference using hierarchical priors. The idea behind this class of algorithms is rather simple, and under certain assumptions (what we will call mean-field approximation later) variational Bayes algorithms can be fairly simple to implement by practitioners who are familiar with the Gibbs sampler. For notational simplicity, assume we have a vector of parameters $\bm \theta$ with support $\Theta$ and data $\bm D$, resulting to the posterior distribution

Assume that this posterior is intractable because the data density $p(\bm D \vert \bm \theta)$ is complex (e.g.\ it is a highly nonlinear function, or it has unidentified parameters), or because the prior $p(\bm \theta)$ is complex (non-conjugate), or because $\theta$ is high-dimensional (in which case the posterior is a high-dimensional function), or due to combinations of the above cases. In such settings, MCMC is not only computationally costly, but it can also become numerically unstable/unreliable.\footnote{For example, consider the case of a high-dimensional nonlinear regression with highly correlated predictors. In this example, unless modifications are introduced such as adaptive tuning, mixing and convergence of standard Gibbs sampler algorithms with hierarchical priors will tend to be slow.}

The idea behind variational Bayes is to introduce a family of simpler, approximate densities over the parameters $\bm \theta$ which is denoted by the set $\mathcal{Q}$. The objective is to find a member of the family $q(\bm \theta) \in \mathcal{Q}$ that is as close as possible to the true posterior. “Closeness” is measured by the Kullback-Leibler (KL) divergence, and the optimal density $q^{\star}(\bm \theta)$, among all densities $q(\bm \theta)$, is the one that minimizes this criterion:

In the above formula the KL measure on the RHS is defined as {

}where all expectations are w.r.t $q(\bm \theta)$, for example, $E_{q(\bm \theta)} \left( \log(p(\bm \theta \vert \bm D)) \right) = \int_{\bm \theta \in \Theta} q(\bm \theta) \log(p(\bm \theta \vert \bm D)) d \bm \theta$. In the last equation we have used the fact that $ E_{q(\bm \theta)} \left(\log(p(\bm D)) \right) = \log(p(\bm D))$ since $p(\bm D)$ does not involve $\bm \theta$. For the same reason the variational Bayes minimization problem is equal to minimizing the difference between the first and the third terms in equation (ref). We can also solve this equation for the log marginal likelihood $\log(p(\bm D))$ to show that

where we define evidence lower bound (ELBO) to be the quantity $ELBO = + E_{q(\bm \theta)} \left(\log \left( p(\bm D \vert \bm \theta)p(\bm \theta) \right) \right) - E_{q(\bm \theta)} \left(\log(q(\bm \theta)) \right) $. The ELBO has this name exactly because it is a lower bound for the log evidence (marginal data density). This is because in the equation above the KL divergence term is non-negative, such that $\log(p(\bm D)) \geq ELBO$. Therefore, the optimal $q^{\star}(\bm \theta)$ can be found by equivalently maximizing the ELBO criterion function.

{ The CAVI algorithm} \\ When latent variables are present, optimizations such as maximizing the ELBO criterion can be implemented using the popular expectation-maximization (EM) algorithm, where the complete log likelihood is computed (E-step) and then it is maximized (M-step). However, in Bayesian inference all parameters are latent (random) variables and as the optimization problem above involves optimizing over the functional $q(\bm \theta)$ and not $\bm \theta$ itself, the EM algorithm is not appropriate. Variational inference instead requires to choose the variational family of distributions $\mathcal{Q}$ and then maximize the ELBO. In most cases this can be done iteratively, with certain schemes that resemble the EM algorithm (but are not identical to EM), and convergence is guaranteed to a local maximum and if the likelihood is log-concave then to a global maximum. The simplest algorithm for maximizing the ELBO is called Coordinate Ascent Variational Inference (CAVI). Its simplicity comes at the cost of certain simplifying assumptions. The first one is that $\mathcal{Q}$ must strictly belong to the exponential family of distributions (e.g.\ the normal satisfies this condition, but the Student's t does not). The second restriction is the use of the mean-field approximation that postulates that the proposed posterior distribution $q(\bm \theta)$ can be decomposed into $M$ independent groups of the form

where the groups could either have $\bm \theta_{m}$ being a scalar or a vector. The estimated variational posteriors will be independent, meaning that the mean-field approximation/factorization implies that $\bm \theta_{m}$ will be a-posteriori uncorrelated with $\bm \theta_{k}$, for $k \neq m$ and $k,m = 1,...,M$. This assumption in several modeling settings can be harmless, but in several others it can become harmful -- we discuss this issue in detail later when we examine variational Bayes inference in a linear regression with variable selection prior.

Under the assumption of the mean field approximation it can be shown Bleietal17 that the optimal densities $q_{m}(\bm \theta_{m})$ satisfy

where $E_{q_{(-m)}(\bm \theta_{(-m)})}()$ means that the expectation is w.r.t all variational densities except $q_{m}(\bm \theta_{m})$. Broadly speaking this formula says that in order to optimize w.r.t. $q_{m}(\bm \theta_{m})$ we need to evaluate the posterior under the assumption that all other parameters $\bm \theta_{(-m)}$ are fixed to their posterior expectation (posterior mean). We would obviously need to iterate through (ref) for each $m=1,...,M$ keeping all other parameters fixed to their posterior means, but it can be shown that such iteration results in increasing the ELBO criterion. If the ELBO hasn't changed from one iteration to the next, the algorithm has converged. This criterion resembles the EM algorithm that converges when the value of the likelihood in subsequent iterations is approximately similar. The fact that $\bm \theta_{m}$ is updated conditional on fixing all other parameters $\bm \theta_{(-m)}$ makes variational Bayes resemble Gibbs sampling inference -- despite the fact that there is no sampling involved. We next derive a CAVI algorithm for a linear regression model with variable selection prior, in order to clearly demonstrate how the mean-field approximation is applied and how the functions in (ref) look like.

{ A variational Bayes approach to variable selection} \\ This subsection follows closely the analysis of Ormerodetal2017, and the reader should consult this paper for extensive discussion and proofs. Consider the KuoMallick1998 regression we explored in (ref) and is of the form

where $\bm \Gamma = diag(\gamma_{1},...,\gamma_{p})$ and $\bm D$ is a diagonal prior covariance matrix, e.g.\ $\bm D = c \times \bm I_{p}$ for some constant $c$. Therefore, according to the information above the joint prior is decomposed into $p \left( \bm \beta, \left\lbrace \gamma_{j} \right \rbrace_{j=1}^{p}, \sigma^{2} \right) = p(\sigma^2)\prod_{j=1}^{p}p(\beta_{j}) \times p(\gamma_{j}) $ meaning that all parameters are a-priori uncorrelated. Such choices are both conceptually and practically fine, first because we don't have prior information on how parameters are correlated, and second because we can construct very powerful shrinkage and variable selection algorithms based on these forms. It is important for the posterior to allow the parameters to be correlated, as this posterior correlation will come from information in the data likelihood. We discussed previously that the mean-field factorization implies that some groups of parameters will be uncorrelated a-posteriori. Ormerodetal2017 look into three different ways of applying the mean field factorization, based on how we want to define the groups $m=1,...,M$, and their implications for posterior inference. These factorizations are the following

The factorization (A) means that application of formula (ref) to the set of parameters $(\bm \beta, \bm \gamma)$ given $\sigma^{2}$ gives

or, similarly, that

where $\lambda = \log\left(\frac{\pi_0}{1-\pi_0}\right)$ and $\kappa=E_{q} \left(1/\sigma^2 \right)$. Therefore, we can easily obtain the conditional variational density of $\bm \beta$ and the (marginal) variational density of $\bm \gamma$ as

where $\bm V_{\bm \gamma}=\left( \kappa \bm \Gamma \bm X^{\prime} \bm X \bm \Gamma \bm \beta + \bm D^{-1} \right)^{-1}$ and $\bm \mu_{\bm \gamma} = \bm V_{\bm \gamma} \bm \Gamma \bm X^{\prime} \bm y$. In the second equation we only have the kernel of $q \left( \bm \gamma\right)$, but we can easily normalize this to integrate to one by dividing with the sum of the density of all possible combinations of $\bm \gamma$ (which is $2^{p}$ in the regression with $p$ covariates). In order to derive the marginal variational posterior density of $\bm \beta$ we need to integrate out the $\bm \gamma$, which is easily done as these are binary indicators. This marginal density is of the form

which is a combinatorial sum over all $2^{p}$ outcomes for the vector $\bm \gamma$. This sum can be evaluated in finite time only for small $p$. However, for small $p$ there are other numerous analytical algorithms that can be used, for example, we can use a $g$-prior and obtain marginal likelihoods analytically for all $2^{p}$ models, in which case there is no point in using variational Bayes. Therefore, this mean-field factorization is not useful.

The mean-field factorization/approximation (B), which was used in CarbonettoStephens12, provides a scalable variational Bayes algorithm where the $q(\beta_{j})$ can be estimated independently and efficiently (by means of parallelization) for each $j=1,...,p$. However, posterior variances of these regression coefficients will tend to be underestimated exactly because of this assumption of posterior independence. Unless the predictors in $\bm X$ are uncorrelated (which is not realistic for economic data, and for large-$p$ settings), the bias in posterior variances can be substantial.

The most reasonable case, which is the choice of Ormerodetal2017, is case (C). Under this factorization, one full implementation of the iterations in (ref) looks like this:

In the equations above we have used some matrices vectors/matrices that are based on $\pi_{j}$, namely $\bm \pi= (\pi_{1},....,\pi_{p})^{\prime}$, $\bm \Pi = diag(\bm \pi)$ and $\bm \Omega = \bm \pi \bm \pi^{\prime} + \bm \Pi(\bm I - \bm \Pi)$. The symbol $\odot $ denotes the Hadamard product. We used the notations that for a general matrix $\bm A$, $\bm A_j$ is the $j$th column of $\bm A$, $\bm A_{(-j)}$ is $\bm A$ with the $j$th column removed, $A_{i,j}$ is the $(i,j)$th entry of $\bm A$, $\bm A_{(-i),j}$ is the vector corresponding to the $j$th column of $\bm A$ with the $i$th component removed.

The formulas look quite similar to the conditional posteriors in the KuoMallick1998 Gibbs samplers, but there is no sampling involved. Instead, when the variational posterior mean and variance of $\bm \beta$ are calculated, $\sigma^{-2}$ is fixed to its posterior mean $\kappa$ and the same for $\gamma$ (the vector $\bm \pi$ and its variants, i.e.\ $\bm \Pi$ and $\bm \Omega$). Given that in the very first iteration $\kappa$ and $\bm \pi$ will be initialized to some random values, a convergence period is required until we end up with final estimates of the posterior moments of all parameters.

{ Further readings} \\ There are numerous papers on variational Bayes inference in computing science problems, for example, natural language processing (text analytics) and Bayesian networks. In statistics and machine learning there has been a consistent effort to establish consistency and other properties of variational Bayes estimates; see for example Giordanoetal2018 andWangBlei2019. With regards to high-dimensional regression and hierarchical priors, the contributions of CarbonettoStephens12, Ormerodetal2017 and Nevilleetal2014 are an excellent starting point. KoopKorobilis2018 provide a variational Bayes algorithm for a dynamic spike and slab prior for models featuring time-varying parameters and stochastic volatility (see also next section).

EM algorithm

We discussed previously how the expectation-maximization (EM) algorithm is not appropriate for variational Bayes inference, since there we are looking to find the “best” density function of $\bm \theta$, rather than a point estimate of our parameters $\bm \beta$. However, the EM algorithm can be used to find the posterior mode of $p(\bm \theta \vert \bm D)$, an inference method known as maximum a-posteriori (MAP) inference. The mode of the posterior under a diffusing (flat) prior distribution is identical to the maximum likelihood estimate, while the MAP estimate under a hierarchical prior corresponds to a penalized likelihood estimator. Therefore, MAP inference -- which was also popularized in computing science -- can be thought of as a bridge between Bayesian and maximum/penalized likelihood inferences where it combines the strengths of both approaches.

There are numerous implementations of MAP inference using the EM algorithm but, unlike the Gibbs sampler, in many cases algorithms are model-specific and cannot generalize easily. With regards to variable selection and shrinkage we indicatively mention the key contributions of CaronDoucet2008, Figueiredo2003 and GriffinBrown2011. A notable recent contribution is the EM variable selection (EMVS) of RockovaGeorge2014. These authors adopt a setting (likelihood and prior) that is identical to GeorgeMcCulloch1993 but they use the EM algorithm as a means of lowering the computational burden of Markov-Chain Monte Carlo methods when estimating posterior distributions over subsets of potential predictors.

Other approximate algorithms

There are several other algorithms for approximate high-dimensional inference. These include parallel MCMC, Hamiltonian Monte Carlo, Approximate Bayesian Computation (ABC), Expectation propagation, and Message Passing. A review of all these classes of algorithms can be found in KorobilisPettenuzzo2020. A few representative works relying on such algorithms are DehaeneBarthelme2018, KimWand2016, Korobilis2021, Liuetal2019, WainwrightJordan2008 and Zouetal2016.

Monte Carlo exercise: Conjugate vs independent hierarchical priors

Should we be using a conditional or unconditional hierarchical prior (see (ref))? In order to investigate this question, we consider simulation by generating data from the regression model $\bm y = \bm X \bm \beta + \bm \epsilon$ with $\bm \epsilon \sim N_n(\bm 0,\sigma^2 \bm I_n)$. We let $n=100$ and $\sigma^2=3$. We construct the true vector of slope parameters $\bm \beta =c \tilde{\bm \beta}$ by assigning values $\{1.5, -1.5, 2, -2, 2.5, -2.5 \}$ to the first 6 elements of $ \tilde{\bm \beta}$ and setting others to zero. We choose a constant $c>0$ to achieve a desired level of signal-to-noise ratio \footnote{ In a general linear regression $\bm y = \bm X \bm \beta+ \bm \epsilon$, the signal-to-noise ratio (SNR) is defined as $SNR=\frac{|| \bm \Sigma_X^{1/2} \bm \beta ||^2 }{ \sigma^2 }$ where $\sigma^2$ is the error variance and $\bm \Sigma_X$ is a $p \times p$ covariance matrix of $\bm X$. $||\bm \Sigma_X^{1/2} \bm \beta||^2 = \bm \beta' \bm \Sigma_X\bm \beta$ measures the overall signal strength. A related quantity is $R^2_{pop}$, the population value of $R^2$, defined as $\frac{SNR}{1+SNR}$. }. The data matrix $\bm X$ is generated from the multivariate normal distribution with mean zero and covariance matrix being an identity matrix. The covariates are standardized for estimation. We examine different values of the number of covariates $p\in \{50, 100,300 \}$ and the signal-to-noise ratio $R^2_{pop}\in \{ 0.4,0.8 \}$.

We consider three shrinkage priors (1) student-t, (2) lasso, and (3) horseshoe and compare performance under the conditional and independent priors. The analysis was repeated 100 times with new covariates and responses generated each time. For each, the metrics recorded were: the estimated value of $\sigma^2$, the bias and MSE of the first 6 elements of the coefficients vectors, the number of false negatives (FN), the number of false positives (FP), and the number of true positives (TP). Posterior means were used as point estimates of the slope coefficients and the error variance. We utilize the post-processing approach of LiPati2017 in order to categorize the covariates into signals and noises.

(ref) summarizes the results. Panel (a) shows results when the signal is relatively strong i.e.\ $R^2_{pop} =0.8$. Both conjugate and independent priors do well in terms of TPs and FNs. However, there are some notable differences. First, the error variance tends to be underestimated when conjugate priors are used. This was in fact pointed out by Moranetal2019. Intuitively, conjugate priors implicitly add $p$ “pseudo-observations” to the posterior (compare (ref) with (ref)) which can result in underestimations of $\sigma^2$ when $\bm \beta$ is sparse. Second, the independent priors tend to have larger bias and MSE of the signals under the high-dimensional case (i.e.\ $p=300$). ParkCasella2008 point out that independent priors can induce bi-modality of the posterior on the slope coefficients. This can make the posterior distributions for $\bm \beta$ more spread than in the conjugate case. We also see that independent priors have larger FPs when $p=300$, which could be a result of this. Panel (b) shows results under relatively weak signal i.e.\ $R^2_{pop} = 0.4$. We see that all methods face difficulty with distinguishing signals with noise (see FNs, FPs, and TPs) and have large bias and MSEs, compared to the case with $R^2_{pop} = 0.8$. However, the general findings on the difference between conjugate and independent priors are the same: conjugate priors tend to underestimate $\sigma^2$ while independent priors tend to have higher bias and MSE of the signals when $p$ is large compared to $n$. We encourage researchers to be aware of these issues when choosing priors and to conduct sensitivity checks.

\FloatBarrier

Bayesian shrinkage and variable selection beyond linear regression

So far we explored variable selection in the high-dimensional linear regression, also known as “large $p$, small $n$” regression. This setting is already flexible enough, as there are various cases of generalized linear models (GLMs) that have conditionally linear forms. The purpose of this section is to demonstrate that there is an even larger list of models where hierarchical priors have immediate applicability. In particular, we explore key applications of hierarchical shrinkage and variable selection priors in vector autoregressions, factor models, time-varying parameter models, high-dimensional confounder selection in models for treatment effects, and Bayesian quantile regression. This list is far from exhaustive\footnote{For example, one application of Bayesian shrinkage and selection that we do not cover in this section, but is of importance in statistics and in finance, is high-dimensional covariance matrix estimation and selection, see WangPillai2013 as an indicative example. Another topic we won't cover here, but is becoming increasingly very important in statistics and econometrics, is Bayesian additive regression trees (BART). For an up-to-date review of the topic see Hilletal2020.} and its only purpose is to illustrate how Bayesian computation simplifies high-dimensional inference in unconventional settings.

Vector autoregressions

The most popular working model for economists using time series variables is the vector autoregression (VAR). VARs are used for the joint modeling of the dynamics of many macroeconomic and financial time series, $\bm Y$, allowing analysts and policy-makers to answer questions regarding dynamic responses of variable $\bm Y_{i}$ to a shock in some other variable $\bm Y_{j}$, $i \neq j$. This is a very important tool especially when variable $\bm Y_{j}$ is controlled by the policy-maker. For example, the central bank controls the short-term interest rate as well as other quantities related to monetary policy, while government controls taxes and fiscal policy in general. Due to the fact that availability of time series observations for macroeconomic and (low-frequency) financial data is limited\footnote{Especially in countries other than the US, where statistical agencies might have available only a handful of decades of data; e.g.\ euro area time series typically begin in 1995 or 1999.}, estimation of macroeconomic VARs by and large relies on Bayesian shrinkage priors. Additionally, parameters in VAR models proliferate at a polynomial rate as the dimensions of the model increases. In univariate linear regression settings, a model with twice as many exogenous predictors has twice as many parameters to estimate. There is not such an analogy in VARs where all variables are endogenous and each variable (and its lagged terms) affects all other variables in the system.

Unlike our previous notation, consider time series observations $t=1,...,T$ and an $n$-dimensional vector of variables $\bm Y_{t}$, that is, $n$ denotes the number of variables of interest (and not the number of observations anymore) with $\bm Y = \left[ \bm Y_{1}^{\prime},...,\bm Y_{T}^{\prime} \right]^{\prime}$ is a $T \times n$ data matrix. The VAR model for $\bm Y_{t}$ with $p$ lags, also denoted as VAR($p$), is of the form

where $\bm c$ is an $n \times 1$ vector of intercepts, $\bm A_{i}$ are $n \times n$ matrices of lagged terms for each $i=1,...,p$, and $\bm E_{t} \sim N(\bm0,\bm \Sigma)$ with $\bm \Sigma$ an $n \times n$ symmetric, semi-positive definite covariance matrix. The VAR is a heavily parametrized model: it has $(1 + np)n$ coefficients $\bm B = \left[\bm c, \bm A_{1}, ..., \bm A_{p} \right]$, plus another $n(n+1)/2$ unique elements in the covariance matrix $\bm \Sigma$. For example, the largest VAR model specified in Koopetal2019 has $n=129$ and $p=13$ which implies that the total number of parameters is in excess of $200,000$.

Vector autoregressions are effectively linear regression models with parameter matrix $\bm B = \left[\bm c, \bm A_{1}, ..., \bm A_{p} \right]$ and data matrix $\bm X_{t} = \left[ \bm 1, \bm Y_{t-1},..., \bm Y_{t-p} \right]$, such that application of hierarchical priors for shrinkage and variable selection is fairly straightforward. The task of sampling from the conditional posterior of the regression coefficients $\bm B$ can be further simplified if the VAR is written in seemingly unrelated regressions (SUR)\footnote{For a thorough and accessible introduction to Bayesian inference in VARs see KoopKorobilis2010.} form

where $vec(\bullet)$ is the operator that stacks the columns of a matrix into a single column vector. That way, $ \bm y = vec \left( \bm Y \right)$ is a $Tn \times 1$ vector where the first $T$ elements are the observations of the first variable, the next $T$ rows correspond to observations of the second variable, and so on up to variable $n$. The measurement matrix $\bm Z = \left(\bm I \otimes \bm X\right)$ is a block-diagonal matrix with the $T \times (1 + np)$ matrix $\bm X$ repeating on its diagonal $n$ times. The formulation above is observationally identical to the one in (ref), since there are no new parameters or data introduced, but it has the benefit that VAR parameters show up as the $(1 + np)n \times 1$ vector $\bm b = vec\left( \bm B\right)$. Therefore, the SUR form in (ref) is identical to a univariate regression model, even though this model has both many covariates but also many observations ($\bm y$ and $\bm Z$ both have $Tn$ rows, instead of $T$ rows in a univariate regression). Therefore, it is straightforward to define any hierarchical prior we desire for the vector of VAR parameters $\bm b = vec\left( \bm B\right)$ and derive conditional posteriors, despite the fact that Gibbs sampling might become quite cumbersome as the dimension $n$ of the VAR increases.\footnote{See for example, Korobilis2013b,Korobilis2016 and KoopKorobilis2016.}

In large $n$ cases and when shrinkage on the covariance matrix $\bm \Sigma$ is needed, Carrieroetal and Koopetal2019 propose to estimate the VAR equation-by-equation. For example, in the formulation of Koopetal2019 one can write (ref) as

where $\bm P$ is a lower triangular matrix with ones on its main diagonal (unitriangular) that satisfies the LDL-decomposition $\bm \Sigma = \bm P \bm D \bm P^{\prime}$. In this case $V_{t} \sim N(\bm0,\bm D)$ where $\bm D$ is a diagonal matrix with variance elements $d_{ii}^{2}$ on its main diagonal, $i=1,...,n$. This formulation is equivalent to (ref) because $ \bm P \bm V_{t} \sim N(\bm 0, \bm P \bm D \bm P^{\prime}) = N(\bm 0,\bm \Sigma) \overset{d}{=} \bm E_{t}$. Since by construction $\bm P$ is invertible, with $\bm P^{-1}$ also a unitriangular matrix, we can write

where $\bm \Gamma = \bm P^{-1}\bm B$, and we have split $\bm P^{-1}$ into an identity matrix and a lower triangular matrix $\widetilde{\bm P}^{-1}$ by means of the equation $ \bm P^{-1} = \bm I + \widetilde{\bm P}^{-1}$. In (ref) we have a VAR on $\bm Y_{t}$ where the covariance matrix elements in the original covariance matrix show up as contemporaneous elements of $\bm Y_{t}$ itself on the right-hand side in the term $ - \widetilde{\bm P}^{-1} \bm Y_{t}$. Notice that in matrix form this is a nonlinear system as $\bm Y_{t}$ shows up both on the left hand side and the right hand side of (ref). However, exactly because $\widetilde{\bm P}^{-1}$ is lower triangular and $\bm V_{t}$ has a diagonal covariance matrix $\bm D$, equation-by-equation estimation is feasible. In particular, each VAR equation $i$ depends on lags of all other equations and contemporaneous terms in the previous $i-1$ equations. This means that estimation of the VAR collapses to estimation of $n$ independent univariate models, such that specification of hierarchical priors and MCMC estimation are also straightforward. The additional benefit from this approach is that hierarchical priors can be specified to the elements of $\widetilde{\bm P}^{-1}$, thus leading to shrinkage or sparse estimation of the original VAR covariance matrix $\bm \Sigma$, similar to the methodology of SmithKohn2002. More details of this approach can be found in Koopetal2019, while variants of this approach have been proposed in Baumeisteretal2020, Carrieroetal, and KorobilisPettenuzzo2019.

Factor model shrinkage and selection

Factor models have a long history in econometrics, and an even longer history in statistics and psychology/psychometrics. For that reason, while there are some popular formulations across different literatures, there are also different variations based on the data and applications. We first establish some key results for the specific case of the so-called static factor model, and subsequently we review some of the most popular factor models used in economics and finance. We end our discussion with strategies for Bayesian shrinkage and variable selection in this class of models.

Consider an $n \times 1$ vector of economic variables $\bm X_{t}$ observed over $t=1,...,T$ (without loss of generality $t$ can measure time series, but it can also be observations on individuals or other cross-sectional units). The dimension $n$ can be inconveniently high\footnote{This description includes both the ultra high-dimensional case where $n$ can be in the order of thousands, or larger, but also the case where $n$ is small but much larger than the number of available observations $T$.}, such that unrestricted estimation of models (e.g.\ linear regressions) using the data $\bm X_{t}$ is infeasible. Our target is to estimate a lower-dimensional $k \times 1$ vector ($k<<n$) of latent variables (factors), that summarizes as much as possible the information contained in $\bm X_{t}$. For that reason we define the following multivariate model

where $\bm \Lambda$ is an $n \times k$ matrix of parameters, $\bm F_{t}$ are the latent variables and $\bm \varepsilon$ is a disturbance term. This is not a regression model, as both $\bm \Lambda$ and $\bm F_{t}$ are latent. For simplicity, we follow LopesWest2004 and make the assumption that $\bm F_{t} \sim N_{k}(\bm 0, \bm I)$, although we can allow the factors to have a more general covariance matrix such that they are correlated with each other. For the disturbance term we assume $\bm \varepsilon_{t} \sim N_{n}(\bm 0, \bm \Sigma)$.

While the model in (ref) looks like a linear regression -- in which case application of hierarchical priors would be straightforward -- this is not the case due to the fact that all terms on the right-hand side of the equation are latent. In many instances, researchers in economics, finance and other disciplines replace $\bm F$ with the first $k$ principal components (PCs). PCs are nonparametric estimates of the factors, that is, they are only approximate estimators of the true factors implied by the likelihood of the model in (ref). Plugging in the place of the latent factors (parameters) $\bm F$ the PC estimates turns the factor model into a regression model and inference is simplified. Conditional on the principal component estimates, $\bm \Lambda$ and $\bm \Sigma$ can be estimated simply via least squares, but also standard Bayesian methods for multivariate regression can be used. The benefit of this two-step approach is that it is simple, both conceptually and computationally, and that principal component estimates always provide a sensible fit. However, in many more complex settings (e.g.\ macroeconomic dynamic factor models or financial factor models with stochastic volatility), the PC only provide a rough approximation, and it might be preferable to use the likelihood function to estimate $\bm F$. In such cases, it is imperative to make sure the factor estimates are unique. Therefore, we discuss first how to uniquely identify the factors, loadings and other parameters, before discussing Bayesian inference using hierarchical priors.

Identification of the factor model

The factor model implies that the conditional covariance matrix of $\bm X$ can be decomposed as

This decomposition illustrates the fact that likelihood-based (maximum likelihood or Bayesian) estimation of the factor model suffers from lack of identification of a unique set of parameter estimates. For example, consider the case where $\bm \Sigma$ is a full matrix, then there are infinite ways to construct the decomposition in (ref). In order to deal with this issue, it is common in factor analysis to set $\bm \Sigma$ to be a diagonal matrix.\footnote{This is known as the exact factor model assumption, as opposed to the class of approximate factor models that allow for “some” weak correlation among the variables in $\bm X$ and a $\bm \Sigma$ covariance that has some non-zero off-diagonal elements. Approximate factor models are typically not estimated with likelihood-based methods, so we don't consider this class of models here.} A consequence of this assumption is that the disturbances $\bm \varepsilon_{j}$ become idiosyncratic to each variable $\bm X_{j}$, $j=1,...,n$, that is, they capture measurement errors and other idiosyncrasies of each variable that are not attributed to its covariation with the remaining $n-1$ variables. Instead, any comovements/commonalities in the $n$ variables $\bm X$ are solely captured by the common component $\bm \Lambda \bm F$.

Having $\bm \Sigma$ diagonal is a big step towards identification in the factor model. As LopesWest2004 mention, $\bm \Omega$ has $n(n+1)/2$ unique elements, therefore, the number of elements in the decomposition of (ref) should not exceed that threshold. The matrix $\bm \Lambda$ has $nk$ elements, and $\bm \Sigma$ being diagonal has $n$ elements, therefore we obtain the inequality $n(n+1)/2 \geq nk + n$, which provides an upper bound on the number of factors one can extract: with $n=5$ variables we can extract $k=2$ factors, and when $n=20$ the maximum number of factors that can be extracted is $k=9$. However, there is a further problem impairing identification of the factor model, and this pertains to separating $\bm \Lambda$ from $\bm F$. Without further restrictions, there are infinite ways of finding such matrices that provide exactly the same values for the common component. Put more formally, if $\bm P$ is an $k \times k$ orthogonal matrix such that $\bm P \bm P^{\prime} = \bm I_{k}$, then the factor model can be rewritten as

where $\bm \Lambda^{\star}$ and $\bm F^{\star}_{t}$ are alternative estimates to $\bm \Lambda$ and $\bm F_{t}$ that provide exactly the same likelihood value (they are observationally equivalent). Given that the variance of the factors is normalized to be one, unique identification of the loadings and the factors requires an additional $k(k-1)/2$ restrictions on the loadings matrix $\bm \Lambda$. A standard restriction that is imposed in this case LopesWest2004 is to restrict $\bm \Lambda$ to be lower triangular, that is, the top $k \times k$ block of this matrix has its $k(k-1)/2$ upper triangular elements equal to zero. This restriction provides local identification up to a rotation of the sign, meaning that we could multiply any column of $\bm \Lambda$ with $-1$ and do the same to the respective column of the factors $\bm F$, and arrive to an observationally equivalent solution. For that reason, GewekeZhou1996 suggest to further assume the $k$ diagonal elements of $\bm \Lambda$ to be restricted to be positive.

When modeling comovements between financial time series, the assumption $\bm F_{t} \sim N_{k}(\bm 0, \bm I)$ is often not empirically relevant, and instead it is assumed that $\bm F_{t} \sim N_{k}(\bm 0, \bm \Sigma^{F})$ with $\bm \Sigma^{F}$ a diagonal matrix, with possibly heteroskedastic elements that capture changing (over time) volatility of financial variables Chibetal2006. In this case further restrictions on $\bm \Lambda$ are needed, and Chibetal2006 choose to fix the diagonal elements of the loadings matrix to be one. BBE2005 extract factors from a large macroeconomic dataset and in their methodology it is imperative for $\bm \Sigma^{F}$ to be a full covariance matrix (it is the covariance matrix of a VAR from which they want to identify shocks and estimate impulse response functions). Therefore, with $\bm \Sigma^{F}$ a full matrix, these authors further restrict the upper $k \times k$ block of the loadings matrix to be the identity matrix.

All these identification restrictions in various applications of the factor model do not come at no cost. Imposing zeros in the loadings matrix means that certain variables are excluded from determining the factors. In the case of BBE2005, in particular, the identity restriction means that the first variable exclusively loads on the first factor, the second variable exclusively on the second factor, and so on. Therefore, the ordering of the variables in $\bm X$ ends up affecting the estimates of $\bm F$, and in their case this restriction turns out to be empirically detrimental.\footnote{While their factors are statistically identified, they do not carry any economic content (i.e.\ they are not “structurally identified”). BBE2005 is one of the few papers that estimates a factor model both with principal components and least squares (the “plug-in” approach explained previously) and with Bayesian inference. Comparing the impulse response functions for some key variables using the two estimation methods, there are marked differences. Whenever PCA has been used to estimate the unknown factors, impulse response functions have the signs and shapes expected by economic theory. When likelihood-based factors have been used, the impulse responses of variables such as inflation degenerate to zero for all horizons; BBE2005.}

Bayesian shrinkage and variable selection in the factor model

Despite the fact that statistical identification of the factor model using zero and sign restrictions on the loadings might contradict evidence in the data, once the factor model is fully identified Bayesian inference becomes straightforward. To see this, we re-write for convenience the factor model including now all relevant identification restrictions that are imposed prior to estimation

LopesWest2004 show that this model is conveniently estimated sequentially via the Gibbs sampler, by sampling from conditional posteriors. The priors are of the form

where $\delta_0(\Lambda_{ij})$ is the Dirac delta function, that is, a point mass function for $\Lambda_{ij}$ that is concentrated at zero. The conditional posteriors are

where $\bm V_{F} = (\bm I + \bm \Lambda^{\prime}\bm \Sigma^{-1} \bm \Lambda)^{-1}$, $\bm V_{L,i} = (\bm D^{-1} + \frac{1}{\Sigma_{ii}}\bm F^{\prime}\bm F)^{-1}$ and $SSE_{i} = (\bm X_{i} - \bm F\bm \Lambda_{i})^{\prime}(\bm X_{i} - \bm F\bm \Lambda_{i})$. The $\vert \bullet$ notation above is used to denote conditioning on other parameters and data.

By updating $\bm \Lambda$ conditional on $\bm F$ and vice-versa, the Gibbs sampler works around the issue that the product of these two parameters shows up in the likelihood function. Of course, this sequential updating of the common component by updating each of $\bm \Lambda$ and $\bm F$ conditional on the other, will inevitably generate unwanted correlation in the Gibbs chain. In order to deal with the sampling inefficiency associated with correlated MCMC draws, Chibetal2006 proposed an alternative Metropolis-Hasting step for updating $\bm \Lambda$ conditional on the factors, while GhoshDunson2009 proposed a parameter-expanded Gibbs sampler RockovaGeorge2016. Notice that the fact that $\bm \Lambda$ has the required zero and sign restrictions imposed prior to estimation, means that every time we sample $\bm F$ conditional on $\bm \Lambda$ the factors will be sampled from a unique, identified posterior distribution. Had we not imposed these restrictions, the Gibbs sampler would still work numerically, but lack of identification means that each sample could correspond to different pairs of solutions for $\bm \Lambda$ and $\bm F$. That is, in the case where $\bm \Lambda$ and $\bm F$ are not separately identified, their product (the common component $\bm \Lambda \bm F$) is always identified. There are certain inference exercises, such as prediction, where it might be the case that identification and interpretation of the factors is not required; see for example the arguments in favor of this approach in BhattacharyaDunson2011 and in Korobilis2020.

An early attempt to full-Bayes inference in factor models is West2003 who developed a variable selection prior in the loadings of the static factor model. Under the LopesWest2004 identification scheme, the extension proposed by West2003 simply involves replacing the prior in (ref) with a spike and slab prior. As long as the identification restrictions are maintained, the presence of the variable selection prior can be used to find further data-based restrictions in the loadings matrix. Carvalhoetal2008 extend the variable ideas in West2003 to create a very sparse static factor model for genome data; see our discussion of their prior in (ref). KnowlesGhahramani2011 further extend these ideas to a spike and slab prior that is semiparametric, utilizing the ability of an Indian Buffet Process prior to allow for infinitely many factors. These authors use a Metropolis-within-Gibbs algorithm for inference. RockovaGeorge2016 propose a similar spike and slab formulation based on the Indian Buffet Process, but unlike KnowlesGhahramani2011 they propose maximum a posteriori (MAP) inference by means of approximating the posterior mode using an EM algorithm. The Bayesian asymptotic theory and posterior contraction rates for the sparse static factor model with continuous spike and slab priors is explored in detail in Patietal2014.

GhoshDunson2009 proposed a heavy-tailed prior on $\bm \Lambda$ (using a normal/inverse gamma mixture prior) which they argue performs better than the default normal/truncated-normal prior in (ref). BhattacharyaDunson2011 proposed a novel multiplicative gamma process prior on the factor loadings that shrinks more aggressively columns of $\bm \Lambda$ that correspond to a higher number of factors. They call their approach a sparse infinite factor model, as it allows to specify a maximum number of factors and the prior is able to determine zero and non-zero loadings, as well as the number of factors. The gamma process prior for the loadings matrix is of the following “global-local shrinkage” form

While the local shrinkage parameter is the same for each element of $\bm \Lambda$, the global shrinkage parameter $\tau_{j} $ is shrinking more aggressively as the index $j$ increases, where $j=1,...,k$ indexes the number of factors. This is because $\tau_{j} $ is a $j$-dimensional product of gamma-distributed random variables. Legramantietal2020 propose a cumulative shrinkage process prior and Srivastavaetal2017 propose a multi-scale generalized double Pareto prior; both these priors are similar in spirit to the BhattacharyaDunson2011 prior in terms of shrinking the loadings towards zero and selecting the appropriate number of factors at the same time.

We close this section by mentioning ongoing research on alternative solutions to the identification problem (rotational interdeterminacy) in factor models, that do not rely on preimposing zero restrictions on the loadings matrix. The expanded parametric forms proposed in papers such as BhattacharyaDunson2011 and Legramantietal2020 discussed above, deal with this issue efficiently. Other approaches include the ex-post processing approaches of Assmann2016 and KaufmannSchumacher2019. FruewirthLopes2018 introduce a generalized lower triangular representation of the factor model and propose a sparsity-inducing prior that overshrinks. While papers like West2003 apply sparsity after imposing zero identifying restrictions, the parameterization of FruewirthLopes2018 allows the prior to impose zeros that are sufficient for identification and inference, thus, not suffering from the rotation problem. Finally, Chanetal2018 propose an invariant parameterization of the static factor model that is based on the singular value decomposition.

Dynamic sparsity and shrinkage

When working in a time series setting the concepts of shrinkage, variable selection, and model averaging need not be static. This is true for economics where there has always been evidence that predictors can be unstable. There is significant theoretical and empirical evidence that when forecasting oil prices, stock prices, consumer prices, exchange rates, and numerous other economic/financial variables there is hardly a single exogenous predictor that can be claimed to be important over a substantial time sample. In practice, we observe “pockets of predictability”, that is, short periods where a specific variable might have predictive information for another variable of interest. This concept of unstable predictors has been popularized since the global financial crisis of 2007-2009, a period when it was obvious that all constant parameter relationships between economic variables completely broke down. Combined with the availability of new Bayesian tools for high-dimensional inference, a large literature has emerged since then that uses terms such as “time-varying sparsity” or “dynamic model averaging” or “time-varying dimension models”. KoopKorobilis2018 provide a detailed discussion of this literature.\footnote{At the same time, in the field of signal processing there is a related literature on “dynamic compressive sensing” for streaming signals (e.g.\ video); see ZinielSchniter13. }

The starting point for imposing dynamic sparsity and dynamic shrinkage is a regression with time-varying parameters (TVPs) and stochastic volatility (henceforth, abbreviated as TVP regression) of the form

where $y_{t}$ is a scalar time series observation for $t=1,...,T$, $\bm X_{t}$ is a $p$-dimensional vector of covariates (that can include an intercept, own lags of $y$ and exogenous predictors), $\bm \beta_{t}$ is a vector of time-varying (or drifting) regression coefficients, and $h_{t}$ is the scalar time-varying (or stochastic) variance/volatility parameter. Additionally, we assume $z_{t} \sim N(0,1)$, $\bm u_{t} \sim N_p( \bm 0, \bm Q)$ with $\bm Q$ a $p \times p$ covariance matrix, and $v_{t} \sim N(0,\delta^{2})$ with $\delta^{2}$ a scalar variance parameter.

As is the case with the constant parameter regression, shrinkage is mainly desirable in the TVPs $\bm \beta_{t}$, but this can take now a dual form: shrinkage towards time-invariance ($\bm \beta_{t}$ becomes a constant parameter) and “traditional” shrinkage towards zero.\footnote{Shrinkage of $h_{t}$ towards a constant variance specification is feasible, but it is not desirable for economic and financial time series data, since we know that economic shocks are very volatile and constant parameter specifications are always inferior (both using in-sample and out-of-sample measures of fit).} Notice that the TVP regression as is specified in equations (ref) - (ref) is a hiearchical model, and (ref) in particular can be though of as a hiearchical prior for $\bm \beta_{t}$ of the form $\bm \beta_{t}\vert \bm \beta_{t-1}, \bm Q \sim N(\bm \beta_{t-1}, \bm Q)$. Seen like this, it is straightforward to assume that $\bm Q$ is diagonal and allow its elements to follow of the hyperpriors we examined previously (e.g.\ Student-t, Laplace etc). However, doing so would only regularize the evolution of $\bm \beta_{t}$ around $\bm \beta_{t-1}$, where in the limit of $\bm Q = \bm 0$ then $\bm \beta_{t}$ becomes a constant parameter. Shrinking $\bm \beta_{t}$ towards zero requires different treatment, and there are numerous ways one can deal with this problem.

Dynamic variable selection or dynamic model averaging can be implemented in this setting by simply placing appropriate hierarchical priors that will allow to test the hypothesis $H_{0}:\beta_{jt}= 0$ vs $H_{1}: \beta_{jt} \neq 0$ for all $j=1,...,p$ and for all $t=1,...,T$. Recall that in “static” Bayesian model averaging the challenge is to average over $2^{p}$ regressions. Therefore, the dynamic version of model averaging implies that one has to average over $2^{p}$ regression models for each $t=1,...,T$. It is not surprising then that many proposed approaches in the literature for dealing with this problem are not based on computationally intensive MCMC algorithms. For example, KoopKorobilis2012 and DanglHalling2012 use variance discounting methods West1997 in order to provide plug-in estimators of $h_{t}$ and $Q$ and estimate a single time-varying parameter regression very quickly. Subsequently, dynamic variable selection and model averaging can be implemented by enumerating and estimating all $2^{p}$ possible models -- as long as $p$ is fairly small (e.g.\ 20 predictors).

In terms of directly introducing shrinkage and sparsity via hierarchical priors, there are numerous ways of doing so in a TVP regression model. Belmonteetal2014 and BittoFruewirthSchnatter2019 place hierarchical priors in an equivalent “non-centered” parameterization of the TVP regression that takes the form

where $\bm d_{t} \sim N_p(\bm 0,\bm I_p)$ and $\bm W$ is a diagonal matrix with elements $w_{j}$, $j=1,...,p$. This formulation is observationally equivalent to the TVP regression of equations (ref) - (ref). As long as the initial condition is $\bm \theta_{0}$ it holds that $\bm \theta + \bm \theta_{t} = \bm \beta_{t}$. This allows to split the time-varying parameter into a constant parameter level $\bm \theta$ (determined by data $\bm X_{t}$) and the additive time-variation around the constant level. Additionally, notice that the state equation is now standardized ($\bm d_{t}$ has unit variance) which can be done by setting $\bm W^{\prime -1}\bm W^{-1} = \bm Q$. Belmonteetal2014 place a Bayesian lasso (Laplace) prior on $\bm \theta$ and on the diagonal elements of $\bm W$. By doing so, they can shrink the total coefficient $\beta_{j,t}$ into a constant parameter $\theta_{j}$ (if $w_{j}\rightarrow 0$), or shrink it to zero (when both $\beta_{j,t}$ and $w_{j}$ are shrunk towards zero). Alternatively, the model can become an unrestricted TVP regression when both $\beta_{j,t}$ and $w_{j}$ are not shrunk towards zero by the Laplace prior.

NakajimaWest2013 convert the TVP regression into a latent threshold dynamic regression of the form

where $\bm S_{t}$ is a $p \times p$ diagonal matrix with element $s_{j,t} = I(\beta_{j,t}\geq d_{j})$. That is, the $s_{j,t}$ are 0/1 indicators that can shrink $b_{j,t}$ either towards zero or towards the unrestricted TVP $\beta_{j,t}$. The threshold value $d_{j}$ can be estimated endogenously such that the data decide which coefficients are zero (or not) at each point in time. Of course, similar to interpretation of spike and slab priors, the condition $I(\beta_{j,t}\geq d_{j})$ is a soft, rather than a hard, thresholding rule, due to the fact that $s_{j,t}$ (in a Bayesian setting) is a random variable. This means that once considering the full uncertainty in the posterior the approach of NakajimaWest2013 provides a class of smooth thresholding models; see also NakajimaWest2013JFE,NakajimaWest2015DSP,NakajimaWest2017BJPS.

RockovaMcAlinn2017 specify a dynamic spike and slab prior of the form

where $\bm \Gamma = diag(\bm \gamma)=diag(\gamma_{1},...,\gamma_{p})$ and $\lambda_0$ and $\lambda_1$ can also have further exponentional prior distributions, converting this prior into a dynamic version of the spike and slab lasso of RockovaGeorge2018. This prior is a spike and slab prior for $\bm \beta_{t}$ but it is only the slab component that incorporates the random walk evolution via the prior mean for $\bm \mu_{t}$. In contrast, KoopKorobilis2018 propose a similar but non-hierarchical prior of the form

where $\bm D_{\tau} = diag(\bm \tau^{2}) = diag(\tau^{2}_1,...,\tau^{2}_p)$ and $c$ is a small constant (set to $c=0.0001$ in KoopKorobilis2018). Again it is trivial to allow $\bm \tau^{2}$ to have its own hyperprior, such that we can combine shrinkage with sparsity in one setting. Finally, KalliGriffin2014 modify (ref) and introduce a normal-gamma mixture evolution for $\bm \beta_{t}$, which can be written in the following hyperprior form

which makes this a normal-gamma-Poisson mixture distribution. While this mixture having a very flexible distributional form, implying interesting shapes for $\bm \beta_{t}$, there might be sensitivity to the choice of the key hyperparameters $(\lambda_{j},\rho_{j}, \mu_{j})$.

All the examples above use the state-space form of the TVP regression and rely on recursive estimation methods, either in the form of the simple Kalman filter or (within the context of simulation methods) forward filtering backward sampling (FFBS) algorithms. However, as noted by Korobilis2021 one can simply discard the prior $\bm \beta_{t}\vert \bm \beta_{t-1}, \bm Q \sim N(\bm \beta_{t-1}, \bm Q)$ and treat the TVP regression as a constant parameter regression. This can be seen if we stack all observations in (ref) and write it as

where $\varepsilon_{t} \sim N(0,h_{t}^2)$. In this form, the TVP regression is a model with $T$ observations and $Tp$ covariates and it can be estimated as a high-dimensional “static” regression with data matrices $\bm y$ and $\mathcal{X}$ as defined in (ref). Korobilis2021 shows that a large class of hierarchical shrinkage priors can be placed on the $Tp \times 1$ parameter vector $\bm B$, and inference can proceed using the regression form in (ref) without the need to rely on state-space methods. Since $\bm B$ has $T$ time copies of parameters on the $p$ predictors in $\bm X_{t}$, more structured shrinkage can be placed by using a group lasso or other similar grouping prior.

Applying a shrinkage or variable selection prior directly to the vector of parameters $\bm \beta_{t}$ means that a certain $\beta_{j,t}$ might be unrestricted in period $s$, then restricted to zero in period $s+1$, then switch back to being unrestricted in $s+2$, and so on, for $s \in \{1,...,T\}$. This is a noisy approach to dynamic shrinkage/sparsity, and more persistent estimates over time might be desirable such that we prevent an important coefficient from becoming sparse just for a period or two, and vice-versa for a sparse coefficient. In many economic data, there is evidence of prolonged regimes where coefficients are either important or not important (e.g.\ macroeconomic recessions vs expansions, or bull vs bear stock markets). In this case, it might be desirable to incorporate the information in $\bm \beta_{t}\vert \bm \beta_{t-1}, \bm Q \sim N(\bm \beta_{t-1}, \bm Q)$ alongside a hierarchical shrinkage prior. A simple way to do this is to write the TVP regression as a static regression for the parameters $\Delta \bm \beta_{t} = \bm \beta_{t} - \bm \beta_{t-1}$. This takes the form{

} where we implicitly assume that $\bm \beta_{0}=0$ such that $\Delta \bm \beta_{1} = \bm \beta_{1}$. The $t$-th equation of the system above can be written as:

that is, equations (ref) and (ref) are observationally equivalent. Under this specification the prior implied by (ref) becomes (in matrix form)

and this prior can now be converted into a hierarchical prior by placing appropriate hyper-prior distributions on $\bm Q$.

Dynamic shrinkage and sparsity is a very active area of research, and there are several other important contributions that we don't explore here due to space constraints. For further readings we direct the reader to Chanetal2012, Irie2019, Kowaletal2019 and UribeLopes2017, among others.

High-dimensional causal inference

Let $y_i$ denote an outcome variable and $T_i$ be some treatment variable. Suppose that the $p$-dimensional vector of cofounders $\bm x_i$ is high-dimensional. The parameter of interest is the treatment effect $\alpha$ in the model below:

A naive post selection approach would be to apply the lasso to the equation above, excluding $\alpha$ from the $\mathcal{l}_{1}$-penalty and then regress $y_i$ on $T_i$ as well as on the selected covariates to estimate and conduct inference about the treatment effect. Any control variable that is highly correlated with $T_i$ but weakly with $y_i$ tends to drop out of the selection in the first stage, and could lead to omitted variable bias in estimating $\alpha$ in the second stage. Bellonietal2014 propose post-double selection to overcome such bias. Hahnetal2018 and Antonellietal2019 offer Bayesian counterparts in linear models, using shrinkage priors.

Using the model (ref) as a benchmark, Hahnetal2018 consider the following system of equations:

The likelihood can be factorized:

With re-parameterization $\left(\alpha, \bm \beta+ \alpha\bm \gamma, \bm \gamma \right)' \to (\alpha, \bm \beta_d,\bm \beta_c)'$, the system can be written as

The authors place independent shrinkage priors over $\bm \beta_c$ and $\bm \beta_d$:

where $C^+(0,1)$ denotes a folded standard Cauchy distribution. This prior is a proxy of the horseshoe prior. Non-informative priors are used for other parameters $\alpha \propto 1$, $\sigma_\epsilon \propto \frac{1}{\sigma_\epsilon}$, and $\sigma_v \propto \frac{1}{\sigma_v}$. They use a slice sampler for posterior sampling. Hahnetal2020 extends the approach to nonparametric case using regression trees.

Antonellietal2019 propose a spike-and-slab lasso prior approach. Their proposed framework is

where $\psi_0(\beta_j; \lambda_0, \sigma^2) =\frac{\lambda_0}{2\sigma}e^{-\lambda_0|\beta_j| / \sigma}$ and $\psi_1(\beta_j; \lambda_1, \sigma^2) =\frac{\lambda_1}{2\sigma}e^{-\lambda_1|\beta_j| / \sigma}$. $\lambda_1$ is fixed to be a small value, say 0.1, so that the prior variance in the slab component is high enough to be uninformative.

The hyperparameter $\lambda_0$ is chosen via empirical Bayes. A new feature that they introduce is the weights $\omega_j$ which are tuning parameters that they use to prioritize variables to be included (i.e.\ $\gamma_j=1$) if they are associated with the treatment. Specifically, they first fit the standard lasso on the model for predicting $T$ given $X$. If the $j$th covariate $x_j$ has non-zero coefficient from the lasso, they set $\omega_j=\delta$ for some $\delta \in (0,1)$. For other variables, $\omega_j=1$. On the one hand, a smaller value of $\delta$ leads to higher inclusion probability and hence more protection against the omitted variable bias. On the other hand, one needs to ensure a small enough inclusion probability for an unimportant variable in the outcome model. To balance the trade off, the authors choose $\delta \in(0,1)$ as the smallest value of $\omega_j$ such that the inclusion probability of $\beta_j=0$ is less than 0.1. See also Antonellietal2020, who introduce how posterior distributions of treatment and outcome models can be used together with doubly robust estimators.

Bayesian quantile regression

A regression specification can be represented more generally using the formulation

where $f(y_i \vert \bm X_{i} )$ is a conditional mean function for $y$ (conditional on covariates $\bm X$). Using this notation, the linear regression model can be recovered if we set $f(y_i \vert \bm X_{i} ) = E(y_i \vert \bm X_{i} ) = \bm X_{i} \bm \beta$, that is, the linear regression only models the (conditional) mean of $y$. The distribution of $y$ is fully determined by the assumptions we make about the disturbance term $\varepsilon$. In many cases it is desirable to use the information in covariates in such a way that the full distribution of $y$ is determined by $\bm X$. While such feature can also be incorporated implicitly in a traditional linear regression setting\footnote{For example, we can assume $\varepsilon_{i} \sim N(0,\sigma_{i}^{2})$, where $\sigma_{i}^{2}$ can be some function of $\bm X_{i}$.}, a structured (and popular) way is to model the conditional quantiles of $y$, $\mathcal{Q}_{r}(y_i \vert \bm X_{i} )$, where $r \in (0,1)$ denotes the quantile of $y$. While the conditional quantile can be modeled using either linear or nonlinear functional forms, the linear form is by far the most widely used.

In this case, we replace in (ref) $f(y_i \vert \bm X_{i} ) = \mathcal{Q}_{r}(y_i \vert \bm X_{i} ) = \bm X_{i} \bm \beta_{r}$ and obtain the following quantile regression specification

The model above is a linear regression for each quantile level $r$. KoenkerBassett1978 show that an estimator of this quantile regression model can be obtained as

where $\rho_{r}(u) = (r - \mathbb{I}(u<r))$ is a loss function. YuMoyeed2001 show that the same estimator $\widehat{\bm \beta}_{r}$ can be obtaining by obtaining the maximum likelihood estimator under the assumption that $\varepsilon_{i}$ is distributed as asymmetric Laplace, i.e.\ if it has the density

where $\sigma_{r}^{2}$ is a scale parameter. Therefore, the contribution of YuMoyeed2001 provides a parametric framework for implementing Bayesian inference.\footnote{Of course here we have similar conceptual issues as with the Bayesian representation of the lasso estimator: while Tibshirani1996 showed that the $\mathcal{l}_{1}$ optimization problem for the lasso is equivalent to the posterior mode of Bayesian regression estimator under a Laplace prior, castillo2015 show that the posterior distribution does not contract at the same rate as the posterior mode. Similarly here, there is an equivalence between quantile regression and maximizing the likelihood under an asymmetric Laplace likelihood as both problems provide unique point solutions. Bayesian inference, in contrast, assumes that coefficients are random variables and (unless one focuses on MAP or MMSE estimation) cannot be obtained as the solution to an optimization problem. In practice, however, it turns out that Bayesian quantile regression estimation using the asymmetric Laplace likelihood is a very flexible model, even if it is not identical to the model introduced by KoenkerBassett1978.} In particular, KozumiKobayashi2011 take advantage of the fact that the asymmetric Laplace likelihood can be written as a normal-exponential mixture of the form

where $z_{i,r} \sim Exp(\sigma^2_{r})$ and $u_{t} \sim N(0,1)$, while $\theta_r,\kappa_{r}^{2}$ are parameters defined as $\theta_{r} = \frac{1-2r}{r(1-r)}$, $\kappa_{r}^{2} = \frac{2}{r(1-r)}$. If we marginalise Equation (ref) over the exponentially distributed $z_{i,r}$ we obtain Equation (ref); see a derivation in KhareHobert2012.

Therefore, the Bayesian quantile regression model has the following representation

where $z_{i,r} \sim Exp(\sigma^2_{r})$ and $u_{t} \sim N(0,1)$ can either be thought of two disturbance terms (similar to frontier models in econometrics), or equivalently $u_{t}$ can be interpreted as the disturbance term and $z_{i,r}$ can be thought of as an unobserved factor (similar to the factors we examined in (ref)). In any case, conditional posteriors are trivial to derive as conditional on all other parameters, the posterior of each parameter of interest has standard form. This is easier to see for instance for the coefficients $\bm \beta_{r}$ by rewriting the model as

where $y_{i}^{\star} = (y_{i} - \theta_{r} z_{i,r})$ and $\sigma_{r}^{\star} = \sqrt{\sigma_{r}^{2} \kappa_{r}^{2} z_{i,r}} $. As long as we condition on $z_{i,r},\sigma_{r}^{2},\kappa_{r}^{2},\theta_{r}$, we can obtain a sample of $\bm \beta_{r}$ (assuming a normal prior), from the standard formulas for a linear regression of $y_{i}^{\star}$ on $\bm X_{i}$ with known variance $\sigma_{r}^{\star}$.

In more detail, we assume the following priors

where for simplicity assume that $\bm D_{\tau,r} = \bm D_{\tau} = \tau \times I_p$ where $\tau$ is fixed and known for all $r$. The conditional posteriors KhareHobert2012 are of the form

where the notation $\vert \bullet $ means “conditioning on other parameters and data”, $\bm V = \left(\bm X^{\prime} \bm U^{-1} \bm X + D_{\tau,r}^{-1} \right)^{-1}$, $\bm U = \left( \sigma_{r}^2 \kappa_r^{2} \right) \times diag \left(z_{1,r},...,z_{n,r} \right)$, $\widetilde{\bm y} = \left(\bm y - \theta_{r} \bm z_{r} \right)$, $y_{i}^{\star} = \left( y_{i} - \bm X_{i} \bm \beta_{r} - \theta_{r} z_{i,r} \right)$, and $a_r=n_{0,r} + \frac{3n}{2}$.

We need to note here a few important points regarding the implementation of this model

After taking the modeling points above into consideration, estimation of the Bayesian quantile regression model is not much more challenging than the normal (Gaussian) linear regression model. All the hierarchical priors described previously can be readily applied to the quantile model, with very minor modifications to the conditional posteriors presented in equations (ref) - (ref). The interesting feature is that because a separate regression needs to be run for each $r$, we can also use shrinkage or sparsity priors to allow different covariates in $\bm X$ to affect different quantiles of the distribution of $y$. Exactly because application of hierarchical priors is so trivial, we won't provide detailed examples and derivations here. KozumiKobayashi2011, while proposing a Gibbs sampler for the Bayesian quantile model, they also show in a subsection how easy it is to adopt the double-exponential (lasso) prior. Other applications include AlhamzawiYu2012, Yuetal2013, Korobilis2017 and RePEc:ecb:ecbwps:20212600, and the reader can consult these papers for modeling details. Limetal2020 is the case of a paper that estimates a Bayesian quantile regression using variational Bayes methods.

Concluding remarks

We have reviewed a wide range of concepts and algorithms for Bayesian sparse and shrinkage estimation. Our focus was on recent contributions in the field, covering the mainly academic publications during the decade 2010-2020, that are increasingly focusing on efficient computation and inference in high-dimensional models. A major contribution of our work is to collect in a single document all these recent contributions, as other reviews and surveys we are aware of provide only very focused reviews of certain hierarchical priors and algorithms. While we believe contributions to the field of Bayesian sparse and shrinkage estimation will keep expanding at a polynomial rate, we do hope that this review will become a useful manual for PhD students and researchers who want an accessible introduction to the field.

\addcontentsline{toc}{section}{References}