Spike-and-Slab Prior for Covariate Selection

Summary

The spike-and-slab prior enables automatic selection of which control time series to include in the regression component of the BSTS model. The “spike” places point mass at zero (a variable is excluded); the “slab” is a weakly informative Gaussian for included variables. A Gibbs sampler updates the binary inclusion vector jointly with regression coefficients, allowing selection from tens or hundreds of candidates without overfitting.

Overview

When constructing a synthetic control, there may be many candidate predictor series. Including all of them would cause overfitting; pre-selecting a subset requires strong domain knowledge. The spike-and-slab prior provides an automatic Bayesian solution.

Main Content

Definition: Spike-and-Slab Prior

Let $\varrho =({\varrho }_1,\dots ,{\varrho }_J)$ be binary inclusion indicators, where ${\varrho }_j=1$ if ${\beta }_j\ne 0$ and ${\varrho }_j=0$ otherwise. The full prior factorizes as:

$$p(\varrho ,\beta ,1/{\sigma }_{\epsilon }^2)=p(\varrho )\cdot p({\sigma }_{\epsilon }^2\mid \varrho )\cdot p({\beta }_{\varrho }\mid \varrho ,{\sigma }_{\epsilon }^2)\text{\hspace{1em}(2.8)}$$

Spike (inclusion probability):

$$p(\varrho )=\prod _{j=1}^J{\pi }_j^{{\varrho }_j}(1-{\pi }_j{)}^{1-{\varrho }_j}\text{\hspace{1em}(2.9)}$$

where ${\pi }_j$ is the prior probability of including predictor $j$.

Slab (Gaussian prior for included coefficients):

$${\beta }_{\varrho }\mid \varrho ,{\sigma }_{\epsilon }^2\sim \mathcal{N}({\mathbf{b}}_{\varrho },{\sigma }_{\epsilon }^2({\Sigma }_{\varrho }^{-1}{)}^{-1})\text{\hspace{1em}(2.10)}$$

Inverse-Gamma prior on error variance:

$$\frac{1}{{\sigma }_{\epsilon }^2}\sim \mathcal{G}\left(\frac{{\nu }_{\epsilon }}{2},\frac{s_{\epsilon }}{2}\right)\text{\hspace{1em}(2.11)}$$

Setting the Inclusion Prior ${\pi }_j$

Rather than setting ${\pi }_j$ individually, the recommended approach is to elicit an expected model size $M$ and set:

$${\pi }_j=\frac{M}{J}$$

This scales naturally with $J$ (total number of predictors) and avoids having to specify a hierarchical prior.

Special cases:

• ${\pi }_j=1$: Force predictor $j$ into the model
• ${\pi }_j=0$: Force predictor $j$ out of the model

Zellner’s g-Prior for the Slab

The precision matrix ${\Sigma }^{-1}$ in equation (2.10) uses a g-prior (Zellner 1986):

$${\Sigma }^{-1}=\frac{g}{n}\left\{w{X}^{\top }X+(1-w)\text{diag}(X^{\top }X)\right\}\text{\hspace{1em}(2.12)}$$

• $g$: number of observations worth of prior weight; default $g=1$
• $w$: mixing weight between $X^{\top }X$ (full correlation structure) and diagonal (independent); default $w=0.5$
• $n$: number of observations

Interpretation: $g$ observations worth of prior information; $g/n\cdot X^{\top }X$ is the prior information matrix. The averaging with the diagonal ensures propriety when $X^{\top }X$ is not positive definite.

Sufficient Statistics for Posterior Sampling

Given the spike-and-slab structure, the posterior sufficient statistics for $(\varrho ,\beta ,{\sigma }_{\epsilon }^2)$ are:

$$V_{\varrho }^{-1}=(X^{\top }X{)}_{\varrho }+{\Sigma }_{\varrho }^{-1},{\tilde{\beta }}_{\varrho }=(V_{\varrho }^{-1}{)}^{-1}(X_{\varrho }^{\top }{\dot{y}}_{1:n}+{\Sigma }_{\varrho }^{-1}{b}_{\varrho })\text{\hspace{1em}(2.13)}$$ $$N={\nu }_{\epsilon }+n,S_{\varrho }=s_{\epsilon }+{\dot{y}}_{1:n}^{\top }{\dot{y}}_{1:n}+b_{\varrho }^{\top }{\Sigma }_{\varrho }^{-1}{b}_{\varrho }-{\tilde{\beta }}_{\varrho }^{\top }{V}_{\varrho }^{-1}{\tilde{\beta }}_{\varrho }$$

These are updated efficiently in the Gibbs sampler by drawing each ${\varrho }_j$ given ${\varrho }_{-j}$ (all others fixed).

Computational Efficiency

Each full-conditional $p({\varrho }_j\mid {\varrho }_{-j},\dots )$ evaluates easily because ${\varrho }_j$ takes only two values. The dimension of matrices in (2.13) is ${\sum }_j{\varrho }_j$ (number of included variables), which is small if the model is truly sparse. Thus even with hundreds of candidates, the algorithm is fast.

Connections

• Used by MCMC Inference for CausalImpact — Gibbs sampler updates $\varrho$
• Part of Bayesian Structural Time-Series Model — applies to the static regression component
• Related to Nonparametric Models Overview — different approach to regularization
• In the advertising application (CausalImpact Empirical Application), expected model size $M=3$ with $R^2\approx 0.8$ expected

See Also

• Bayesian Structural Time-Series Model — full model
• MCMC Inference for CausalImpact — how this prior is sampled
• The Horseshoe Prior — continuous global-local alternative to spike-and-slab for sparsity