\section{Introduction} \label{Intro}

High-dimensional data are now routinely analyzed. In these settings, one often wants to impose a low-dimensional structure such as sparsity. For example, in astronomy and other image processing contexts, there may be thousands of noisy observations of image pixels, but only a small number of these pixels are typically needed to recover the objects of interest \cite{JohnstoneSilverman2004,JohnstoneSilverman2005}. In genetic studies, scientists routinely observe tens of thousands of gene expression data points, but only a few genes may be significantly associated with a phenotype of interest. For example, \cite{WellcomeTrust2007} has confirmed that only seven genes have a non-negligible association with Type I diabetes. Among practitioners, the main objectives in these scenarios are typically: a) identification (or \textit{variable selection}) of the non-negligible variables, and b) \textit{estimation} of their effects.

A well-studied model for sparse recovery in the high-dimensional statistics literature is the normal linear regression model,
\begin{equation} \label{linearregression}
	\boldsymbol{y} = \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon}, \hspace{.5cm} \boldsymbol{\varepsilon} \sim \mathcal{N}_n (\boldsymbol{0}, \sigma^2 \boldsymbol{I}_n ),
\end{equation}
where $\boldsymbol{y} \in \mathbb{R}^{n}$ is a vector of $n$ responses, $\boldsymbol{X} = [ \boldsymbol{x}_1, \ldots, \boldsymbol{x}_p] \in \mathbb{R}^{n \times p}$ is a design matrix of $p$ potential covariates, $\boldsymbol{\beta} = (\beta_1, \ldots, \beta_p)^T$ is a $p$-dimensional vector of unknown regression coefficients, and $\boldsymbol{\varepsilon}$ is the noise vector. When $p > n$, we often assume that most of the elements in $\boldsymbol{\beta}$ are zero or negligible. Under this setup, there have been a large number of methods proposed for selecting and estimating the active coefficients in $\boldsymbol{\beta}$. In the frequentist framework, penalized likelihood approaches such as the least absolute shrinkage and selection operator (LASSO) \cite{Tibshirani1996} are typically used to achieve sparse recovery for $\boldsymbol{\beta}$. In the Bayesian framework, spike-and-slab priors are a popular approach for sparse modeling of $\boldsymbol{\beta}$.

The spike-and-slab LASSO (SSL), introduced by \cite{RockovaGeorge2018}, forms a continuum between these penalized likelihood and spike-and-slab constructs. The spike-and-slab LASSO methodology has experienced rapid development in recent years, and its scope now extends well beyond the normal linear regression model \eqref{linearregression}. The purpose of this paper is to offer a timely review of the SSL and its many variants. We first provide a basic review of frequentist and Bayesian approaches to high-dimensional variable selection and estimation under the model \eqref{linearregression}. We then review the spike-and-slab LASSO and provide an overview of its attractive properties and techniques to implement it within the context of normal linear regression. Next, we review the methodological developments of the SSL and some of its theoretical developments.

\section{Variable selection in high dimensions: Frequentist and Bayesian strategies} \label{FrequentistVsBayes}

We first review the frequentist penalized regression and the Bayesian spike-and-slab frameworks before showing how the spike-and-slab LASSO bridges the gap between them.

\subsection{Penalized likelihood approaches} \label{PenalizedLikelihood}

In the frequentist high-dimensional literature, there have been a variety of penalized likelihood approaches proposed to estimate $\boldsymbol{\beta}$ in \eqref{linearregression}. A variant of the penalized likelihood approach estimates $\boldsymbol{\beta}$ with
\begin{equation} \label{penalizedlikelihood}
	\widehat{\boldsymbol{\beta}} = \argmax_{\boldsymbol{\beta} \in \mathbb{R}^{p}} - \frac{1}{2} \lVert \boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta} \rVert_2^2 + \textrm{pen}_{\lambda} (\boldsymbol{\beta}),
\end{equation}
where $\textrm{pen}_{\lambda} (\boldsymbol{\beta})$ is a penalty function indexed by penalty parameter $\lambda$. Most of the literature has focused on penalty functions which are separable, i.e. $\textrm{pen}_{\lambda} (\boldsymbol{\beta}) = \sum_{j=1}^{p} \rho_{\lambda} (\beta_j)$. In particular, the popular least absolute shrinkage and selection operator (LASSO) penalty of \cite{Tibshirani1996} uses the function $\rho_{\lambda} (\beta_j) = - \lambda \lvert \beta_j \rvert$. Besides the LASSO and its many variants \cite{BelloniChernozhukovWang2011,SunZhang2012,ZhangZhang2014,Zou2006,ZouHastie2005}, other popular choices for $\rho_{\lambda} (\cdot)$ include non-concave penalty functions, such as the smoothly clipped absolute deviation (SCAD) penalty \cite{FanLi2001} and the minimax concave penalty (MCP) \cite{Zhang2010}. All of the aforementioned penalties threshold some coefficients to zero, thus enabling them to perform variable selection and estimation simultaneously. In addition, SCAD and MCP also mitigate the well-known estimation bias of the LASSO.

Any penalized likelihood estimator \eqref{penalizedlikelihood} also has a Bayesian interpretation in that it can be seen as a posterior mode under an independent product prior $p(\boldsymbol{\beta} \,|\: \lambda) = \prod_{j=1}^{p} p(\beta_j \,|\: \lambda)$, where $\textrm{pen}_{\lambda} (\boldsymbol{\beta}) = \log p (\boldsymbol{\beta} \,|\: \lambda) = \sum_{j=1}^{p} \log p (\beta_j \,|\: \lambda)$. In particular, the solution to the LASSO is equivalent to the posterior mode under a product of Laplace densities indexed by hyperparameter, $\lambda$:
\begin{align} \label{bayesianLASSO}
	p ( \boldsymbol{\beta} \,|\: \lambda) = \prod_{j=1}^{p} \frac{\lambda}{2} e^{- \lambda \lvert \beta_j \rvert }.
\end{align}
This prior, known as the Bayesian LASSO, was first introduced by \cite{ParkCasella2008}. In \cite{ParkCasella2008}, both fully Bayes and empirical Bayes procedures were developed to tune the hyperparameter $\lambda$ in \eqref{bayesianLASSO}. The fully Bayes approach of placing a prior on $\lambda$, in particular, renders the Bayesian LASSO penalty \textit{non}-separable. Thus, the fully Bayesian LASSO has the added advantage of being able to share information across different coordinates. Despite this benefit, \cite{RockovaGeorge2016Abel} showed that the fully Bayesian LASSO cannot simultaneously adapt to sparsity \textit{and} avoid the estimation bias issue of the original LASSO. In addition, \cite{GhoshTangGhoshChakrabarti2016} showed that the univariate Bayesian LASSO often undershrinks negligible coefficients, while overshrinking large coefficients. Finally,  \cite{CastilloSchmidtHieberVanDerVaart2015} also proved that the posterior under the Bayesian LASSO contracts at a suboptimal rate. In Sections \ref{SpikeAndSlabLASSO}-\ref{Illustration}, we will illustrate how the \textit{spike-and-slab LASSO} mitigates these issues.

In addition to the spike-and-slab LASSO, other alternative priors have also been proposed to overcome the limitations of the Bayesian LASSO \eqref{bayesianLASSO}. These priors, known as global-local shrinkage (GL) priors, place greater mass around zero and have heavier tails than the Bayesian LASSO. Thus, GL priors shrink small coefficients more aggressively towards zero, while their heavy tails prevent overshrinkage of large coefficients. Some examples include the normal-gamma prior \cite{GriffinBrown2010}, the horseshoe prior \cite{CarvalhoPolsonScott2010}, the generalized double Pareto prior \cite{ArmaganDunsonLee2013}, the Dirichlet-Laplace prior \cite{BhattacharyaPatiPillaiDunson2015}, and the normal-beta prime prior \cite{BaiGhosh2021}. We refer the reader to \cite{BhadraDattaPolsonWillard2019} for a detailed review of GL priors.

\subsection{Spike-and-slab priors} \label{sspriors}
In the Bayesian framework, variable selection under the linear model \eqref{linearregression} arises directly from probabilistic considerations and has frequently been carried out through placing spike-and-slab priors on the coefficients of interest. The spike-and-slab prior was first introduced by \cite{MitchellBeauchamp1988} and typically has the following form,
\begin{equation} \label{pointmassspikeandslab}
	\begin{array}{rl}
		p (\boldsymbol{\beta} \,|\: \boldsymbol{\gamma}, \sigma^2 ) = & \displaystyle \prod_{j=1}^{p} \left[ (1-\gamma_j) \delta_0 (\beta_j) + \gamma_j p (\beta_j \,|\: \sigma^2 ) \right], \\
		p (\boldsymbol{\gamma} \,|\: \theta) = & \displaystyle \prod_{j=1}^{p} \theta^{\gamma_j} (1-\theta)^{1-\gamma_j}, \hspace{.5cm} \theta \sim p(\theta), \\
		\sigma^2 \sim & p(\sigma^2),
	\end{array}
\end{equation}
where $\delta_0$ is a point mass at zero used to model the negligible entries (the ``spike''), $p( \beta_j \,|\: \sigma^2)$ is a diffuse and/or heavy-tailed density (rescaled by the variance $\sigma^2$) to model the non-negligible entries (the ``slab''), $\boldsymbol{\gamma}$ is a binary vector that indexes the $2^p$ possible models, and $\theta \in (0,1)$ is a mixing proportion. The error variance $\sigma^2$ is typically endowed with a conjugate inverse gamma prior or an improper Jeffreys prior, $p (\sigma^2) \propto \sigma^{-2}$. With a well-chosen prior on $\theta$, this prior \eqref{pointmassspikeandslab} also automatically favors parsimonious models in high dimensions, thus avoiding the curse of dimensionality.

The point-mass spike-and-slab prior \eqref{pointmassspikeandslab} is often considered ``theoretically ideal,'' or a ``gold standard'' for sparse Bayesian problems \cite{CarvalhoPolsonScott2009,PolsonSun2019,Rockova2018}. In high dimensions, however, exploring the full posterior over the entire model space using point-mass spike-and-slab priors \eqref{pointmassspikeandslab} can be computationally prohibitive, in large part because of the combinatorial complexity of updating the discrete indicators $\boldsymbol{\gamma}$. There has been some work to mitigate this issue by using either shotgun stochastic search (SSS) \cite{BottoloRichardson2010,HansDobraWest2007} or variational inference (VI) \cite{RaySzabo2020} to quickly identify regions of high posterior probability.

As an alternative to the point-mass spike-and-slab prior, fully continuous spike-and-slab models have been developed. In these continuous variants, the point-mass $\delta_0$ in \eqref{pointmassspikeandslab} is replaced by a continuous density that is heavily concentrated about zero. The first such continuous relaxation was made by \cite{GeorgeMcCulloch1993}, who used a normal density with very small variance for the spike and a normal density with very large variance for the slab. Specifically, the prior for $\boldsymbol{\beta}$ in \cite{GeorgeMcCulloch1993} is
\begin{equation} \label{spikeandslabnormals}
	p(\boldsymbol{\beta} \,|\: \boldsymbol{\gamma}, \sigma^2) = \prod_{j=1}^{p} \left[ (1-\gamma_j) \mathcal{N}(0, \sigma^2 \tau_0^2) + \gamma_j \mathcal{N}(0, \sigma^2 \tau_1^2) \right],
\end{equation}
where $0 < \tau_0^2 \ll \tau_1^2$. \cite{GeorgeMcCulloch1993} developed a stochastic search variable selection (SSVS) procedure based on posterior sampling with Markov chain Monte Carlo (MCMC) and thresholding the posterior inclusion probabilities, $\Pr (\gamma_j = 1 \,|\: \boldsymbol{y}), j=1, \ldots, p$. In practice, the ``median thresholding'' rule \cite{BarbieriBerger2004}, i.e. $\Pr (\gamma_j = 1 \,|\: \boldsymbol{y}) > 0.5, j=1, \ldots, p$, is often used to perform variable selection.  \cite{IshwaranRao2005, NarisettyHe2014} further extended the model \eqref{spikeandslabnormals} by rescaling the variances $\tau_0^2$ and $\tau_1^2$ with sample size $n$ in order to better control the amount of shrinkage for each individual coefficient.

To further reduce the computational intensiveness of SSVS, a deterministic optimization procedure called EM variable selection (EMVS) was developed by \cite{RockovaGeorge2014}. The EMVS procedure employs \eqref{spikeandslabnormals} as the prior for $\boldsymbol{\beta}$ and uses an EM algorithm to target the posterior mode for $(\boldsymbol{\beta}, \theta, \sigma )$. (\cite{RockovaGeorge2014} also consider continuous spike-and-slab models where the slab, $\mathcal{N}(0, \sigma^2 \tau_1^2)$, is replaced with a polynomial-tailed density, such as a Student's $t$ or a Cauchy distribution, to prevent overshrinkage of the non-negligible entries in $\boldsymbol{\beta}$).   Compared to SSS and SSVS, EMVS has been shown to more rapidly and reliably identify those sets of higher probability submodels which may be of most interest \cite{RockovaGeorge2014}. Recently, \cite{KimGao2019} proposed a general algorithmic framework for Bayesian variable selection with graph-structured sparsity which subsumes the EMVS algorithm as a special case. Like  SSVS, these algorithms also require thresholding the posterior inclusion probabilities to perform variable selection. Letting $(\boldsymbol{\widehat{\beta}}, \widehat{\theta}, \widehat{\sigma})$ denote the posterior mode for $(\boldsymbol{\beta}, \theta, \sigma)$, \cite{RockovaGeorge2014} recommend using median thresholding, $\Pr(\gamma_j = 1 \,|\: \boldsymbol{\widehat{\beta}}, \widehat{\theta}, \widehat{\sigma}) > 0.5$, for selection.

\section{The spike-and-slab LASSO} \label{SpikeAndSlabLASSO}

Having reviewed the penalized likelihood and spike-and-slab paradigms for sparse modeling in the normal linear regression model \eqref{linearregression}, we are now in a position to review the spike-and-slab LASSO (SSL) of \cite{RockovaGeorge2018}. The SSL forms a bridge between these two parallel developments, thereby combining the strengths of both approaches into a single procedure. Throughout this section and Section \ref{Computing}, we assume that $\boldsymbol{y}$ has been centered at zero to avoid the need for an intercept and that the design matrix $\boldsymbol{X}$ has been centered and standardized so that $\lVert \boldsymbol{x}_j \rVert_2^2 = n$ for all $1 \leq j \leq p$.

\subsection{Prior specification}

The spike-and-slab LASSO prior is specified as
\begin{equation} \label{SSLprior}
	\begin{array}{rl}
		p(\boldsymbol{\beta}\,|\:\boldsymbol{\gamma} ) = & \displaystyle \prod_{j=1}^{p} \left[ (1-\gamma_j) \psi(\beta_j \,|\: \lambda_0) + \gamma_j \psi (\beta_j \,|\: \lambda_1) \right], \\
		p(\boldsymbol{\gamma} \,|\: \theta) = & \displaystyle \prod_{j=1}^{p} \left[ \theta^{\gamma_j} (1-\theta)^{1-\gamma_j} \right], \\
		\theta \sim &  \mathcal{B}eta(a,b),
	\end{array}
\end{equation}
where $\psi ( \beta \,|\: \lambda ) = (\lambda / 2) e^{-\lambda \lvert \beta \rvert}$ denotes the Laplace density with scale parameter $\lambda$.   Figure \ref{fig:laplace} depicts the Laplace density for two different choices of scale parameter. We see that for large $\lambda$ ($\lambda=20$), the density is very peaked around zero, while for small $\lambda$ ($\lambda = 1$), it is diffuse. Therefore, in our prior \eqref{SSLprior}, we typically set $\lambda_0 \gg \lambda_1$, so that $\psi ( \cdot \,|\: \lambda_0)$ is the ``spike'' and $\psi ( \cdot\,|\:\lambda_1 )$ is the ``slab.''

\begin{figure}[t!]
	\centering
	\includegraphics[width=0.7\textwidth]{laplace_central_region.pdf}
	\caption{Plot of the central region for the Laplace density with two different choices of scale parameter. }
	\label{fig:laplace}
\end{figure}

The original SSL model of \cite{RockovaGeorge2018} assumed known variance $\sigma^2=1$. \cite{MoranRockovaGeorge2018} extended the SSL to the unknown variance case. As $\sigma^2$ is typically unknown, we consider the hierarchical formulation in \cite{MoranRockovaGeorge2018} in this paper and place an independent Jeffreys prior on $\sigma^2$,
\begin{align*}
	p(\sigma^2) \propto \sigma^{-2}.
\end{align*}
Note that unlike the mixture of normals \eqref{spikeandslabnormals}, we do \textit{not} scale the Laplace priors in $p( \boldsymbol{\beta} \,|\: \gamma)$ by $\sigma^2$. \cite{MoranRockovaGeorge2018} showed that such scaling severely underestimates the variance $\sigma^2$ when $\boldsymbol{\beta}$ is sparse or when $p > n$, thus making the model prone to overfitting.

By choosing $\lambda_1 = \lambda_0$ in \eqref{SSLprior}, we obtain the familiar LASSO $\ell_1$ penalty. On the other hand, if $\lambda_0 \rightarrow \infty$, we obtain the ``theoretically ideal'' point-mass spike-and-slab \eqref{pointmassspikeandslab} as a limiting case. Thus, a feature of the SSL prior is its ability to induce a nonconcave continuum between the penalized likelihood and (point-mass) spike-and-slab constructs.

Since it is a mixture of two Laplace distributions, the SSL prior \eqref{SSLprior} can be seen as a two-group refinement of the LASSO's $\ell_1$ penalty on the coefficients. Thus, the posterior mode for $p(\boldsymbol{\beta} \,|\: \boldsymbol{y})$ under \eqref{SSLprior} is \textit{exactly} sparse and can be used to perform simultaneous variable selection and parameter estimation. This automatic modal thresholding property offers an advantage over previous spike-and-slab formulations \eqref{pointmassspikeandslab}-\eqref{spikeandslabnormals} which do not give exactly sparse estimates of the coefficients and which typically require \textit{post hoc} thresholding of the posterior inclusion probabilities for selection.

It is well-known that the original LASSO \cite{Tibshirani1996} suffers from estimation bias, wherein coefficients with large magnitude are overshrunk. One may wonder what advantages the SSL \eqref{SSLprior} confers over penalized likelihood approaches such as the adaptive LASSO, SCAD, or MCP penalties \cite{FanLi2001, Zhang2010,Zou2006} which are designed to mitigate the bias problem of the LASSO. In what follows, we discuss two major advantages of the SSL. First, we demonstrate that the SSL mixes two LASSO ``bias'' terms \textit{adaptively} in such a way that either a very large amount of shrinkage is applied if $\lvert \beta_j \rvert$ is small or a very small amount of shrinkage is applied if $\lvert \beta_j \rvert$ is large. This is in contrast to the adaptive LASSO \cite{Zou2006} and similar penalties which assign \textit{fixed} coefficient-specific penalties and thus do not gear the coefficient-specific shrinkage towards these extremes. Second, the prior on $\theta$ in \eqref{SSLprior} ultimately renders the coordinates in $\boldsymbol{\beta}$ \textit{dependent} in the marginal prior $p(\boldsymbol{\beta})$ and the SSL penalty \textit{non}-separable. This provides the SSL with the additional ability to borrow information across coordinates and adapt to ensemble information about sparsity.

\subsection{Selective shrinkage and self-adaptivity to sparsity}

As noted in Section \ref{PenalizedLikelihood}, any sparsity-inducing Bayesian prior can be recast {  in the penalized likelihood framework by treating} the logarithm of the marginal prior $\log p(\boldsymbol{\beta})$ as a penalty function. The SSL penalty is defined as
\begin{equation} \label{SSLpenalty}
	\textrm{pen} (\boldsymbol{\beta} ) = \log \left[ \frac{p(\boldsymbol{\beta})}{p(\boldsymbol{0}_p )} \right],
\end{equation}
where the penalty has been centered at $\boldsymbol{0}_p$, the $p$-dimensional zero vector, so that $\textrm{pen}({\boldsymbol{0}_p})= 0$ \cite{RockovaGeorge2018}. Using \eqref{SSLpenalty} and some algebra, the log posterior under the SSL prior (up to an additive constant) can be shown to be
\begin{equation} \label{logposteriorSSL}
	L(\boldsymbol{\beta}, \sigma^2) = -\frac{1}{2 \sigma^2} \lVert \boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta} \rVert_2^2 - (n+2) \log \sigma + \sum_{j=1}^{p} \textrm{pen}(\beta_j \,|\: \theta_j),
\end{equation}
where for $j=1, \ldots, p$,
\begin{equation} \label{singletonpenalty}
	\textrm{pen}(\beta_j \,|\: \theta_j) = - \lambda_1 \lvert \beta_j \rvert + \log[ p_{\theta_j}^{\star} (0 ) / p^{\star}_{\theta_j} (\beta_j ) ],
\end{equation}
with
\begin{equation} \label{conditionalinclusionprob}
	p^{\star}_{\theta_j} ( \beta_j ) = \frac{\theta_j \psi  (\beta_j \,|\: \lambda_1)}{\theta_j \psi (\beta_j \,|\: \lambda_1) + (1-\theta_j) \psi (\beta_j \,|\: \lambda_0)},
\end{equation}
and
\begin{equation} \label{thetaj}
	\theta_j = E[\theta \,|\: \boldsymbol{\beta}_{\setminus j}] = \int \theta p (\theta \,|\: \boldsymbol{\beta}_{\setminus j}) d \theta.
\end{equation}
When $p$ is large, \cite{RockovaGeorge2018} noted that $\theta_j$ is very similar to $E [\theta\,|\: \boldsymbol{\beta}] = \int \theta p (\theta \,|\: \boldsymbol{\beta}) d \theta$ for every $j = 1, \ldots, p$. Thus, for practical purposes, we replace the individual $\theta_j$'s in \eqref{logposteriorSSL}-\eqref{thetaj} with a single $\widehat{\theta} = E [ \theta \,|\: \boldsymbol{\beta}] $ going forward.

The connection between the SSL and penalized likelihood methods is made clearer when considering the derivative of each singleton penalty $\textrm{pen}(\beta_j \,|\: \widehat{\theta })$ in \eqref{singletonpenalty}. This derivative corresponds to an implicit bias term \cite{RockovaGeorge2018} and is given by
\begin{equation} \label{penaltyderivative}
	\frac{ \partial \textrm{pen} ( \beta_j \,|\: \widehat{\theta} )}{\partial \lvert \beta_j \rvert } = - \lambda^{\star}_{\hat{\theta}} ( \beta_j ),
\end{equation}
where
\begin{equation} \label{lambdastar}
	\lambda^{\star}_{\hat{\theta}} (\beta_j) = \lambda_1 p^{\star}_{\hat{\theta}} ( \beta_j ) + \lambda_0 [1 - p^{\star}_{\hat{\theta}} (\beta_j)].
\end{equation}
The Karush-Kuhn-Tucker (KKT) conditions yield the following necessary condition for the global mode $\widehat{\boldsymbol{\beta}}$:
\begin{equation} \label{necessarycond}
	\widehat{\beta}_j = \frac{1}{n} \left[ \lvert z_j \rvert - \sigma^2 \lambda^{\star}_{\hat{\theta}} (\widehat{\beta}_j) \right]_{+} \textrm{sign}(z_j), \hspace{.5cm} j=1, \ldots, p,
\end{equation}
where $z_j = \boldsymbol{x}_j^T (\boldsymbol{y} - \sum_{k \neq j} \boldsymbol{x}_k \widehat{\beta}_k )$. Notice that the condition \eqref{necessarycond} resembles the soft-thresholding operator for the LASSO, except that it contains an adaptive penalty term $\lambda_{\hat{\theta}}^{\star}$ for \textit{each} coefficient. In particular, the quantity \eqref{lambdastar} is a weighted average of the two regularization parameters, $\lambda_1$ and $\lambda_0$, and the weight $p_{\hat{\theta}}^{\star} (\beta_j)$. Thus, \eqref{lambdastar}-\eqref{necessarycond} show that the SSL penalty induces an \textit{adaptive} regularization parameter which applies a different amount of shrinkage to each coefficient, unlike the original LASSO which applies the same shrinkage to every coefficient.

It is worth looking at the term $p_{\hat{\theta}}^{\star} (\beta_j)$ more closely. In light of \eqref{conditionalinclusionprob}, this quantity can be viewed as a conditional probability that $\beta_j$ was drawn from the slab distribution rather than the spike distribution, having seen the regression coefficient $\beta_j$. We have $p_{\hat{\theta}}^{\star} (\beta_j) = \Pr (\gamma_j = 1 \,|\: \beta_j, \hat{\theta})$, where
\begin{equation} \label{pstarexpanded}
	p_{\hat{\theta}}^{\star} (\beta_j) = \frac{1}{1 + \frac{(1-\hat{\theta})}{\hat{\theta}} \frac{\lambda_0}{\lambda_1} \exp \left[ - \lvert \beta_j \rvert ( \lambda_0 - \lambda_1 ) \right]}
\end{equation}
is an \textit{exponentially increasing} function in $\lvert \beta_j \rvert$. From \eqref{pstarexpanded}, we see that the functional $p_{\hat{\theta}}^{\star}$ has a sudden increase from near-zero to near-one. Therefore, $p_{\hat{\theta}}^{\star} (\beta_j)$ gears $\lambda_{\hat{\theta}}^{\star}$ in \eqref{lambdastar} towards the extreme values $\lambda_1$ and $\lambda_0$, depending on the size of $\lvert \beta_j \rvert$. Assuming that $\lambda_1$ is sufficiently small (and hence, the slab $\psi (\beta_j \,|\: \lambda_1)$ is sufficiently diffuse), this allows the large coefficients to escape the overall shrinkage effect, in sharp contrast to the single Laplace distribution \eqref{bayesianLASSO}, where the bias issue remains even if a prior is placed on $\lambda$ \cite{RockovaGeorge2016Abel}.

Apart from its selective shrinkage property, a second key benefit of the SSL model \eqref{SSLprior} is its \textit{self-adaptivity} to the sparsity pattern of the data through the prior on the mixing proportion $\theta$, $p(\theta) \sim \mathcal{B}eta(a,b)$. As mentioned previously, this prior ultimately renders the SSL penalty \textit{non}-separable. Fully separable penalty functions, such as those described in Section \ref{PenalizedLikelihood}, are limited by their inability to adapt to common features across model parameters because they treat these parameters independently. In contrast, treating $\theta$ (the expected proportion for non-negligible coefficients in $\boldsymbol{\beta}$) as random, allows for automatic adaptivity to different levels of sparsity. As shown in \eqref{conditionalinclusionprob}-\eqref{thetaj} (and replacing $\theta_j$ with $\hat{\theta}$ and $\boldsymbol{\beta}_{\setminus j}$ with $\boldsymbol{\beta}$), the mixing weight $p^{\star}_{\hat{\theta}}$ is obtained by averaging $p_{\theta}^{\star} (\cdot)$ over $p(\theta \,|\: \boldsymbol{\beta})$, i.e. $p_{\hat{\theta}}^{\star} (\beta) = \int_{0}^{1} p_{\theta}^{\star} (\beta_j) p(\theta \,|\: \boldsymbol{\beta}) d \theta$. It is through this averaging that the SSL penalty \eqref{SSLpenalty} is given an opportunity to borrow information across coordinates and learn about the underlying level of sparsity in $\boldsymbol{\beta}$.

For the hyperparameters $(a,b)$ in the beta prior on $\theta$, \cite{RockovaGeorge2018} recommended the default choice of $a=1, b=p$. By Lemma 4 of \cite{RockovaGeorge2018}, this choice ensures that $E[\theta \,|\: \widehat{\boldsymbol{\beta}}] \sim \widehat{p}_{\gamma}/p$, where $\widehat{p}_{\gamma}$ is the number of nonzero coefficients in $\boldsymbol{\widehat{\beta}}$. Further, this choice of hyperparameters results in an automatic multiplicity adjustment \cite{ScottBerger2010} and ensures that $\theta$ is small (or that most of the coefficients belong to the spike) with high probability. Thus, the SSL also favors parsimonious models in high dimensions and avoids the curse of dimensionality.

\subsection{The spike-and-slab LASSO in action}

Before delving into the implementation details of the SSL, we perform a small simulation study to illustrate the benefits of the adaptive shrinkage of SSL versus the non-adaptive shrinkage of the LASSO. We simulated data of $n=50$ observations with $p=12$ predictors generated as four independent blocks of highly correlated predictors. More precisely, $n$ rows of our design matrix $\boldsymbol{X}$ were generated independently from a $\mathcal{N}_p (\boldsymbol{0}, \boldsymbol{\Sigma})$ distribution with block diagonal covariance matrix $\boldsymbol{\Sigma} = \textrm{bdiag} ( \widetilde{\boldsymbol{\Sigma}}, \ldots, \widetilde{\boldsymbol{\Sigma}})$, where $\widetilde{\boldsymbol{\Sigma}} = \{\widetilde{\sigma}_{ij} \}_{i,j=1}^{3}, \widetilde{\sigma}_{ij} = 0.9$  if $i \neq j$ and $\widetilde{\sigma}_{ii} = 1$. The response was generated from $\boldsymbol{y} \sim \mathcal{N}_n ( \boldsymbol{X} \boldsymbol{\beta}_0, \boldsymbol{I} )$, with $\boldsymbol{\beta}_0 = (1.3, 0, 0, 1.3, 0, 0, 1.3, 0, 0, 1.3, 0, 0)'$. Note that only $x_1$, $x_4$, $x_7$, and $x_{10}$ are non-null in this true model.

\begin{figure}[t!]
	\centering
	\includegraphics[width=0.49\textwidth]{EXAMPLE_NEW_SSL.pdf}
	\includegraphics[width=0.49\textwidth]{EXAMPLE_NEW_LASSOa.pdf}
	\caption{The coefficient paths of $\widehat{\boldsymbol{\beta}}^{\textrm{SSL}}$ (left panel) and $\widehat{\boldsymbol{\beta}}^{\textrm{LASSO}}$ (right panel) as $\lambda_0$ and $\lambda$ respectively are increased. The connected points off the horizontal axis are the nonzero estimates. The points along the horizontal axis are the zero values where the negligible estimates disappear.}
	\label{fig:coefficientpaths}
\end{figure}

We fit both the SSL and the LASSO of \cite{Tibshirani1996} to this model. Figure \ref{fig:coefficientpaths} displays the coefficient paths for both SSL and LASSO as the spike parameter $\lambda_0$ in the SSL and the regularization parameter $\lambda$ in the LASSO are increased. For the SSL, the spike parameter $\lambda_1 = 0.01$ is fixed throughout. Both the SSL and LASSO begin at $\lambda_0 = \lambda = 0$ with the same 12 (nonzero) ordinary least squares (OLS) estimates for $\boldsymbol{\beta}_0$. However, as $\lambda_0$ increases for the SSL,  the eight smaller OLS estimates are gradually shrunk to zero by the SSL's spike. Meanwhile, the four large estimates are held steady by the SSL's slab, eventually stabilizing at values close to their OLS estimates. The SSL correctly selects the four nonzero coefficients in the true model, demonstrating its self-adaptivity to the true sparsity pattern of the data.

In contrast, Figure \ref{fig:coefficientpaths} also shows that as the LASSO's single penalty parameter $\lambda$ increases, \textit{all} twelve estimates are gradually shrunk to zero. This is because without a slab distribution to help hold the large values steady, the LASSO eventually shrinks all estimates to zero for a large enough $\lambda$. Additionally, due to the order in which the 12 estimates have been thresholded to zero, no value of $\lambda$ yields the correct subset selection $\{ x_1, x_4, x_y, x_{10} \}$. In particular, if $\lambda$ is chosen from cross-validation, the LASSO selects a subset of variables with four false positives.

Our small simulation study illustrates the advantage of the SSL over the LASSO. Specifically, because the LASSO applies the same amount of shrinkage to all regression coefficients, it may estimate a null model if its regularization parameter $\lambda$ is too large. The SSL's two-group refinement of the LASSO penalty helps to mitigate this problem by facilitating selective shrinkage. In Section \ref{Illustration}, we further illustrate the strong empirical performance of the SSL in high-dimensional settings when $p > n$.

\section{Computational details} \label{Computing}

We now turn our attention to implementation of the SSL model \eqref{SSLprior} under the normal linear regression model \eqref{linearregression}. The method described in Section \ref{MAPAlgorithm} is implemented in the publicly available \textsf{R} package \texttt{SSLASSO} \cite{SSLASSOpackage}. However, we also describe an alternative implementation approach in Section \ref{EMAlgorithm}, which is amenable to situations outside of the Gaussian likelihood.

\subsection{Coordinate-wise optimization} \label{MAPAlgorithm}

As mentioned in Section \ref{SpikeAndSlabLASSO} and shown in \eqref{necessarycond}, the (global) posterior mode under the SSL prior \eqref{SSLprior} is exactly sparse, while avoiding the excessive bias issue for large coefficients. Therefore, we can obtain estimates for $\boldsymbol{\beta}$ by targeting the posterior mode.

Marginalizing out $\boldsymbol{\gamma}$ in \eqref{SSLprior} gives the prior for $\boldsymbol{\beta}$ (conditional on $\theta$),
\begin{equation} \label{SSLconditionalontheta}
	p(\boldsymbol{\beta} \,|\: \theta) = \prod_{j=1}^{p} \left[ (1-\theta) \psi (\beta_j \,|\: \lambda_0) + \theta \psi (\beta_j \,|\: \lambda_1) \right].
\end{equation}
Using this reparametrization, \cite{MoranRockovaGeorge2018,RockovaGeorge2018} developed a highly efficient coordinate ascent algorithm to quickly target the mode for $(\boldsymbol{\beta}, \sigma^2)$.

Since the SSL is a non-convex method, the KKT conditions only give a necessary condition \eqref{necessarycond} for $\widehat{\boldsymbol{\beta}}$ to be a global mode, but not a sufficient one. When $p > n$ and $\lambda_0 \gg \lambda_1$, the posterior will typically be multimodal. Nevertheless, it is still possible to obtain a refined characterization of the global mode. Building upon theory developed by \cite{ZhangZhang2012}, \cite{MoranRockovaGeorge2018,RockovaGeorge2018} gave necessary \textit{and} sufficient conditions for $\widehat{\boldsymbol{\beta}}$ to be a \textit{global} mode. By Theorems 3-4 in \cite{RockovaGeorge2018} and Propositions 4-5 in \cite{MoranRockovaGeorge2018}, the global mode under the SSL prior \eqref{SSLprior} is a blend of soft-thresholding \textit{and} hard-thresholding, namely
\begin{equation} \label{globalmode}
	\widehat{\beta}_j = \frac{1}{n} \left[ | z_j | - \sigma^2 \lambda_{\hat{\theta}}^{\star} (\widehat{\beta}_j) \right]_{+} \textrm{sign}(z_j) \mathbb{I}( |z_j| > \Delta ),
\end{equation}
where $z_j = \boldsymbol{x}_j^T (\boldsymbol{y} - \sum_{k \neq j} \boldsymbol{x}_k \widehat{\beta}_k)$ and $\Delta \equiv \inf_{t >0} [nt/2 - \sigma^2 \textrm{pen}(t \,|\:\hat{\theta} ) / t]$. In \cite{MoranRockovaGeorge2018}, an approximation for $\Delta$ is given by
\begin{align*}
	\Delta = \left\{
	\begin{array}{ll}
		\sqrt{2 n \sigma^2 \log [ 1 / p_{\hat{\theta}}^{\star} (0)]} + \sigma^2 \lambda_1\ & \textrm{if } g_{\hat{\theta}}(0) > 0, \\
		\sigma^2 \lambda_{\hat{\theta}}^{\star} (0) & \textrm{otherwise},
	\end{array}
	\right.
\end{align*}
where $g_{\theta}(x) = [ \lambda_{\theta}^{\star}(x) - \lambda_1]^2 + (2n / \sigma^2) \log p_{\theta}^{\star} (x)$. The generalized thresholding operator \eqref{globalmode} allows us to eliminate many suboptimal local modes from consideration through the threshold $\Delta$. This refined characterization also facilitates a highly efficient coordinate ascent algorithm \cite{MazumderFriedmanHastie2011} to find the global mode, which we now detail.

After initializing $(\Delta^{(0)}, \boldsymbol{\beta}^{(0)}, \theta^{(0)}, \sigma^{2(0)})$, the coordinate ascent algorithm iteratively updates these parameters until convergence. The update for the threshold $\Delta$ at the $t^{th}$ iteration is
\begin{align*}
	\Delta^{(t)} = \left\{
	\begin{array}{ll}
		\sqrt{2 n \sigma^{2(t-1)} \log[1 / p_{\hat{\theta}^{(t-1)}}^{\star}(0)]} + \sigma^{2(t-1)} \lambda_1 & \textrm{if } g_{\theta^{(t-1)}} (0) > 0, \\
		\sigma^{2(t-1)} \lambda_{\theta^{(t-1)}}^{\star}(0) & \textrm{otherwise}.
	\end{array}
	\right.
\end{align*}
Next, $\boldsymbol{\beta}$ is updated as
\begin{align*}
	\beta_j^{(t)} \leftarrow \frac{1}{n} \left( |z_j| - \lambda_{\hat{\theta}^{(t-1)}}^{\star} ( \widehat{\beta}_j^{(t-1)} ) \right)_{+} \textrm{sign}(z_j) \mathbb{I}( |z_j| > \Delta^{(t)}).
\end{align*}
Using the approximation for $E[\theta \,|\: \widehat{\boldsymbol{\beta}}]$ in Lemma 4 of \cite{RockovaGeorge2018}, the update for $\hat{\theta}$ is
\begin{align*}
	\hat{\theta}^{(t)} \leftarrow \frac{a + \widehat{p}_{\gamma}^{(t)}}{a+b+p},
\end{align*}
where $\widehat{p}_{\gamma}^{(t)}$ is the number of nonzero entries in $\boldsymbol{\beta}^{(t)}$. Finally, the update for $\sigma^2$ is
\begin{align*}
	\sigma^{2(t)} \leftarrow \frac{ \lVert \boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta}^{(t)} \rVert_2^2}{n+2}.
\end{align*}

\subsection{Dynamic posterior exploration} \label{DynamicPosteriorExploration}

The performance of the SSL model depends on good choices for the hyperparameters $(\lambda_0, \lambda_1)$ in \eqref{SSLprior}. To this end, \cite{RockovaGeorge2018} recommend a ``dynamic posterior exploration'' strategy in which the slab hyperparameter $\lambda_1$ is held fixed at a small value and the spike hyperparameter $\lambda_0$ is gradually increased along a ladder of increasing values, $\{ \lambda_0^{1}, \ldots, \lambda_0^{L} \}$. The algorithm is not very sensitive to the specific choice of $\lambda_1$, provided that the slab is sufficiently diffuse. For each $\lambda_0^{s}$ in the ladder for the spike parameters, we reinitialize $(\Delta^{(0)}, \boldsymbol{\beta}^{(0)}, \theta^{(0)}, \sigma^{2(0)})$ using the MAP estimates for these parameters from the previous spike parameter $\lambda_0^{s-1}$ as a ``warm start.''

This sequential reinitialization strategy allows the SSL to more easily find the global mode. In particular, when $(\lambda_1 - \lambda_0)^2 < 4$ and $\sigma^2$ is fixed, the objective \eqref{logposteriorSSL} is convex. The intuition here is to use the solution to the convex problem as a ``warm'' start for the non-convex problem (when $\lambda_0 \gg \lambda_1$). As we increase $\lambda_0$, the posterior becomes ``spikier,'' with the spikes absorbing more and more of the negligible parameters. Meanwhile, keeping $\lambda_1$ fixed at a small value allows the larger coefficients to escape the pull of the spike. For large enough $\lambda_0$, the algorithm will eventually stabilize so that further increases in $\lambda_0$ do not change the solution. In Section \ref{Illustration}, we illustrate this with plots of the SSL solution paths.

Additionally, as noted by \cite{MoranRockovaGeorge2018}, some care must also be taken when updating $\sigma^2$. When $p > n$ and $\lambda_0 \approx \lambda_1$, the model can become saturated, causing the residual variance to go to zero. To avoid this suboptimal mode at $\sigma^2 = 0$, \cite{MoranRockovaGeorge2018} recommend fixing $\sigma^2$ until the $\lambda_0$ value in the ladder at which the algorithm starts to converge in less than 100 iterations. Then, $\boldsymbol{\beta}$ and $\sigma^2$ are simultaneously updated for the next largest $\lambda_0$ in the sequence. The complete algorithm for coordinate-wise optimization with dynamic posterior exploration is given in Section 4 of the supplementary material in \cite{MoranRockovaGeorge2018}. This algorithm is implemented in the \textsf{R} package \texttt{SSLASSO}.

\subsection{EM implementation of the spike-and-slab LASSO} \label{EMAlgorithm}

The coordinate ascent algorithm of Section \ref{MAPAlgorithm} specifically appeals to the theoretical framework of \cite{ZhangZhang2012} to search for the global SSL mode $\boldsymbol{\widehat{\beta}}$.  An alternative approach, also proposed by  \cite{RockovaGeorge2018},  is to use an EM algorithm in the vein of EMVS \cite{RockovaGeorge2014}.    Again treating the latent variables $\boldsymbol{\gamma}$ in \eqref{SSLprior} are treated as ``missing'' data, this EM implementation of the SSL proceeds as follows.

In the E-step at the $t$th iteration, we compute $E[\tau_j \,|\: \boldsymbol{y}, \boldsymbol{\beta}^{(t-1)}, \theta^{(t-1)}, \sigma^{2(t-1)}] =  p_{\theta^{(t-1)}}^{\star} (\beta_j^{(t-1)})$, where $p_{\theta}^{\star}$ is as in \eqref{conditionalinclusionprob}. The M-step then iterates through the following updates:
\begin{align*}
	& \boldsymbol{\beta}^{(t)} \leftarrow \displaystyle \argmax_{\boldsymbol{\beta} \in \mathbb{R}^{p}} \left\{ - \frac{1}{2} \lVert \boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta} \rVert_2^2 - \sum_{j=1}^{p} \sigma^{2(t-1)} \lambda_{\theta^{(t-1)}}^{\star}(\beta_j^{(t-1)}) \lvert \beta_j \rvert \right\}, \\
	& \theta^{(t)} \leftarrow \frac{ \sum_{j=1}^{p} p_{\theta^{(t-1)}}^{\star} (\beta_j^{(t)}) + a - 1}{a+b+p - 2},
	\\
	& \sigma^{2(t)} \leftarrow \frac{ \lVert \boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta}^{(t)} \rVert_2^2}{n+2},
\end{align*}
where $\lambda_{\theta}^{\star} (\beta) = \lambda_1 p_{\theta}^{\star} (\beta) + \lambda_0 [ 1- p_{\theta}^{\star} (\beta) ]$. Note that the update for $\boldsymbol{\beta}^{(t)}$ is an adaptive LASSO regression with weights $\sigma^2 \lambda_{\theta}^{\star}$ and hence can be solved very efficiently using coordinate descent algorithms \cite{FriedmanHastieTibshirani2010}. Like EMVS \cite{RockovaGeorge2014}, the dynamic posterior exploration strategy detailed in Section \ref{DynamicPosteriorExploration} can be used to find a more optimal mode for $(\boldsymbol{\beta}, \sigma^2)$.

This EM approach can be straightforwardly adapted for other statistical models where the SSL prior \eqref{SSLprior} is used (such as the methods described in Section \ref{MethodologicalExtensions}) but where the likelihood function differs and the theory of \cite{ZhangZhang2012} is not applicable.
Similar to the coordinate ascent algorithm described in Section \ref{MAPAlgorithm}, this EM algorithm may be sensitive to the initialization of $(\boldsymbol{\beta}^{(0)}, \theta^{(0)}, \sigma^{2(0)} )$. The dynamic posterior exploration strategy described earlier can partly help to mitigate this issue, since the posterior starts out relatively flat when $\lambda_0 \approx \lambda_1$ but becomes ``spikier'' as $\lambda_0$ increases. By the time that the spikes have reappeared, the ``warm start'' solution from the previous $\lambda_0$ in the ladder should hopefully be in the basin of dominant mode. Other strategies such as running the algorithm for a wide choice of starting values or deterministic annealing can also aid in adding robustness against poor initializations \cite{McLachlanBasford1988,RockovaGeorge2014,UedaNakano1998}.

\section{Uncertainty quantification} \label{UncertaintyQuantification}

While the algorithms described in Section \ref{Computing} can be used to rapidly target the modes of the SSL posterior, providing a measure of uncertainty for our estimates is a challenging task. In this section, we outline two possible strategies for the task of uncertainty quantification. The first is based on debiasing the posterior mode. The second involves posterior simulation.

\subsection{Debiasing the posterior mode}

One possible avenue for uncertainty quantification is to use debiasing \cite{BaiMoranAntonelliChenBoland2019, JavanmardMontanari2018, VanDeGeerBuhlmannRitovDezeure2014,ZhangZhang2014}. Let $\widehat{\boldsymbol{\Sigma}} = \boldsymbol{X}^T \boldsymbol{X} / n$ and let $\widehat{\boldsymbol{\Theta}}$ be an approximate inverse of $\widehat{\boldsymbol{\Sigma}}$. Note that when $p > n$, $\boldsymbol{X}$ is singular, so $\widehat{\boldsymbol{\Sigma}}^{-1}$ does not necessarily exist. However, we can still obtain a sparse estimate of the precision matrix $\widehat{\boldsymbol{\Theta}}$ for the rows of $\boldsymbol{X}$ by using techniques from the graphical models literature, e.g. the nodewise regression procedure in \cite{MeinshausenBuhlmann2006} or the graphical lasso \cite{FriedmanHastieTibshirani2007}. We define the quantity $\widehat{\boldsymbol{\beta}}_d$ as
\begin{equation}
	\widehat{\boldsymbol{\beta}}_d = \widehat{\boldsymbol{\beta}} + \widehat{\boldsymbol{\Theta}} \boldsymbol{X}^T (\boldsymbol{y} - \boldsymbol{X} \widehat{\boldsymbol{\beta}})/n.
\end{equation}
where $\widehat{\boldsymbol{\beta}}$ is the MAP estimator of $\boldsymbol{\beta}$ under the SSL model. By \cite{VanDeGeerBuhlmannRitovDezeure2014}, this quantity $\widehat{\boldsymbol{\beta}}_d$ has the following asymptotic distribution:
\begin{equation} \label{asymptoticdist}
	\sqrt{n}(\widehat{\boldsymbol{\beta}}_d - \boldsymbol{\beta}) \sim \mathcal{N}(\boldsymbol{0}, \sigma^2 \widehat{\boldsymbol{\Theta}} \widehat{\boldsymbol{\Sigma}} \widehat{\boldsymbol{\Theta}}^T).
\end{equation}
For inference, we replace the population variance $\sigma^2$ in \eqref{asymptoticdist} with the modal estimate $\widehat{\sigma}^2$ from the SSL model. Let $\widehat{\beta}_{dj}$ denote the $j$th coordinate of $\widehat{\boldsymbol{\beta}}_d$. We have from \eqref{asymptoticdist} that the $100(1-\alpha)  \%$ asymptotic pointwise confidence intervals for $\beta_{j}, j = 1, \ldots, p$, are
\begin{align} \label{confidenceintervals}
	[ \widehat{\beta}_{dj} - c(\alpha, n, \widehat{\sigma}^2), \widehat{\beta}_{dj} + c(\alpha, n, \widehat{\sigma}^2) ],
\end{align}
where $c(\alpha, n, \widehat{\sigma}^2) := \Phi^{-1} (1-\alpha/2) \sqrt{ \widehat{\sigma}^2 ( \widehat{\boldsymbol{\Theta}} \widehat{\boldsymbol{\Sigma}} \widehat{\boldsymbol{\Theta}}^T )_{jj} / n}$ and $\Phi(\cdot)$ denotes the cumulative distribution function of $\mathcal{N}(0,1)$.

Note that the posterior modal estimate $\widehat{\boldsymbol{\beta}}$ under the SSL prior already has much less bias than the LASSO estimator \cite{Tibshirani1996}. Therefore, the purpose of the debiasing procedure above is mainly to obtain an estimator with an asymptotically normal distribution from which we can construct asymptotic pointwise confidence intervals. While this procedure is asymptotically valid, \cite{AntonelliParmigianiDominici2019} showed through numerical studies that constructing confidence intervals based on asymptotic arguments may provide coverage below the nominal level in finite samples, especially small samples. Therefore, it may be more ideal to use the actual SSL posterior $p(\boldsymbol{\beta} \,|\: \boldsymbol{y})$ for inference.

\subsection{Posterior sampling for the spike-and-slab LASSO}
Fully Bayesian inference with the SSL can be carried out via posterior simulation.
However, posterior sampling under spike-and-slab priors has continued to pose challenges.
One immediate strategy for sampling from the SSL posterior is the SSVS algorithm of \cite{GeorgeMcCulloch1993}, described in Section \ref{sspriors}. One can regard the Laplace distribution  as a scale mixture of Gaussians with an exponential mixing distribution \cite{ParkCasella2008} and perform a variant of SSVS. Recently, several clever computational tricks have been suggested that avoid costly matrix inversions needed by SSVS by using linear solvers \cite{BhattacharyaChakrabortyMallick2016}, low-rank approximations \cite{JohndrowOrensteinBhattacharya2020}, or by disregarding correlations between active and inactive coefficients \cite{NarisettyShenHe2019}. These techniques can be suitably adapted for fast posterior sampling of the SSL as well.

Intrigued by the speed of SSL mode detection, \cite{NieRockova2020} explored the possibility of turning SSL into approximate posterior sampling by performing MAP optimization on many independently perturbed datasets.  Building on Bayesian bootstrap ideas, they introduced a method for approximate sampling called Bayesian bootstrap spike-and-slab LASSO (BB-SSL) which scales linearly with both $n$ and $p$.  Beyond its scalability, they show that BB-SSL has strong theoretical support, matching the convergence rate of the original posterior in sparse normal-means and in high-dimensional regression.

\section{Illustrations} \label{Illustration}

In this section, we illustrate the SSL's potential for estimation, variable selection, and prediction on both simulated and real high-dimensional data sets.

\subsection{Example on synthetic data} \label{Simulations}

For our simulation study, we slightly modified the settings in \cite{MoranRockovaGeorge2018}. We set $n = 100$ and $p = 1000$ in \eqref{linearregression}. The design matrix $\boldsymbol{X}$ was generated from a multivariate Gaussian distribution with mean $\boldsymbol{0}_p$ and a block-diagonal covariance matrix $\boldsymbol{\Sigma} = \textrm{bdiag} ( \widetilde{\Sigma}, \ldots, \widetilde{\Sigma} )$, where $\widetilde{\Sigma} = \{ \widetilde{\sigma} \}_{i,j=1}^{50}$, with $\widetilde{\sigma}_{ij} = 0.9$ if $i \neq j$ and $\widetilde{\sigma}_{ii} = 1$. The true vector of regression coefficients $\boldsymbol{\beta}_0$ was constructed by assigning regression coefficients $\{-3.5, -2.5, -1.5, 1.5, 2.5, 3.5 \}$ to 6 entries located at the indices $\{ 1, 51, 101, 151, 201, 251 \}$ and setting the remaining coefficients equal to zero. Hence, there were 20 independent blocks of 50 highly correlated predictors, where the first six blocks contained only one active predictor. We then generated the response $\boldsymbol{y}$ using \eqref{linearregression}, where the error variance was set as $\sigma^2 = 3$.

We compared the SSL with the LASSO \cite{Tibshirani1996}, SCAD \cite{FanLi2001}, and MCP \cite{Zhang2010}. The SSL method was applied using the \textsf{R} package \texttt{SSLASSO}. The competing methods were applied using the \textsf{R} package \texttt{ncvreg}. We repeated our experiment 500 times with new covariates and responses generated each time. For each experiment, we recorded the mean squared error (MSE) and mean prediction error (MPE), defined as
\begin{align*}
	\textrm{MSE} = \frac{1}{p} \lVert \widehat{\boldsymbol{\beta}} - \boldsymbol{\beta}_0 \rVert_2^2 \hspace{.3cm} \textrm{ and } \hspace{.3cm}
	\textrm{MPE} = \frac{1}{n} \lVert \boldsymbol{X} ( \widehat{\boldsymbol{\beta}} - \boldsymbol{\beta}_0 ) \rVert_2^2.
\end{align*}
We also kept track of $\widehat{p}_{\gamma}$, or the size of the model selected by each of these methods. Finally, we recorded the false discovery rate (FDR), the false negative rate (FNR), and the Matthews correlation coefficient (MCC) \cite{Matthews1975}, defined respectively as
\begin{align*}
	\textrm{FDR} & = \frac{\textrm{FP}}{\textrm{TN}+\textrm{FP}}, \hspace{.5cm} \textrm{FNR} = \frac{\textrm{FN}}{\textrm{TP}+\textrm{FN}}, \\
	\textrm{MCC} & = \frac{ \textrm{TP} \times \textrm{TN} - \textrm{FP} \times \textrm{FN}}{ \sqrt{ ( \textrm{TP} + \textrm{FP} )(\textrm{TP} + \textrm{FN})(\textrm{TN} + \textrm{FP})(\textrm{TN} + \textrm{FN})}},
\end{align*}
where TP, TN, FP, and FN denote the number of true positives, true negatives, false positives, and false negatives respectively. The MCC is a correlation coefficient between the predicted set of significant coefficients and the actual set of nonzero coefficients \cite{Matthews1975}. MCC has a range of -1 to 1,  with -1 indicating completely incorrect selection (i.e. TP=TN=0) and 1 indicating completely correct variable selection (i.e. FP=FN=0). Models with MCC closer to 1 have higher selection accuracy. \textsf{R} code to reproduce these experiments is available in the online supplementary material.

\begin{table}[t!]
	\centering
	\begin{tabularx}{\textwidth}{@{}*7{>{\raggedright\arraybackslash}X}@{}}
		\hline
		& MSE & MPE & $\widehat{p}_{\gamma}$ & FDR & FNR & MCC \\
		\hline
		\hline
		SSL & \textbf{0.0067} (0.0076) & \textbf{0.701} (0.542) & \textbf{6.05} (0.271) & \textbf{0.0012} (0.0010) & 0.187 (0.160) & \textbf{0.809} (0.162) \\
		\hline
		LASSO & 0.011 (0.0045) & 1.14 (0.303) & 33.38 (5.38) & 0.028 (0.0055) & \textbf{0.083} (0.109) & 0.387 (0.062) \\
		\hline
		SCAD & 0.011 (0.012) & 0.985 (0.691) & 12.74 (3.98) & 0.0081 (0.0043) & 0.225 (0.187) & 0.554 (0.178) \\
		\hline
		MCP & 0.020 (0.016) & 1.55 (0.849) & 11.31 (2.70) & 0.0077 (0.0031) & 0.395 (0.211) & 0.447 (0.173) \\
		\hline
	\end{tabularx}
	\caption{MPE, MPE, estimated model size, FDR, FNR, and MCC for SSL, LASSO, SCAD, and MCP.  The results are averaged across 500 replications. In parentheses, we report the empirical standard errors.}
	\label{Table:1}
\end{table}

Table \ref{Table:1} reports our results averaged across the 500 replications. We see that the SSL had the lowest average MSE and MPE, in addition to selecting (on average) the most parsimonious model. The LASSO (along with SCAD) had the second lowest MSE, but it tended to select far more variables than the other methods, leading to the highest FDR. In contrast, SSL had the lowest FDR and the highest MCC, indicating that the SSL had the best overall variable selection performance of all the methods. Our simulation study demonstrates that SSL achieves both parsimony \textit{and} accuracy of estimation and selection.

Figure \ref{fig:dynamicposteriorexploration} illustrates the benefits of the dynamic posterior exploration approach outlined in Section \ref{DynamicPosteriorExploration}. Specifically, Figure \ref{fig:dynamicposteriorexploration} plots the solution paths for the regression coefficients from one of our experiments as the spike hyperparameter $\lambda_0$ increases. We see that the SSL solution stabilizes fairly quickly (when $\lambda_0$ is less than 20), so that further increases in $\lambda_0$ do not change the solution. This demonstrates that dynamic posterior exploration offers a viable alternative to cross-validation. The \textsf{R} package \texttt{SSLASSO} provides the functionality to generate plots of these solution paths.

\begin{figure}[t!]
	\centering
	\includegraphics[width=\textwidth]{SSL_solution_paths.pdf}
	\caption{The solution paths for the SSL along the ladder of spike parameters $\lambda_0$. We see that the solution stabilizes after a certain point, so that further increases in $\lambda_0$ do not change the solution.  The points along the horizontal axis are the zero values where the negligible estimates disappear. }
	\label{fig:dynamicposteriorexploration}
\end{figure}

\subsection{Bardet-Beidl syndrome gene expression study} \label{RealDataAnalysis}

We now analyze a microarray data set consisting of gene expression measurements from the eye tissue of 120 laboratory rats\footnote{Data accessed from the Gene Expression Omnibus \texttt{ www.ncbi.nlm.nih.gov/geo} (accession no. GSE5680).}. The data was originally studied by \cite{Scheetz06} to investigate mammalian eye disease. In this data, the goal is to identify genes which are associated with the gene TRIM32. TRIM32 has previously been shown to cause Bardet-Biedl syndrome \cite{chiang06}, a disease affecting multiple organs including the retina.

The original data consists of 31,099 probe sets. For our analysis, we included only the 10,000 probe sets with the largest variances in expression (on the log scale). This resulted in $n = 120$ and $p = 10,000$. We then fit the model \eqref{linearregression} with an SSL penalty.  We compared the SSL approach to LASSO, SCAD, and MCP.

To assess predictive accuracy, we randomly split the data set into 90 training observations and 30 test observations. We then fit the models on the training set and used the estimated  $\widehat{\boldsymbol{\beta}}_{\textrm{train}}$ to compute the mean squared prediction error (MSPE) on the left-out test set,
\begin{equation*}
	\textrm{MSPE} = \frac{1}{30} \sum_{i=1}^{30} (y_{i, \textrm{test}} - \boldsymbol{x}_{i, \textrm{test}}^T \widehat{\boldsymbol{\beta}}_{\textrm{train}})^2,
\end{equation*}
where $(\boldsymbol{x}_{i, \textrm{test}}, y_{i, \textrm{test}}), i = 1, \ldots, 30,$ are the observations in the test set. We repeated this process 100 times and took the average MSPE.

Table \ref{Table:2} shows the results for our analysis, as well as the number of selected probe sets when we fit the different models to the complete data set. SSL had the lowest out-of-sample MSPE, indicating the highest predictive power. MCP selected the most parsimonious model, with only nine probe sets out of the 10,000 selected. However, MCP also had a much higher MSPE than the other methods. On the other hand, SSL selected 28 probe sets (compared to 32 and 44 for LASSO and SCAD respectively) and still achieved the lowest MSPE. Our analysis illustrates that on this particular data set, SSL achieved both the best predictive performance and parsimony.

\begin{table}[t!]
	\centering
	\begin{tabular}{l c c}
		\hline
		& MSPE & Number of selected probe sets \\
		\hline
		SSL & \textbf{0.011} & 28 \\
		LASSO & 0.012 & 32 \\
		SCAD & 0.015 & 44 \\
		MCP & 3.699 &  9 \\
		\hline
	\end{tabular}
	\caption{Average MSPE and the number of selected probe sets for the Bardet-Beidl Syndrome data analysis.}
	\label{Table:2}
\end{table}

Of the 28 probe sets selected by SSL as being significantly associated with TRIM32, 14 of them had identifiable gene symbols. These genes were SCGB1A1, CELF1, ASXL3, FGFR2, MOBP, TGM7, SLC39A6, DDX58, TFF2, CLOCK, DUS4L, HTR5B, BIK, and SLC16A6.  In particular, according to \texttt{https://www.genecards.org} \footnote{Accessed from \texttt{https://www.genecards.org/cgi-bin/carddisp.pl?gene=TRIM32} on October 11, 2020.}, SCGB1A1 is known to be an interacting protein for the TRIM32 gene. The other associations that we found may be useful for researchers in studying the genetic factors contributing to Bardet-Biedl syndrome.

\section{Methodological extensions} \label{MethodologicalExtensions}

While we have focused on the normal linear regression model \eqref{linearregression} in Sections \ref{FrequentistVsBayes}-\ref{Illustration}, the spike-and-slab LASSO methodology has now been adopted in a variety of other statistical applications. In this section, we survey some of the extensions of the SSL to models beyond the normal linear regression framework.

\textbf{Generalized linear models (GLMs).} GLMs allow for a flexible generalization of the normal linear regression model \eqref{linearregression} which can accommodate categorical and count data, in addition to continuous variables. Letting $\boldsymbol{x}_i \in \mathbb{R}^{p}$ denote a vector of covariates for the $i$th observation, GLMs assume that the mean of the response variable is related to the linear predictor via a link function,
\begin{equation} \label{linkfunction}
	E(y_i \,|\: \boldsymbol{x}_i) = h^{-1} ( \boldsymbol{x}_i^T \boldsymbol{\beta} ),
\end{equation}
and that the data distribution is expressed as
\begin{equation} \label{GLMlikelihood}
	p( \boldsymbol{y} \,|\: \boldsymbol{X}, \boldsymbol{\beta}, \varphi ) = \prod_{i=1}^{n} p(y_i \,|\: \boldsymbol{x}_i, \boldsymbol{\beta}, \varphi),
\end{equation}
where $\varphi$ is a dispersion parameter and the distribution $p(y_i \,|\: \boldsymbol{x}_i, \boldsymbol{\beta}, \varphi )$ can take various forms, including normal, binomial and Poisson distributions. Obviously, the normal linear model \eqref{linearregression} is a special case of \eqref{linkfunction} with the identity link function $h(u)=u$.

\cite{TangShenZhangYi2017GLM} extended the SSL \eqref{SpikeAndSlabLASSO} to GLMs, including binary regression and Poisson regression, by placing the SSL prior \eqref{SSLprior} on the coefficients vector $\boldsymbol{\beta} \in \mathbb{R}^{p}$ in \eqref{linkfunction} and developing an EM algorithm to perform MAP estimation for $\boldsymbol{\beta}$. For inference with grouped variables in GLMs, \cite{TangShenLiZhangWenQianZhuangShiYi2018} further employed group-specific sparsity parameters $\theta_g$ for each group of variables, instead of a single $\theta$, as in \eqref{SpikeAndSlabLASSO}.

\textbf{Survival analysis.} The SSL has also proven to be useful for predicting censored survival outcomes and detecting and estimating the effects of relevant covariates. Cox proportional hazards models are the most widely used method for studying the relationship between a censored survival response and an explanatory variable $\boldsymbol{x}_i \in \mathbb{R}^{p}$ \cite{KleinMoeschberger2003}. This model assumes that the hazard function of survival time $t$ takes the form,
\begin{equation} \label{CoxRegression}
	h(t \,|\: \boldsymbol{x}_i ) = h_0(t) \exp(\boldsymbol{x}_i^T \boldsymbol{\beta}).
\end{equation}
\cite{TangShenZhangYi2017Cox} introduced the spike-and-slab LASSO Cox model which endows the coefficients $\boldsymbol{\beta}$ in \eqref{CoxRegression} with the SSL prior \eqref{SSLprior}. They developed an EM coordinate ascent algorithm to fit SSL Cox models. \cite{TangLeiZhangYiGuoChenShenYi2019} further introduced the GssLASSO Cox model which incorporates grouping information by endowing each group of coefficients with a group-specific sparsity parameter $\theta_g$ instead of the single $\theta$ of \eqref{SSLprior}.

\textbf{Sparse factor analysis and biclustering.} Factor models aim to explain the dependence structure among high-dimensional observations through a sparse decomposition of a $p \times p$ covariance matrix $\boldsymbol{\Omega}$ as $\boldsymbol{B} \boldsymbol{B}^T + \boldsymbol{\Sigma}$ where $\boldsymbol{B}$ is a $p \times K$ factor loadings matrix with $K \ll p$ and $\boldsymbol{\Sigma} = \textrm{diag}(\sigma_1^2, \ldots, \sigma_p^2)$. A generic latent factor model is
\begin{equation} \label{factormodel}
	\boldsymbol{y}_i = \boldsymbol{B} \boldsymbol{\eta}_{i} + \boldsymbol{\varepsilon}_i, \hspace{.5cm} \boldsymbol{\varepsilon}_i \sim \mathcal{N}_p (\boldsymbol{0}, \boldsymbol{\Sigma}),
\end{equation}
where $\boldsymbol{y}_i$ is a $p$-dimensional continuous response and $\boldsymbol{\eta}_i \sim \mathcal{N}_K (\boldsymbol{0}, \boldsymbol{I}_K)$ are unobserved latent factors. Many existing factor analysis approaches entail prespecification of the unknown factor cardinality $K$ and \textit{post hoc} rotations of the original solution to sparsity.  For the factor model \eqref{factormodel}, \cite{RockovaGeorge2016} endowed the entries of the loading matrix $\boldsymbol{B}$ with independent SSL priors,
\begin{align*}
	p(\beta_{jk} \,|\: \gamma_{jk}, \lambda_{0k}, \lambda_1) = (1-\gamma_{jk}) \psi(\beta_{jk} \,|\: \lambda_{0k} ) + \gamma_{jk}\psi(\beta_{jk}\,|\: \lambda_{1}).
\end{align*}
However, instead of endowing each of the indicators $\gamma_{jk}$ with the usual beta-Bernoulli prior as in \eqref{SSLprior}, \cite{RockovaGeorge2016} endowed these with the Indian buffet process (IBP) prior \cite{GriffithsGhahramani2011}, which avoids the need to prespecify $K$. Further, \cite{RockovaGeorge2016} developed a parameter-expanded EM (PXL-EM) algorithm which employs \textit{intermediate} orthogonal rotations rather than post hoc rotations. In addition to obtaining a sparse solution, the PXL-EM algorithm also converges much faster than the vanilla EM algorithm and offers robustness against poor initializations.

For the problem of biclustering, i.e. identifying clusters using only subsets of their associated features, \cite{MoranRockovaGeorge2019} utilized the factor model \eqref{factormodel} in which \textit{both} the factors $\boldsymbol{\eta} = [ \boldsymbol{\eta}_1^T, \ldots, \boldsymbol{\eta}_n^T ] \in \mathbb{R}^{n \times K}$ and the loadings are sparse. To achieve a doubly sparse representation, \cite{MoranRockovaGeorge2019} placed an SSL prior coupled with an IBP prior on the factors and an SSL prior coupled  with a beta-Bernoulli prior on the loadings. An EM algorithm with a variational step was developed to implement spike-and-slab LASSO biclustering.

\textbf{Graphical models.} Suppose we are given data $\boldsymbol{Y} = (\boldsymbol{y}_1, \ldots, \boldsymbol{y}_n)^T$, where the $\boldsymbol{y}_i$'s are assumed to be iid $p$-variate random vectors distributed as $\mathcal{N}_p (\boldsymbol{0}, \boldsymbol{\Omega}^{-1})$ and $p>n$. In this setting, off-diagonal zero entries $\omega_{ij}$ encode conditional independence between variables $i$ and $j$.	To obtain a sparse estimate of $\boldsymbol{\Omega} = ( \omega_{i,j} )_{i,j}$, \cite{GanNarisettyLiang2018} introduced the following prior on $\boldsymbol{\Omega}$:
\begin{align} \label{graphicalmodelprior}
	p(\boldsymbol{\Omega}) = \prod_{i < j} \left[ (1-\theta) \psi(\omega_{ij}\,|\: \lambda_0 ) + \theta \psi(\omega_{ij} \,|\: \lambda_1 ) \right] \prod_{i=1}^{p} [ \tau e^{-\tau \omega_{ii}} ] \mathbb{I}( \Omega \succ 0 ) \mathbb{I}( \lVert \boldsymbol{\Omega} \rVert_2 \leq B ),
\end{align}
for some $\tau > 0, B > 0$. Here, $\boldsymbol{\Omega} \succ 0$ denotes that $\boldsymbol{\Omega}$ is positive-definite and $\lVert \boldsymbol{\Omega} \rVert_2$ denotes the spectral norm of $\boldsymbol{\Omega}$. The prior on $\boldsymbol{\Omega}$ \eqref{graphicalmodelprior} entails independent exponential priors on the diagonal entries and SSL priors on the off-diagonal entries. A similar prior formulation was considered in \cite{DeshpandeRockovaGeorge2019}, except \cite{DeshpandeRockovaGeorge2019} did not constrain $\boldsymbol{\Omega}$ to lie in the space of $p \times p$ matrices with uniformly bounded spectral norm.  \cite{GanNarisettyLiang2018} showed that constraining the parameter space for $\boldsymbol{\Omega}$ in such a way ensures that a) the corresponding optimization problem for the posterior mode is \textit{strictly} convex, and b) the posterior mode is a symmetric positive definite matrix. \cite{GanNarisettyLiang2018} developed an EM algorithm to estimate the posterior mode of $p(\boldsymbol{\Omega} \,|\: \boldsymbol{Y})$.

\hspace{.2cm} The SSL prior was also extended to perform \textit{joint} estimation of multiple related Gaussian graphical models by \cite{LiMcCormickClark2019}. \cite{LiMcCormickClark2019}
leveraged similarities in the underlying sparse precision matrices and developed an EM algorithm to perform this joint estimation.

\textbf{Seemingly unrelated regression models.} In seemingly unrelated regression models, multiple correlated responses are regressed on multiple predictors. The multivariate linear regression is an important case. Letting $\boldsymbol{y}_i \in \mathbb{R}^{q}$ be the vector of $q$ responses and $\boldsymbol{x}_i \in \mathbb{R}^{p}$ be the vector of $p$ covariates, this model is
\begin{equation} \label{multivariateregression}
	\boldsymbol{y}_i = \boldsymbol{x}_i^T \boldsymbol{B} + \boldsymbol{\varepsilon}_i, \hspace{.5cm} \boldsymbol{\varepsilon}_i \sim \mathcal{N}_q ( \boldsymbol{0}, \boldsymbol{\Omega}^{-1} ),
\end{equation}
\cite{DeshpandeRockovaGeorge2019} introduced the \textit{multivariate spike-and-slab LASSO} (mSSL) to perform joint selection and estimation from the $p \times q$ matrix of regressors $\boldsymbol{B}$ \textit{and} the precision matrix $\boldsymbol{\Omega}$. To obtain a sparse estimate of $(\boldsymbol{B}, \boldsymbol{\Omega})$, \cite{DeshpandeRockovaGeorge2019} placed the SSL prior \eqref{SSLprior} on the individual entries $\beta_{jk}, 1 \leq j \leq p, 1 \leq k \leq q$ in $\boldsymbol{B}$ and a product prior similar to \eqref{graphicalmodelprior} on $\boldsymbol{\Omega}$ (except \cite{DeshpandeRockovaGeorge2019} did not constrain $\boldsymbol{\Omega}$ to have bounded spectral norm). An expectation/conditional maximization (ECM) algorithm was developed to perform this joint estimation.

\textbf{Causal inference.} In observational studies, we are often interested in estimating the causal effect of a treatment $T$ on an outcome $y$, which requires proper adjustment of a set of potential confounders $\boldsymbol{x} \in \mathbb{R}^{p}$. When $p > n$, direct control for all potential confounders is infeasible and standard methods such as propensity scoring \cite{RosenbaumRubin1983} often fail. In this case, it is crucial to impose a low-dimensional structure on the confounder space. Given data $(y_i, T_i, \boldsymbol{x}_i), i=1, \ldots, n,$ where $T_i$ is the treatment effect,  \cite{AntonelliParmigianiDominici2019} estimated the (homogeneous) average treatment effect (ATE) $\Delta(t_1, t_2 ) = E(Y(t_1)-Y(t_2))$ by utilizing the model,
\begin{equation} \label{homogeneoustreatmentefectmodel}
	y_i \,|\: T_i, \boldsymbol{x}_i, \beta_0, \beta_t, \boldsymbol{\beta}, \sigma^2 \sim \mathcal{N}(0, \beta_0 + \beta_t T_i+ \boldsymbol{x}_i ^T \boldsymbol{\beta}, \sigma^2 ).
\end{equation}
Under \eqref{homogeneoustreatmentefectmodel}, the ATE is straightforwardly estimated as $\Delta (t_1, t_2) = (t_1 - t_2) \widehat{\beta}_t$. \cite{AntonelliParmigianiDominici2019} endowed the coefficients of the confounders $\boldsymbol{\beta}$ with the SSL prior \eqref{SSLprior}. In addition, \cite{AntonelliParmigianiDominici2019} also weighted the sparsity parameter $\theta$ in \eqref{SSLprior} by raising $\theta$ to a power $w_j, j=1, \ldots, p$, for each covariate, in order to better prioritize variables belonging to the slab (i.e. $\gamma_j = 1$) if they are also associated with the treatment. \cite{AntonelliParmigianiDominici2019} further extended the model \eqref{homogeneoustreatmentefectmodel} to the more general case of heterogeneous treatment effects.

\textbf{Regression with grouped variables.} Group structure arises in many statistical applications. For example, in genomics, genes within the same pathway may form a group and act in tandem to regulate a biological system. For regression with grouped variables, we can model the response $\boldsymbol{y}$ as
\begin{equation} \label{groupmodel}
	\boldsymbol{y} = \displaystyle \sum_{g=1}^{G} \boldsymbol{X}_g \boldsymbol{\beta}_g + \boldsymbol{\varepsilon}, \hspace{.5cm} \boldsymbol{\varepsilon} \sim \mathcal{N}_n ( \boldsymbol{0}, \sigma^2 \boldsymbol{I}_n),
\end{equation}
where $\boldsymbol{\beta}_g \in \mathbb{R}^{m_g}$ is a coefficients \textit{vector} of length $m_g$, and $\boldsymbol{X}_g$ is an $n \times m_g$ covariate matrix corresponding to group $g = 1, \ldots G$. Under model \eqref{groupmodel}, it is often of practical interest to select non-negligible groups and estimate their effects. To this end, \cite{BaiMoranAntonelliChenBoland2019} introduced the \textit{spike-and-slab group lasso} (SSGL). To regularize groups of coefficients, the SSGL replaces the univariate Laplace densities in the univariate SSL \eqref{SSLprior} with \textit{multivariate} Laplace densities. The SSGL prior is
\begin{equation} \label{ssgroupLASSO}
	\begin{array}{rl}
		p ( \boldsymbol{\beta} \,|\: \theta ) = & \displaystyle \prod_{g=1}^{G} \left[ (1- \theta) \boldsymbol{\Psi} ( \boldsymbol{\beta}_g \,|\: \lambda_0 ) + \theta \boldsymbol{\Psi} ( \boldsymbol{\beta}_g \,|\: \lambda_1 ) \right],  \\
		\theta \sim & \mathcal{B}eta(a,b),
	\end{array}
\end{equation}
where $\boldsymbol{\Psi} (\boldsymbol{\beta}_g \,|\: \lambda ) \propto \lambda^{m_g} e^{- \lambda \lVert \boldsymbol{\beta}_g \rVert_2 }$ and $\lambda_0 \gg \lambda_1$. The SSGL \eqref{ssgroupLASSO} is a two-group refinement of an $\ell_2$ penalty on groups of coefficients. Accordingly, the posterior mode under the SSGL thresholds entire groups of coefficients to zero, while simultaneously estimating the effects of nonzero groups and circumventing the estimation bias of the original group lasso \cite{YuanLin2006}. \cite{BaiMoranAntonelliChenBoland2019} developed an efficient blockwise-coordinate ascent algorithm to implement the SSGL model.

\textbf{Nonparameteric additive regression.} The advent of the SSGL prior \eqref{ssgroupLASSO} paved the way for the spike-and-slab lasso methodology to be extended to nonparametric problems. \cite{BaiMoranAntonelliChenBoland2019} introduced the \textit{nonparametric spike-and-slab lasso} (NPSSL) for sparse generalized additive models (GAMs). Under this model, the response surface is decomposed into the sum of univariate functions,
\begin{align} \label{GAM}
	y_i = \sum_{j=1}^{p} f_j(x_{ij}) + \varepsilon_i, \hspace{.5cm} \varepsilon_i \overset{iid}{\sim} \mathcal{N}(0, \sigma^2).
\end{align}
In \cite{BaiMoranAntonelliChenBoland2019}, each of the $f_j$'s is approximated using a basis expansion, or a linear combination of basis functions $\mathcal{B}_j = \{ g_{j1}, \ldots, g_{jd} \}$, i.e.
\begin{equation} \label{basisexpansion}
	f_j(x_{ij}) \approx \sum_{k=1}^{d} g_{jk} ( x_{ij}) \beta_{jk}.
\end{equation}
Under sparsity, most of the $f_j$'s in \eqref{GAM} are assumed to be $f_j = 0$. This is equivalent to assuming that most of the weight vectors $\boldsymbol{\beta}_j = (\beta_{j1}, \ldots, \beta_{jd})^T$, $j=1, \ldots, p$, in \eqref{basisexpansion} are equal to $\boldsymbol{0}_d$. The NPSSL is implemented by endowing the basis coefficients $\boldsymbol{\beta} = (\boldsymbol{\beta}_1^T, \ldots, \boldsymbol{\beta}_p^T)^T$ with the SSGL prior \eqref{ssgroupLASSO} to simultaneously select and estimate nonzero functionals. In addition, \cite{BaiMoranAntonelliChenBoland2019} also extended the NPSSL to identify and estimate the effects of nonlinear interaction terms $f_{rs} (X_{ir}, X_{is}), r \neq s$.

\textbf{Functional regression.} The spike-and-slab lasso methodology has also been extended to functional regression, where the response $y(t)$ is a function that \textit{varies} over some continuum $T$ (often time) A very popular model in this framework is the nonparametric varying coefficient model,
\begin{equation} \label{vcmodel}
	y_i(t) = \sum_{k=1}^{p} x_{ik}(t) \beta_k(t) + \varepsilon_i(t), \hspace{.5cm} t \in T,
\end{equation}
where $y_i(t)$ and $x_{ik}(t)$ are time-varying responses and covariates respectively and $\varepsilon_i(t)$ is a zero-mean stochastic process which captures the within-subject temporal correlations for the $i$th subject. Under \eqref{vcmodel}, the $\beta_k(t)$'s are smooth functions of time (possibly $\beta_k(t) = 0$ for all $t \in T$), and our primary interest is in estimation and variable selection from the $\beta_k(t)$'s, $k=1, \ldots, p$.

\hspace{.3cm} \cite{BaiBolandChen2019} introduced the nonparametric varying coefficient spike-and-slab lasso (NVC-SSL) to simultaneously select and estimate the smooth functions $\beta_k(t), k=1, \ldots, p$. Similarly as with the NPSSL, these functions are approximated using basis expansions of smoothing splines, and the basis coefficients are endowed with the SSGL prior \eqref{ssgroupLASSO}. Unlike GAMs, however, the NVC-SSL model does \textit{not} assume homoscedastic, independent error terms. Instead, the NVC-SSL model accounts for within-subject temporal correlations in its estimation procedure.

\textbf{False discovery rate control with missing data}. Sorted L-One Penalized Estimator (SLOPE) is an elaboration of the LASSO tailored to false discovery control by assigning more penalty to the larger coefficients \cite{BodganVanDenBergSabattiSuCandes2015}.  SSL, on the other hand, penalizes large coefficients less and its false discovery rate is ultimately determined by a combination of the  prior inclusion weight $\theta$ and penalties $\lambda_1$ and $\lambda_0$. \cite{JiangBogdanJosseMiasojedowRockovaTBG2019} propose a hybrid procedure called adaptive Bayesian SLOPE, which effectively combines the SLOPE method (sorted $l_1$ regularization) together with the SSL method in the context of variable selection with missing covariate values.
As with SSL, the coefficients are regarded as arising from a hierarchical model consisting of two groups: (1) the spike for the inactive and (2) the slab for the active. However, instead of assigning spike priors for each covariate, they propose a joint ``SLOPE'' spike prior which takes into account ordering of coefficient magnitudes in order to control for false discoveries.

\section{Theoretical properties} \label{SSLTheory}

In addition to its computational tractability and its excellent finite-sample performance, the spike-and-slab LASSO has also been shown to provide strong theoretical guarantees.  Although this paper focuses mainly on methodology, we briefly outline a few of the major theoretical developments for spike-and-slab LASSO methods.

A common theme in Bayesian asymptotic theory is the study of the learning rate of posterior point estimates (such as the mean, median or mode) and/or of the full posterior. Working under the frequentist assumption of a ``true'' underlying model, the aim under the former is to study the \textit{estimation} rate of point estimators under a given risk function, such as expected squared error loss. From a fully Bayes perspective, one may also be interested in the \textit{posterior contraction rate}, or the speed at which the \textit{entire} posterior contracts around the truth. In both cases, {  the frequentist minimax estimation rate is a useful benchmark, since the posterior cannot contract faster than this rate} \cite{GhosalGhoshVanDerVaart2000}.

In a variety of contexts, including the Gaussian sequence model, sparse linear regression, and graphical models, the SSL global posterior mode has been shown to achieve the {  minimax estimation rate} \cite{GanNarisettyLiang2018,Rockova2018,RockovaGeorge2016Abel,RockovaGeorge2018}. From a fully Bayesian perspective, the \textit{entire} posterior under SSL or SSL-type priors has \textit{also} been shown to achieve (near) optimal posterior contraction rates in the contexts of the Gaussian sequence model, linear regression, regression with grouped variables, nonparametric additive regression, and functional regression \cite{BaiBolandChen2019, BaiMoranAntonelliChenBoland2019, NieRockova2020, Rockova2018, RockovaGeorge2016Abel,RockovaGeorge2018}. It is not necessarily the case that the posterior mode and the full posterior contract at the same rate \cite{CastilloMismer2018,CastilloSchmidtHieberVanDerVaart2015}. These theoretical results thus show that the SSL is optimal from \textit{both} penalized likelihood \textit{and} fully Bayesian perspectives.

\section{Discussion} \label{Discussion}

In this paper, we have reviewed the spike-and-slab LASSO \eqref{SSLprior}. The SSL forms a continuum between the {  penalized likelihood LASSO and the Bayesian point-mass spike-and-slab frameworks}, borrowing strength from both constructs while addressing limitations of each. First, the SSL employs a \textit{non}-separable penalty that self-adapts to ensemble information about sparsity and that performs selective shrinkage. Second, the SSL is amenable to fast maximum \textit{a posteriori} finding algorithms which can be implemented in a highly efficient, scalable manner. Third, the posterior mode under the SSL prior can automatically perform both variable selection and estimation. Finally, uncertainty quantification for the SSL can be attained by either debiasing the posterior modal estimate or by utilizing efficient approaches to posterior sampling. Beyond linear regression, the spike-and-slab LASSO methodology is broadly applicable to a wide number of statistical problems, including generalized linear models, factor analysis, graphical models, and nonparametric regression.

{
	\section*{Acknowledgments}
	This work was supported by funding from the University of South Carolina College of Arts \& Sciences, the James S. Kemper Foundation Faculty Research Fund at the University of Chicago Booth School of Business, and NSF Grants DMS-1916245, DMS-1944740.
}

\bibliographystyle{apa}