Structured Shrinkage Priors

TL;DR

This paper introduces structured shrinkage priors for regression coefficients that incorporate known structural relationships while maintaining sparsity, using elliptical slice sampling for efficient computation.

Contribution

It develops a flexible class of priors that combine structure and sparsity, addressing computational challenges with elliptical slice sampling.

Findings

Priors effectively encode structure and sparsity in regression models.

Elliptical slice sampling enables practical inference with these priors.

Structured priors improve model interpretability and performance.

Abstract

In many regression settings the unknown coefficients may have some known structure, for instance they may be ordered in space or correspond to a vectorized matrix or tensor. At the same time, the unknown coefficients may be sparse, with many nearly or exactly equal to zero. However, many commonly used priors and corresponding penalties for coefficients do not encourage simultaneously structured and sparse estimates. In this paper we develop structured shrinkage priors that generalize multivariate normal, Laplace, exponential power and normal-gamma priors. These priors allow the regression coefficients to be correlated a priori without sacrificing elementwise sparsity or shrinkage. The primary challenges in working with these structured shrinkage priors are computational, as the corresponding penalties are intractable integrals and the full conditional distributions that are needed to approximate the posterior mode or simulate from the posterior distribution may be non-standard. We overcome these issues using a flexible elliptical slice sampling procedure, and demonstrate that these priors can be used to introduce structure while preserving sparsity.

Tables

Table 1. Table 1 : Properties of specific SHP priors.

Property	SNG	SPB	SPN
Generalizes an Independent Shrinkage Prior	✓	✓	✓
Univariate Marginal is an Independent Shrinkage Prior	✓	✓	✓
Generalizes a Laplace Prior	$c=1$	$q=1$
Generalizes a Normal Prior	$c\to \infty$	$q\to 2$
Infinite Spike or Pole at Zero		$c\le 1/2$	✓
Quadratic Scale Log Full Conditional			✓
Arbitrary Covariance Structure Achievable			✓

Table 2. Table 2 : True positive (TP) and false positive (FP) rates of approximate 90% intervals based on comparison to “ground truth” 90% confidence intervals that do not include zero based on simple logistic regression estimates computed from all 240 trials.

	SNO	SNG				SPN	SPB
		$c=10$	$c=1$	$c=0.5$	$c=0.3$		$q=0.65$	$q=0.78$	$q=1$	$q=1.75$
TP%	18.18	18.18	13.64	0.00	0.00	0.00	18.18	13.64	11.36	9.09
FP%	17.68	17.68	10.98	3.05	0.00	0.61	15.85	10.98	6.71	6.10

Table 3. Table 3 : Minimum effective sample sizes for 50,000 samples collected from 20 chains of 13,500 samples each, discarding initial 1 , 000 1 000 1,000 samples from each chain as burn-in and thinning by a factor of 5 5 5 .

	SNO	SNG
		$c=0.3$	$c=0.5$	$c=1$	$c=10$
$𝜷$ Only	$24,648.14$	$976.34$	$1,474.14$	$1,980.6$	$25,063.38$
All Parameters	$18,975.88$	$583.2$	$797.4$	$594.44$	$19,371.21$
	SPN	SPB
		$q=0.65$	$q=0.78$	$q=1$	$q=1.75$
$𝜷$ Only	$14,075.65$	$1,300.39$	$1,519.00$	$2,342.22$	$26,764.76$
All Parameters	$1,839.08$	$423.62$	$329.33$	$365.43$	$218.41$

Equations

$$\displaystyle\hat{\beta}=\text{argmin}_{\beta}\frac{1}{2\phi^{2}}\left(\hat{\beta}_{OLS}-\beta\right)^{2}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\text{log}\left(\int p\left(\beta|s,\omega\right)p\left(s|\theta\right)ds\right),$$

$$\displaystyle\hat{\boldsymbol{\beta}}=\text{argmin}_{\boldsymbol{\beta}}\frac{1}{2\phi^{2}}\left|\left|\hat{\boldsymbol{\beta}}_{OLS}-\boldsymbol{\beta}\right|\right|^{2}_{2}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\text{log}\left(\int p\left(\boldsymbol{\beta}|\boldsymbol{s},\boldsymbol{\Omega}\right)p\left(\boldsymbol{s}|\boldsymbol{\theta}\right)d\boldsymbol{s}\right).$$

$$p(\mathbf{s}|\mathbf{\beta },\mathbf{\Omega },\mathbf{\theta })\propto$$

$$p({\mathbf{s}}^{\prime }|\mathbf{\beta },\mathbf{\Omega },\mathbf{\theta })\propto$$

$$\displaystyle p\left(\boldsymbol{s}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}|\boldsymbol{\beta},\boldsymbol{\Omega},\boldsymbol{\theta}\right)\propto$$

$$\displaystyle\text{exp}\left\{-\frac{1}{2}\left(\boldsymbol{1}/\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left(\boldsymbol{V}^{1/2}\boldsymbol{s}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+\boldsymbol{m}\right)\right)^{\prime}\left(\boldsymbol{\Omega}^{-1}\circ\left(\boldsymbol{\beta}\boldsymbol{\beta}^{\prime}\right)\right)\left(\boldsymbol{1}/\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left(\boldsymbol{V}^{1/2}\boldsymbol{s}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+\boldsymbol{m}\right)\right)\right)\times$$

$$\displaystyle p\left(\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f\left(\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{V}^{1/2}\boldsymbol{s}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+\boldsymbol{m}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\right)\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}|\boldsymbol{\theta}\right)$$

$$\displaystyle p\left(\boldsymbol{s}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{\prime\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\boldsymbol{u}\text{sin}(\boldsymbol{\pi})+\boldsymbol{r}\text{cos}(\boldsymbol{\pi})|\boldsymbol{\beta},\boldsymbol{\Omega},\boldsymbol{\theta}\right)\prod_{j=1}^{p}\text{exp}\left\{-\frac{1}{2}\left(r_{j}^{2}+u_{j}^{2}\right)\right\}/\text{exp}\left\{-\frac{s_{j}{{}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}}^{2}}{2}\right\}.$$

$$\mathbf{z}|\mathbf{X},\mathbf{y},\mathbf{s},\mathbf{\Omega }$$

$$\mathbf{s}|\mathbf{X},\mathbf{y},\mathbf{z},\mathbf{\Psi }$$

$$\mathbb{E}\left[\mathbf{\beta }\right]=\mathbb{E}\left[\mathbf{s}\right]\circ \mathbb{E}\left[\mathbf{z}\right]\text{ and }\mathbf{\Sigma }=\mathbb{E}\left[\mathbf{s}{\mathbf{s}}^{\prime }\right]\circ \mathbb{E}\left[\mathbf{z}{\mathbf{z}}^{\prime }\right].$$

$$y_i\stackrel{indep.}{\sim }\text{Binomial}((1+\text{exp}\{-{\gamma }_{jk}-{\beta }_{jk}{x}_{ijk}\}{)}^{-1}),$$

$$\frac{\mathbb{E}\left[{\beta }_j^2{\beta }_k^2\right]}{{\sigma }_{jj}{\sigma }_{kk}}$$

$$\frac{\mathbb{E}\left[{\beta }_j^2{\beta }_k{\beta }_l\right]}{{\sigma }_{jj}\sqrt{{\sigma }_{kk}{\sigma }_{ll}}}$$

$$s_j^2\stackrel{d}{=}\frac{\Gamma \left(\frac{1}{2\alpha }\right){\xi }_j^{\frac{1-\alpha }{\alpha }}}{2\Gamma \left(\frac{3}{2\alpha }\right)}\text{, }{\xi }_j\stackrel{i.i.d.}{\sim }\text{gamma}\left(\text{shape}=\frac{1+\alpha }{2\alpha },\text{rate}=f\left({\delta }_j|\alpha \right)\right)\text{ and }$$

$$p\left({\delta }_j|\alpha \right)\propto f{\left({\delta }_j|\alpha \right)}^{\frac{\alpha -1}{2\alpha }},$$

$$p\left(\mathbf{\beta }\right)$$

$$p\left({\beta }_1\right)$$

$$=\int p\left({\beta }_1/s_1\right)p\left(s_1\right)/\left|s_1\right|\underset{(\ast )}{\underbrace{\left(\int p\left({\mathbf{\beta }}_{-1}/{\mathbf{s}}_{-1}|{\beta }_1/s_1\right)\left(\prod _{j=2}^p\frac{p\left(s_j\right)}{\left|s_j\right|}\right)d{s}_2\dots d{s}_{p}d{\beta }_2\dots d{\beta }_p\right)}}d{s}_1.$$

$${\int }_{-\infty }^{\infty }\frac{1}{\left|s\right|}\text{exp}\left\{-\alpha \left(s^2-\gamma s\right)\right\}ds=$$

$${\int }_0^{\gamma }\frac{1}{s}\text{exp}{\left\{-\alpha s\right\}}^{s-\gamma }ds+{\int }_{\gamma }^{\infty }\frac{1}{s}\text{exp}{\left\{-\alpha s\right\}}^{s-\gamma }ds.$$

$${\int }_0^{\gamma }\frac{1}{s}\text{exp}{\left\{-\alpha s\right\}}^{s-\gamma }ds$$

$$=\left(\text{ln}\left(\gamma \right)-{\text{lim}}_{a\to 0+}\text{ln}\left(a\right)\right)=+\infty .$$

$${\int }_{-\infty }^{\infty }\frac{1}{\left|s\right|}\text{exp}\left\{-\alpha \left(s^2-\gamma s\right)\right\}ds=$$

$${\int }_{\gamma }^0-\frac{1}{s}\text{exp}{\left\{-\alpha s\right\}}^{s-\gamma }ds+{\int }_0^{\infty }\frac{1}{s}\text{exp}{\left\{-\alpha s\right\}}^{s-\gamma }ds.$$

$${\int }_{\gamma }^0-\frac{1}{s}\text{exp}{\left\{-\alpha s\right\}}^{s-\gamma }ds$$

$$={\int }_0^{-\gamma }\frac{1}{s}ds=+\infty .$$

$${\int }_{-\infty }^{\infty }\frac{1}{\left|s\right|}\text{exp}\left\{-\alpha s^2\right\}ds=$$

$${\int }_0^{1/\sqrt{\alpha }}\frac{1}{s}\text{exp}\left\{-\alpha s^2\right\}ds$$

$$\ge {\int }_0^{1/\sqrt{\alpha }}\frac{1-\sqrt{\alpha }s}{s}ds$$

$$={\int }_0^{1/\sqrt{\alpha }}\frac{1}{s}ds-{\int }_0^{1/\sqrt{\alpha }}\alpha ds$$

$$=\text{ln}\left(1/\sqrt{\alpha }\right)-{\text{lim}}_{a\to 0+}\text{ln}\left(a\right)-\sqrt{\alpha }=+\infty .$$

$$p\left(\left(0,{\mathbf{b}}_{-1}\right)|\mathbf{\Psi },\mathbf{\Omega }\right)$$

$$=\int$$

$$\underset{(\ast )}{\underbrace{{\int }_{-\infty }^{\infty }\frac{1}{\left|s_1\right|}\text{exp}\left\{-\frac{1}{2}\left(s_1^2{\left({\mathbf{\Psi }}^{-1}\right)}_{11}-2s_1{\left({\mathbf{\Psi }}^{-1}\right)}_{-1,1}^{\prime }{\mathbf{s}}_{-1}\right)\right\}d{s}_1}}d{\mathbf{s}}_{-1}.$$

$$p\left({\beta }_1=0|\mathbf{\Psi },\mathbf{\Omega }\right)$$

$$=\int$$

$$\underset{(\ast )}{\underbrace{{\int }_{-\infty }^{\infty }\frac{1}{\left|s_1\right|}\text{exp}\left\{-\frac{1}{2}\left(s_1^2{\left({\mathbf{\Psi }}^{-1}\right)}_{11}-2s_1{\left({\mathbf{\Psi }}^{-1}\right)}_{-1,1}^{\prime }{\mathbf{s}}_{-1}\right)\right\}d{s}_1}}d{\mathbf{s}}_{-1}d{\mathbf{\beta }}_{-1}.$$

$${\int }_0^{\infty }(s^2{)}^{c-1/2-1}\text{exp}\left\{-c{s}^2\right\}d{s}^2={\int }_0^{1/c}(s^2{)}^{c-1/2-1}\text{exp}\left\{-c{s}^2\right\}d{s}^2+{\int }_{1/c}^{\infty }(s^2{)}^{c-1/2-1}\text{exp}\left\{-c{s}^2\right\}d{s}^2.$$

Full text

Keywords: Bayesian lasso, Lasso, multivariate Laplace, multivariate normal-scale mixture, sparsity, shrinkage.

1 Introduction

Shrinkage prior-based penalized estimates of regression coefficients are ubiquitous and useful. When we observe an $n\times 1$ vector of responses $\boldsymbol{y}$ and an $n\times p$ matrix of regressors $\boldsymbol{X}$ and the data are high dimensional, i.e. $p$ is large relative $n$, traditional regression methods can fail. They may produce estimates of the $p\times 1$ vector of regression coefficients $\boldsymbol{\beta}$ have prohibitively large variance or are not unique because the data provide relatively little information about the unknown regression coefficients.

Using a prior for $\boldsymbol{\beta}$ that reflects our a priori knowledge, we can obtain better estimates of $\boldsymbol{\beta}$. When our a priori knowledge involves similarities and differences among elements of $\boldsymbol{\beta}$, we might assume a structured mean-zero multivariate normal prior with covariance matrix $\boldsymbol{\Sigma}$. Alternatively, when our a priori knowledge involves magnitudes of elements of $\boldsymbol{\beta}$, we might assume a sparsity inducing mean-zero independent Laplace prior with variance $\sigma^{2}$ in order to encode information about the origin and tail behavior of $\boldsymbol{\beta}$. This is useful when $\boldsymbol{\beta}$ is expected to be sparse, as the posterior mode of $\boldsymbol{\beta}$ under this prior corresponds to the $\ell_{1}$ or lasso penalized estimate of $\boldsymbol{\beta}$ (Tibshirani, 1996; Park and Casella, 2008).

When our a priori knowledge involves both similarities and differences among and magnitudes of elements of $\boldsymbol{\beta}$, it would be desirable to assume a structured sparsity inducing shrinkage prior. As an example, consider the analysis of brain-computer interface data. Brain-computer interfaces (BCIs) are used to detect changes in subjects’ cognitive state from contemporaneous electroencephalography (EEG) measurements at different channels corresponding to physical locations on the subject’s skull, which can be collected non-invasively at high temporal resolution (Makeig et al., 2012; Wolpaw and Wolpaw, 2012). We consider the P300 speller, a specific BCI device which is designed to detect when a subject is viewing a specified target letter (Forney et al., 2013). For an individual subject, we consider $20$ indicators $\boldsymbol{y}$ for whether or not the subject was viewing a specified target letter during trial $i$ and $20$ contemporaneous EEG measurements $\boldsymbol{x}_{jk}$ from time point $j$ and channel $k$, for $j=1,\dots,208$ time points and $k=1,\dots,8$ channels. A total of $240$ indicators and $240$ contemporaneous EEG measurements are available, but we consider a subset of $20$ indicators and contemporaneous EEG measurements because performance with such little data is of special interest and especially challenging. Scientifically, we expect to observe a P300 wave in contemporaneous EEG measurements during trials when the target letter is present, which is characterized by a sharp rise and then dip before returning to equilibrium. We expect that the wave will begin shortly after the target letter is shown, and will be observed earlier and more clearly on some channels than others.

Estimated regression coefficients that describe the relationship between EEG measurements from time point $j$ and channel $k$ and whether or not the subject was viewing the target letter can be obtained by regressing the indicators $\boldsymbol{y}$ on an intercept and measurements for each time point and channel $\boldsymbol{x}_{jk}$ using a logistic regression model, which assumes $y_{i}\stackrel{{\scriptstyle indep.}}{{\sim}}\text{Binomial}((1+\text{exp}\{-\gamma_{jk}-\beta_{jk}x_{ijk}\})^{-1})$. Because EEG data are noisy, making the P300 wave difficult to observe, regression coefficients $\beta_{jk}$ may be poorly estimated. It can be desirable to assume a model for the regression coefficients that incorporates the scientific context suggesting that $\boldsymbol{\beta}$ should be sparse and structured, as only $\beta_{jk}$ that correspond to time points where the P300 wave occurs are expected to be nonzero and pairs of coefficients $(\beta_{jk},\beta_{j^{\prime}k^{\prime}})$ corresponding to nearby time points or channels are expected to be similar. We observe evidence of such behavior in Figure 1, which shows estimated logistic regression coefficients $\hat{\beta}_{logit,jk}$. Consecutive channels correspond to neighboring locations; the appendix includes a map of channel locations. Estimated regression coefficients corresponding to consecutive time points or neighboring channels tend to share the same sign.

This suggests a need for structured shrinkage priors that (i) can incorporate a priori knowledge of structure by allowing elements of $\boldsymbol{\beta}$ to covary with nondiagonal prior covariance matrix $\boldsymbol{\Sigma}=\mathbb{V}[\boldsymbol{\beta}]$, (ii) can incorporate a priori knowledge of sparsity by allowing elementwise shrinkage of $\boldsymbol{\beta}$, and (iii) span the range of common structured and sparsity inducing priors by generalizing multivariate normal and independent Laplace priors. However, existing approaches fail to satisfy all three of these criteria. We can see this by considering the normal scale-mixture representations of many common priors, which are used in Bayesian approaches to sparse regression (Polson and Scott, 2010). Letting ‘$\circ$’ be the Hadamard elementwise product, $\boldsymbol{z}\sim\text{normal}(\boldsymbol{0},\boldsymbol{\Omega})$ and $\boldsymbol{s}$ be a vector of stochastic scales that are independent of $\boldsymbol{z}$, a prior distribution for $\boldsymbol{\beta}$ has a normal scale-mixture representation if there exists a density $p(\boldsymbol{s}|\boldsymbol{\theta})$ such that $\boldsymbol{\beta}\stackrel{{\scriptstyle d}}{{=}}\boldsymbol{s}\circ\boldsymbol{z}$. These priors are interpretable from a data generating perspective and have easy-to-compute moments. Specifically the prior covariance matrix of $\boldsymbol{\beta}$ that encodes the a priori knowledge of structure is $\boldsymbol{\Sigma}=\mathbb{E}[\boldsymbol{s}\boldsymbol{s}^{\prime}]\circ\boldsymbol{\Omega}$.

Much literature focuses on the unstructured case where $\boldsymbol{s}$ is a vector of independent, identically distributed elements and $\boldsymbol{\Omega}\propto\boldsymbol{I}_{p}$. This includes the Laplace prior, bridge/exponential power priors, normal-gamma priors, Dirichlet-Laplace priors, and horseshoe priors (Park and Casella, 2008; Griffin and Brown, 2010; Carvalho et al., 2010; Bhattacharya et al., 2015). These priors can incorporate a priori knowledge of sparsity by modeling a separate stochastic scale for every element of $\boldsymbol{\beta}$ and choosing distributions for $\boldsymbol{s}$ that yield possibly sparse posterior modes. However, they do not allow elements of $\boldsymbol{\beta}$ to covary. Because the marginal covariance matrix $\boldsymbol{\Sigma}$ is an elementwise product of $\mathbb{E}[\boldsymbol{s}\boldsymbol{s}^{\prime}]$ and $\boldsymbol{\Omega}$, elements of $\boldsymbol{\beta}$ are uncorrelated a priori under these priors. The same limitation afflicts structured shrinkage priors which model elements of $\boldsymbol{s}$ as correlated but continue to assume that $\boldsymbol{\Omega}\propto\boldsymbol{I}_{p}$ (van Gerven et al., 2010; Kalli and Griffin, 2014; Wu et al., 2014; Zhao et al., 2016; Kowal et al., 2017).

A different strand of literature includes all elliptically contoured prior distributions. It sets $\boldsymbol{s}=s\boldsymbol{1}_{p}$ and models elements of $\boldsymbol{z}$ as potentially correlated with covariance matrix $\boldsymbol{\Omega}$. When $s^{2}$ is assumed to be exponentially distributed and $\boldsymbol{\Omega}\propto\boldsymbol{I}_{p}$, the prior corresponding to the group lasso introduced in Yuan and Lin (2006) is obtained. The more general multivariate Laplace prior introduced by van Gerven et al. (2009) is obtained by assuming exponentially distributed $s^{2}$ and allowing arbitrary $\boldsymbol{\Omega}$. These priors can incorporate a priori knowledge of structure in two ways: by treating elements of $\boldsymbol{\beta}$ as grouped through their shared stochastic scale $s$ and, when $\boldsymbol{\Omega}$ is allowed to take on arbitrary values, by encoding a priori knowledge of structure through specification of $\boldsymbol{\Omega}$ that does not satisfy $\boldsymbol{\Omega}\propto\boldsymbol{I}_{p}$, as $\boldsymbol{\Sigma}\propto\boldsymbol{\Omega}$. However, as these priors only have a single shrinkage factor $s$, they shrink all elements of $\boldsymbol{\beta}$ jointly and posterior modes based on these priors will only produce sparse estimates of $\boldsymbol{\beta}$ for which $\boldsymbol{\beta}=\boldsymbol{0}$ (Simon et al., 2013). Additionally, these priors do not generalize their independent counterparts. For instance, the multivariate Laplace prior with $\boldsymbol{\Omega}=\omega^{2}\boldsymbol{I}_{p}$ does not correspond to the independent Laplace prior, and the corresponding penalty is the group lasso as opposed to the lasso penalty. This is not appropriate when only some elements of $\boldsymbol{\beta}$ are expected to be sparse.

At least one set of priors can incorporate a priori knowledge of structure and sparsity by also modeling elements of the inverse covariance matrix of $\boldsymbol{\Omega}$. This includes the prior distribution that yields the fused lasso penalty $\lambda_{1}||\boldsymbol{\beta}||_{1}+\lambda_{2}\sum_{i-j=1}|\beta_{j}-\beta_{i}|$ and priors that correspond to more general structured penalties of the form $\boldsymbol{\beta}^{\prime}\boldsymbol{Q}^{-1}\boldsymbol{\beta}+\lambda||\boldsymbol{\beta}||_{1}$ where $\boldsymbol{Q}^{-1}$ is positive semidefinite (Kyung et al., 2010; Ng and Abugharbieh, 2011; de Brecht and Yamagishi, 2012). The penalties corresponding to these priors are very popular, as they yield computationally feasible posterior mode optimization problems. Their main limitation is that relating the prior parameters $\lambda_{1}$ and $\lambda_{2}$ or $\boldsymbol{Q}^{-1}$ and $\lambda$ to the prior moments of $\boldsymbol{\beta}$ is prohibitively challenging. This makes understanding exactly how flexible these priors are and specifying values or priors for these parameters difficult, as it is unclear how to relate the kind of a priori knowledge we might have to the prior parameter values.

One last class of relevant priors are the structured normal-gamma priors introduced in Griffin and Brown (2012b) and Griffin and Brown (2012a), which are obtained by assuming $\boldsymbol{\beta}\stackrel{{\scriptstyle d}}{{=}}\boldsymbol{C}\left(\boldsymbol{s}\circ\boldsymbol{z}\right)$, where $\boldsymbol{C}$ is a $p\times q$ rectangular matrix with $p<q$ and $s^{2}_{j}$ and $z_{j}$ are independent gamma and normal random variables, respectively. Both elementwise and structured shrinkage can be simultaneously encouraged by setting $\boldsymbol{C}=\left[\boldsymbol{I}_{p},\boldsymbol{D}\right]$, where $\boldsymbol{D}$ is a $p\times\left(q-p\right)$ matrix with columns that encode groups of elements of $\boldsymbol{\beta}$ that should be penalized jointly. However, the theory that justifies the use of these structured priors is specific to generalizing independent normal-gamma priors.

We construct a class of “structured Hadamard product” (SHP) priors that can incorporate a priori knowledge of structure and sparsity by allowing non-diagonal $\boldsymbol{\Omega}$ and, accordingly, $\boldsymbol{\Sigma}$. The main challenge in using such priors is computational; their development has been limited as a result. The marginal priors for $\boldsymbol{\beta}$ are intractable integrals and correspond to penalties $-\text{log}(\int p(\boldsymbol{\beta}|\boldsymbol{\Omega},\boldsymbol{s})p\left(\boldsymbol{s}|\boldsymbol{\theta}\right)d\boldsymbol{s})$ with no simple closed form. Their use requires computationally demanding Markov Chain Monte Carlo (MCMC) algorithms. Two exceptions are Finegold and Drton (2011), which develops a multivariate $t$-distribution by assuming independent inverse-gamma squared scales $s^{2}_{j}$, and Roy et al. (2021), which develops the structured product normal (SPN) prior by assuming normal scales \boldsymbol{s}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sim\text{normal}\left(\boldsymbol{0},\boldsymbol{\Psi}\right).

This paper proceeds as follows. In Section 2, we describe the SPN prior and introduce the novel structured normal-gamma (SNG) and structured power/bridge (SPB) priors which generalize the independent normal-gamma priors in Caron and Doucet (2008) and Griffin and Brown (2010) and power/bridge priors in Frank and Friedman (1993) and Polson et al. (2014), respectively.

We do not consider a structured generalization of the horseshoe prior because the corresponding prior covariance matrix $\boldsymbol{\Sigma}=\mathbb{E}[\boldsymbol{s}\boldsymbol{s}^{\prime}]\circ\boldsymbol{\Omega}$ is not finite, which makes it difficult to understand how a priori knowledge of structure is incorporated.

We discuss properties of these priors in Section 3. Several properties and how they differ across the SPN, SNG, and SPB priors are previewed in Table 1. In Section 4, we describe how the elliptical slice sampling method of Murray et al. (2010) can be used to overcome computational issues, regardless of the distribution assumed for elements of $\boldsymbol{s}$, and discuss estimation of hyperparameters. We focus on problems where the log-likelihood of the data given $\boldsymbol{\beta}$ and any additional nuisance parameters $\boldsymbol{\phi}$, denoted by $-h(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\beta},\boldsymbol{\phi})$, can be written as conditionally quadratic in $\boldsymbol{\beta}$, i.e. $\text{exp}\{-h(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\beta},\boldsymbol{\phi})\}\propto_{\boldsymbol{\beta}}\text{exp}\{-\frac{1}{2}(\boldsymbol{\beta}^{\prime}\boldsymbol{A}\boldsymbol{\beta}-2\boldsymbol{\beta}^{\prime}\boldsymbol{c})\}$ for some positive definite matrix $\boldsymbol{A}$ and real valued vector $\boldsymbol{c}$. This includes linear regression models and certain latent variable representations of logistic and negative binomial regression models (Polson et al., 2013). In Section 5, we use SHP priors to analyze the data depicted in Figure 1. A discussion follows in Section 6.

2 Structured Shrinkage Priors

We define several “structured Hadamard product” (SHP) prior distributions for $\boldsymbol{\beta}$ of the form $\boldsymbol{\beta}=\boldsymbol{s}\circ\boldsymbol{z}$, where ‘$\circ$’ is the elementwise Hadamard product and $\boldsymbol{z}\sim\text{normal}(\boldsymbol{0},\boldsymbol{\Omega})$ and $\boldsymbol{s}$ is a vector of stochastic scales that are independent of $\boldsymbol{z}$. Because the scales of elements of $\boldsymbol{s}$ are not separately identifiable from diagonal elements of $\boldsymbol{\Omega}$, we parametrize $\boldsymbol{s}$ such that $\mathbb{E}[s^{2}_{j}]=1$ for $j=1,\dots,p$. These priors are mean $\boldsymbol{0}$ and have prior variance $\boldsymbol{\Sigma}=\mathbb{E}[\boldsymbol{s}\boldsymbol{s}^{\prime}]\circ\boldsymbol{\Omega}$.

Structured Product Normal (SPN) Prior:

When $\boldsymbol{s}\sim\text{normal}(\boldsymbol{0},\boldsymbol{\Psi})$ and $\boldsymbol{\Psi}$ is a positive definite matrix with diagonal elements equal to $1$, the SPN prior is obtained. The parameters of the prior distribution for $\boldsymbol{s}$, denoted by $\boldsymbol{\theta}$, correspond to the off-diagonal elements of $\boldsymbol{\Psi}$. Elements of the hyperparameters $\boldsymbol{\Omega}$ and $\boldsymbol{\Psi}$ can be related to prior moments of $\boldsymbol{\beta}$. Diagonal elements of $\boldsymbol{\Omega}$ correspond to prior variances of $\boldsymbol{\beta}$, and off diagonal elements of $\boldsymbol{\Omega}$ and $\boldsymbol{\Psi}$ determine covariances and fourth-order prior cross moments of $\boldsymbol{\beta}$, as shown in the appendix. Originally discussed as a sparsity inducing penalty in Hoff (2017) and later implemented as a prior distribution in Roy et al. (2021), the SPN prior is uniquely computationally simple to work with as the full conditional distributions of $\boldsymbol{z}$ and $\boldsymbol{s}$ are both multivariate normal distributions when the log-likelihood $-h(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\beta},\boldsymbol{\phi})$ is quadratic or conditionally quadratic in $\boldsymbol{\beta}$. This is described in greater detail in Section 4. The SPN prior is also appealing as it is the only prior we define that has correlations among elements of $\boldsymbol{s}$. To reduce the number of freely varying parameters and to facilitate use of the SPN prior in settings where it is challenging to estimate or formulate prior opinions about fourth-order prior cross moments of $\boldsymbol{\beta}$, we also define a special case of the SPN prior which we call the symmetric SPN (sSPN) prior. The sSPN prior requires that all elements of $\boldsymbol{\Psi}$ be positive and that the correlation matrix corresponding to $\boldsymbol{\Omega}$ has elements that are equal in magnitude to the elements of $\boldsymbol{\Psi}$, i.e. $|\omega_{ij}/\sqrt{\omega_{ii}\omega_{jj}}|=\psi_{ij}$. Under this constraint, $\boldsymbol{\Psi}$ is a deterministic function of $\boldsymbol{\Omega}$ and $\boldsymbol{\theta}$ is an empty vector.

Structured Normal-Gamma (SNG) Prior:

When the squared stochastic scales $s^{2}_{j}$ are independent gamma random variables $s^{2}_{j}\sim\text{gamma}(c,c)$ with $\mathbb{E}[s^{2}_{j}]=1$ for fixed $c\in(0,\infty)$ and the stochastic scales $s_{j}=\sqrt{s^{2}_{j}}$ are strictly positive, the structured normal-gamma (SNG) prior is obtained. The SNG prior generalizes the normal-gamma prior of Griffin and Brown (2010), which is obtained by setting $\boldsymbol{\Omega}\propto\boldsymbol{I}_{p}$. The fixed shape parameter $c$ parametrizes the prior class and $\boldsymbol{\theta}$ is an empty vector. The value chosen for $c$ determines the prior fourth order moments of $\boldsymbol{\beta}$; SNG priors with smaller values of $c$ have lighter tails.

Structured Power/Bridge (SPB) Prior:

When the squared stochastic scales $s^{2}_{j}$ are independently distributed according to a polynomially tilted positive $\alpha$-stable distribution with index of stability $\alpha=q/2$ and $\mathbb{E}[s^{2}_{j}]=1$ for fixed $q\in(0,2)$ and the stochastic scales $s_{j}=\sqrt{s^{2}_{j}}$ are strictly positive, the structured power/bridge (SPB) prior is obtained. The SPB prior generalizes the bridge or exponential power prior discussed in Polson et al. (2014). The fixed shape parameter $q$ parametrizes the prior class and $\boldsymbol{\theta}$ is an empty vector. Working with this prior is especially computationally challenging because the polynomially tilted positive $\alpha$-stable density $p(s_{j}|\boldsymbol{\theta})$ is not available in closed form. Fortunately, a polynomially tilted positive $\alpha$-stable random variable can be represented as a rate mixture of generalized gamma random variables (Devroye, 2009). A more detailed description of this representation and how it enables computation under the SPB prior is provided in the appendix. As with the SNG prior, the value chosen for $q$ for the SPB prior determines the prior fourth order moments of $\boldsymbol{\beta}$; SPB priors with larger values of $q$ have lighter tails.

Relationships Between Priors:

Figure 2 shows relationships between the SPN, SNG, and SPB priors, the independent Laplace prior, and the multivariate normal prior. When $c=q=1$, the SNG and SPB priors are equivalent and generalize the independent Laplace prior. In the limit as $c\rightarrow\infty$ or $q\rightarrow 2$, the SNG and SPB priors generalize the multivariate normal prior. The SPN prior does not generalize the Laplace prior or the multivariate normal prior, but is equivalent to the SNG prior with $c=1/2$ when $\boldsymbol{\Psi}$ and $\boldsymbol{\Omega}$ are diagonal.

3 Properties

3.1 Univariate Prior Properties

Introducing structure does not alter the marginal prior distributions of elements $\beta_{j}$. For instance, if we assume a structured prior that generalizes the independent Laplace prior for $\boldsymbol{\beta}$, e.g. an SNG prior with $c=1$ or an SPB prior with $q=1$, the marginal prior distribution of any element $\beta_{j}$ is Laplace. This follows from the stochastic representation of elements $\boldsymbol{\beta}$ under these priors, $\boldsymbol{\beta}=\boldsymbol{s}\circ\boldsymbol{z}$. Recall that $\boldsymbol{s}$ is a vector of stochastic scales, $\boldsymbol{z}\sim\text{normal}(\boldsymbol{0},\boldsymbol{\Omega})$, $\boldsymbol{s}$ and $\boldsymbol{z}$ are independent of each other, and that all three SHP priors are obtained by assuming different distributions for $\boldsymbol{s}$. If we consider an individual element $\beta_{j}=s_{j}z_{j}$ under the SNG or SPB priors, introducing structure does not affect $s_{j}$ at all and does not affect the marginal distribution of $z_{j}$, as the marginal distribution of an element of a correlated normal vector is still normal. Similarly, if we consider an individual element $\beta_{j}=s_{j}z_{j}$ under the SPN prior, introducing structure does not affect the marginal distributions of $s_{j}$ and $z_{j}$, again because the marginal distribution of an element of a correlated normal vector is still normal.

Because introducing structure does not alter the marginal prior distributions of elements $\beta_{j}$, the marginal prior distributions of elements $\beta_{j}$ retain the same sparsity inducing properties of the corresponding independent priors. For instance, independent normal-gamma priors with $c\leq 1/2$ are known to have an infinite spike or pole at $b_{j}=0$, which has been viewed in the literature as a sufficient condition for the recovery of sparse signals (Carvalho et al., 2010; Griffin and Brown, 2010; Bhattacharya et al., 2015). Because introducing structure does not alter the marginal prior distributions of elements $\beta_{j}$, the marginal prior distributions of elements $\beta_{j}$ under SNG priors with $c\leq 1/2$ and under the SPN prior will also have an infinite spike or pole at $b_{j}=0$. Proofs are provided in the appendix.

3.2 Joint Prior Properties

3.2.1 Range of Achievable $\boldsymbol{\Sigma}$

Introducing structure while retaining elementwise shrinkage can come at a cost. Specifically, under the SNG, SPB and SPN priors, preserving elementwise shrinkage limits how correlated elements of $\boldsymbol{\beta}$ can be. Recall that under all three SHP priors, the prior covariance matrix of $\boldsymbol{\beta}$ is $\boldsymbol{\Sigma}=\mathbb{E}[\boldsymbol{s}\boldsymbol{s}^{\prime}]\circ\boldsymbol{\Omega}$. Under the SNG and SPB priors, $\mathbb{E}[\boldsymbol{s}\boldsymbol{s}^{\prime}]$ is constant for fixed values of $c$ or $q$. Diagonal elements of $\mathbb{E}[\boldsymbol{s}\boldsymbol{s}^{\prime}]$ are equal to $1$, whereas off-diagonal elements are less than one in absolute value. Thus, $\mathbb{E}[\boldsymbol{s}\boldsymbol{s}^{\prime}]$ shrinks off-diagonal elements of $\boldsymbol{\Omega}$, reducing dependence. When $\boldsymbol{\beta}\in\mathbb{R}^{2}$, we can explicitly calculate the maximum marginal prior correlation $\rho$ under the SNG and SPB priors as a function of $c$ or $q$. Under the SNG prior, the maximum correlation is equal to $c^{-1}(\Gamma(c+1/2)/\Gamma(c))^{2}$ and under the SPB prior, the maximum correlation is equal to $(\pi/2)(\Gamma(2/q)/\sqrt{\Gamma(1/q)\Gamma(3/q)})^{2}$. When $q=c=1$ and both priors are equivalent, the maximum correlation is equal to $\pi/4\approx 0.79$.

We plot the maximum correlation as a function of kurtosis, a measure of tail behavior, under both priors in Figure 3. We observe greater reductions in the maximum correlation when the kurtosis is higher, and under the SNG prior relative to a SPB prior with equal kurtosis. The restricted range of $\boldsymbol{\Sigma}$ under the SNG and SPB priors is similar to the restricted range of the variance-covariance matrix of the alternative multivariate $t$-distribution introduced in Finegold and Drton (2011). Intuitively, the restricted range of $\boldsymbol{\Sigma}$ relates to the conflict that arises between the properties of the marginal joint density $p(\boldsymbol{\beta}|\boldsymbol{\Omega},\boldsymbol{\theta})$ needed to preserve elementwise shrinkage, specifically concentration of the density along the axes, and the properties of the marginal joint density needed to encourage structure, specifically concentration of the density along a line when $p=2$.

In contrast, the unrestricted SPN prior can accommodate an arbitrary prior covariance $\boldsymbol{\Sigma}$. Given a positive semidefinite prior covariance $\boldsymbol{\Sigma}$, it is always possible to find at least one pair of positive semidefinite covariance matrices $\boldsymbol{\Omega}$ and $\boldsymbol{\Psi}$ that satisfy $\boldsymbol{\Sigma}=\boldsymbol{\Omega}\circ\boldsymbol{\Psi}$ (Styan, 1973). The sSPN prior is less flexible. It is easy to simulate a positive semidefinite covariance matrix $\boldsymbol{\Sigma}$ for which the corresponding values of $\boldsymbol{\Omega}$ or $\boldsymbol{\Psi}$ satisfying $|\omega_{ij}|=|\psi_{ij}|$ are not positive semidefinite, but challenging to explicitly characterize the class of covariance matrices $\boldsymbol{\Sigma}$ that correspond to non-positive semidefinite values of $\boldsymbol{\Omega}$ or $\boldsymbol{\Psi}$.

3.2.2 Copulas

Even when all three SHP priors share the same prior covariance matrix $\boldsymbol{\Sigma}$, the induced dependence structures vary widely. We compare the induced dependence structures by considering $\boldsymbol{\beta}\in\mathbb{R}^{2}$ with unit marginal variances and correlation $\rho=0.5$, and examining corresponding copula densities. Let $F^{SNG,1}_{j}(\beta_{j})$ refer to the marginal prior CDF of $\beta_{j}$ corresponding to one of the SHP prior distributions. We can always write $\beta_{j}\stackrel{{\scriptstyle d}}{{=}}F^{-1}_{j}(u_{j})$, where $u_{j}$ is a random variable with uniform margins. The joint distribution of the $p\times 1$ vector $\boldsymbol{u}$ is called the copula of $\boldsymbol{\beta}$, and it characterizes the induced dependence structure. Even if the inverse CDF’s $F^{-1}_{j}(u_{j})$ are not known, the copula density can be approximated by simulating values of $\boldsymbol{\beta}$, transforming simulated values of $\beta_{1}$ and $\beta_{2}$ into percentiles $u_{1}$ and $u_{2}$, and computing a kernel bivariate density estimate of the percentiles. Figure 4 shows numerical approximations to copula densities for several SHP priors with unit marginal prior variances and marginal prior correlation $\rho=0.5$.

The first four panels in the top and bottom rows show copula densities under SNG and SPB priors with increasingly heavy tails, and the final panels on the top and bottom show copula densities under SPN priors with \color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{12}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\rho^{0.5} and \color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{12}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\rho^{0.1}. The parameters of the SPB priors have been chosen to ensure that the kurtosis, a measure of tail behavior, of each SPB prior is equal to the kurtosis of the SNG prior above it,  e.g. the SNG prior with $c=0.3$ has the same kurtosis as the SPB prior with $q=0.65$.

The dependence structures induced by the SNG and SPB priors are similar. As $c$ or $q\rightarrow 0$ and the SNG and SPB priors become less normal and heavier tailed, the copulas display increasingly strong orthant dependence. This means that as $c$ or $q\rightarrow 0$, the priors concentrate more strongly around values of $\boldsymbol{\beta}$ that have the same sign. The SNG prior displays especially strong orthant dependence and appears to converge to a uniform distribution over the positive and negative orthants as $c\rightarrow 0$. Additionally, as $c$ or $q\rightarrow 0$ the copula densities become more concentrated around the axes where at least one element of $\boldsymbol{\beta}$ is nearly or exactly equal to zero, suggesting that these priors are still encouraging elementwise shrinkage.

Like the SNG and SPB priors, the SPN priors concentrate in the upper-right and lower-left corners, where $\beta_{1}$ and $\beta_{2}$ have the same sign and are large in magnitude. However in comparison to the SNG and SPB priors, both SPN priors also concentrate strongly at the origin, where $\beta_{1}=\beta_{2}=0$, and less strongly along the axes, where $\beta_{1}=0$ or $\beta_{2}=0$. Both SPN priors also concentrate more about values of $\boldsymbol{\beta}$ that are similar in magnitude and opposite in sign. The extent of this depends on how \color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{12}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0} is chosen. Compared to the sSPN prior which sets \color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{12}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0} and \color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\psi_{12}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0} symmetrically such that \color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{12}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\psi_{12}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}, the SPN prior with \color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{12}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\rho^{0.1} and \color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\psi_{12}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\rho^{0.9} concentrates more strongly about values of $\boldsymbol{\beta}$ that are equal in magnitude but opposite in sign. Intuitively, this is due to the fact that $\boldsymbol{\beta}$ is made up of one strongly correlated component and one very weakly correlated component, both of which can take any value on $\mathbb{R}$ when \color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{12}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\rho^{0.1} and \color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\psi_{12}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\rho^{0.9}.

3.2.3 Conditional Prior Distributions

We can explore how introducing structure can inform shrinkage of elements of $\boldsymbol{\beta}$ by examining the induced conditional prior distributions. Figure 5 shows conditional prior distributions $p(\beta_{1}|\beta_{2},\boldsymbol{\Omega},\boldsymbol{\theta})$, for $\boldsymbol{\beta}\in\mathbb{R}^{2}$ with unit marginal prior variance, marginal prior correlations $\rho=\{0,0.5\}$, and $\beta_{2}=\{0,2\}$. Again, the first four panels in the top and bottom rows show conditional prior distributions under SNG and SPB priors with increasingly heavy tails, and the final panels on the top and bottom show conditional prior distributions under SPN priors with \color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{12}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\rho^{0.5} and \color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{12}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\rho^{0.1}. The parameters of the SPB priors have been chosen to ensure that the kurtosis, a measure of tail behavior, of each SPB prior is equal to the kurtosis of the SNG prior above it,  e.g. the SNG prior with $c=0.3$ has the same kurtosis as the SPB prior with $q=0.65$.

Introducing structure via a positive prior correlation between $\beta_{1}$ and $\beta_{2}$ can allow us to use knowledge of $\beta_{2}$ to better estimate $\beta_{1}$. We see that sparsity inducing SHP priors concentrate more strongly about zero than their independent counterparts when $\beta_{2}=0$, and shift their mass towards positive nonzero values of $\beta_{1}$ when $\beta_{2}=2$. The conditional priors reflect differing amounts of origin versus tail dependence. The nearly normal SNG and SPB priors on the far left display relatively more tail dependence, as the structured conditional priors for $\beta_{1}$ differ more from their independent counterparts when $\beta_{2}=2$ than when $\beta_{2}=0$. At the other end of the spectrum, the heavier-than-Laplace tailed SNG and SPB priors and the sSPN prior display more origin dependence than tail dependence, as the structured conditional priors for $\beta_{1}$ concentrate much more strongly about $\beta_{1}=0$ when $\beta_{2}=0$ than their independent counterparts, and still have a mode near $\beta_{1}=0$ even when $\beta_{2}=2$. This is especially striking under the SNG priors with $c\leq 1/2$ and the sSPN prior, which retain a very sharp peak at $\beta_{1}=0$ even when $\beta_{2}=2$. The asymmetric SPN prior with \color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{12}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\rho^{0.1} is in between; the full conditional distribution of $\beta_{1}$ concentrates much more strongly about $\beta_{1}=0$ when $\beta_{2}=0$ and shifts markedly to towards $\beta_{1}=2$ when $\beta_{2}=2$.

3.2.4 Joint Marginal Prior Contours

The joint marginal prior contours, which can be interpreted as contours of a penalty function, provide a powerful tool for studying the properties of different priors. Figure 6 shows Monte Carlo approximations of the contours of the joint marginal priors $\int p(\boldsymbol{\beta}|\boldsymbol{s},\boldsymbol{\Omega})p(\boldsymbol{s}|\boldsymbol{\theta})d\boldsymbol{s}$, for $\boldsymbol{\beta}\in\mathbb{R}^{2}$ with marginal prior variance-covariance matrix $\boldsymbol{\Sigma}=(1-\rho)\boldsymbol{I}_{2}+\rho\boldsymbol{1}_{2}\boldsymbol{1}_{2}^{\prime}$ and marginal prior correlation $\rho=0.5$.

Again, the first four panels in the top and bottom rows show contours of SNG and SPB priors with increasingly heavy tails. The parameters of the SPB priors have been chosen to ensure that the kurtosis, a measure of tail behavior, of each SPB prior is equal to the kurtosis of the SNG prior above it,  e.g. the SNG prior with $c=0.3$ has the same kurtosis as the SPB prior with $q=0.65$. The rightmost panels show contours of SPN priors. The top panel is a symmetric sSPN prior with \color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{12}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\rho^{0.5} and \color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\psi_{12}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\rho^{0.5}, whereas the bottom panel is an asymmetric SPN prior with \color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{12}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\rho^{0.1} and \color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\psi_{12}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\rho^{0.9}.

All three SHP priors encourage similar estimates of $\beta_{1}$ and $\beta_{2}$ by pushing contours away from the origin when $\beta_{1}$ and $\beta_{2}$ have the same sign and pushing the contours towards the origin when $\beta_{1}$ and $\beta_{2}$ have opposite signs. The SPN, SNG with $c\leq 1$ and SPB with $q\leq 1$ priors encourage sparse estimates of $\boldsymbol{\beta}$ by retaining discontinuities of the log marginal prior on the axes. The contours do not necessarily keep the same shape as the value of prior density changes. Contours closer to the origin are more similar to their independent counterparts, with relatively more encouragement of sparsity than structure, whereas contours farther from the origin tend to encourage relatively more structure and less sparsity. This is especially evident under the SPN prior with $\rho_{\omega}=\rho^{0.1}$. It is clear from the contours that these priors are not log-concave when $\rho\neq 0$ and correspond to non-convex penalties for SNG priors with $c\leq 1$, SPB priors with $q\leq 1$, and SPN and sSPN priors. Under the SNG priors with $c\leq 1/2$ and the SPN priors, the contours do not cross the axes. The joint marginal distribution of $\boldsymbol{\beta}$ under the SNG prior with $c\leq 1/2$ has has an infinite spike or pole along the axes, where at least one element is equal to [math].

Proposition 3.1

For the SNG prior with c\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\leq\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1/2 and any $b_{j}=0$, then $p(\boldsymbol{b}|c,\boldsymbol{\Omega})=+\infty$.

This is a multivariate analogue of the univariate marginal prior’s spike or infinite pole at $b=0$ under the SNG prior with $c\leq 1/2$. The joint marginal SPN prior behaves similarly.

Proposition 3.2

For the SPN prior, if any $b_{j}=0$, then $p(\boldsymbol{b}|\boldsymbol{\Psi},\boldsymbol{\Omega})=+\infty$.

Proofs are given in the appendix. This suggests that the SPN and SNG priors with c\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\leq\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1/2 may recover sparse signals well (Carvalho et al., 2010; Griffin and Brown, 2010; Bhattacharya et al., 2015). However the presence of an infinite spike or pole at $b_{j}=0$ also makes the already challenging problem of computing the posterior mode intractable, as the log marginal prior is not only nonconvex but also has infinitely many modes.

3.3 Posterior Properties

One last perspective on the properties of all three SHP priors can be gained by examining how the posterior mode of $\boldsymbol{\beta}$ relates to the unpenalized OLS estimate $\hat{\beta}_{OLS}$. This can offer an intuitive understanding of the properties of estimators obtained under the three different priors. In the sparsity inducing penalty literature this is often characterized by the thresholding function, which is defined in the context of linear model for $\boldsymbol{y}$ with a single covariate $\boldsymbol{x}$ that satisfies $\boldsymbol{x}^{\prime}\boldsymbol{x}=1$ as

[TABLE]

where $\phi^{2}$ refers to assumed variance of deviations of $\hat{\beta}_{OLS}-\beta$. For example, (1) gives the soft-thresholding function when $\beta$ has a mean-zero Laplace distribution.

Because the advantage of working with any of the three three SHP priors is the ability to introduce a priori dependence across elements of vector valued $\boldsymbol{\beta}$, examining the corresponding univariate thresholding functions for scalar $\beta$ given by (1) does not fully explore how the posterior mode of $\boldsymbol{\beta}$ relates to the unpenalized OLS estimate $\hat{\boldsymbol{\beta}}$. Accordingly, we define and examine a bivariate thresholding function. We continue to assume a linear model for $\boldsymbol{y}$, but assume an orthogonal design matrix comprised of two covariates $\boldsymbol{X}$ that satisfies $\boldsymbol{X}^{\prime}\boldsymbol{X}=\text{diag}\left\{\boldsymbol{1}_{2}\right\}$ instead of a single covariate $\boldsymbol{x}$. The bivariate thresholding function relates the posterior mode of both $\beta_{1}$ and $\beta_{2}$ to the noise variance $\phi^{2}$, prior parameters $\boldsymbol{\Omega}$ and $\boldsymbol{\theta}$, and OLS estimates $\hat{\beta}_{OLS,1}$ and $\hat{\beta}_{OLS,2}$, and is given by

[TABLE]

The bivariate thresholding function can be approximated using a Gibbs-within-EM algorithm described in the Section 4. Figure 7 shows approximate bivariate thresholding functions for $\boldsymbol{\beta}\in\mathbb{R}^{2}$ with unit marginal prior variances, marginal prior correlation $\rho=0.5$, noise variance $\phi^{2}=0.1$, $\hat{\beta}_{OLS,2}\in\{-0.5,0,0.5,1\}$ and $\hat{\beta}_{OLS,1}\in\left[0,1\right]$, computed from $1,001,000$ Gibbs sampler iterations, with the first $1,000$ samples discarded as burn-in.

Examination of the bivariate thresholding functions confirms that the SPN prior, SNG prior with $c\leq 1$ and SPB prior with $c\leq 1$ can yield possibly sparse posterior mode estimates of $\boldsymbol{\beta}$, and that the introduction of structure encourages or discourages sparse estimates $\hat{\beta}_{1}$ depending not only on the observed value of $\hat{\beta}_{OLS,1}$ but also the observed value of $\hat{\beta}_{OLS,2}$. There are a few interesting trends that relate to previously identified properties of the SPN, SNG and SPB priors. Both SPN priors shrink $\hat{\beta}_{1}$ towards [math] less aggressively when $\hat{\beta}_{OLS,2}=-0.5$ than when $\hat{\beta}_{OLS,2}=0$, which reflects the tendency of the SPN priors to encourage estimates of $\boldsymbol{\beta}$ with similar absolute magnitudes. Additionally, the value of $\hat{\beta}_{OLS,2}$ appears to affect the estimate $\hat{\beta}_{1}$ more when $\hat{\beta}_{OLS,1}$ is smaller under the SPN prior, SNG prior with $c\leq 1$ and SPB prior with $c\leq 1$. This is consistent with what we observed when examining the copulas and conditional prior distributions; as the tails become heavier, the Laplace and heavier-than Laplace prior distributions display more dependence at the origin than the tails. Last, the SNG prior with $c=0.3$ not only induces especially strong shrinkage of $\hat{\beta}_{1}$ when $\hat{\beta}_{OLS,2}$ is small, but also induces aggressive inflation of $\hat{\beta}_{1}$ when $\hat{\beta}_{OLS,1}$ is small and $\hat{\beta}_{OLS,2}$ is very large. This makes sense in the context of the prior conditional distributions under the SNG and SPB priors with $c=0.3$ and $q=0.65$. Although these two priors have equally heavy tails, the SNG prior is much more concentrated at the origin and has heavier tails.

4 Computation

4.1 Posterior Approximation

The posterior mode of $\boldsymbol{\beta}$ maximizes the integral $\int p(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\beta},\boldsymbol{\phi})p(\boldsymbol{\beta}|\boldsymbol{\Omega},\boldsymbol{s})p(\boldsymbol{s}|\boldsymbol{\theta})d\boldsymbol{s}$ over $\boldsymbol{\beta}$ and can be thought of as a penalized estimate of $\boldsymbol{\beta}$, where the penalty is given by $-\text{log}(\int p(\boldsymbol{\beta}|\boldsymbol{\Omega},\boldsymbol{s})p(\boldsymbol{s}|\boldsymbol{\theta})d\boldsymbol{s})$. This integral is generally intractable when $\boldsymbol{\Omega}$ is not diagonal, but can be maximized using an MCMC within Expectation Maximization (EM) algorithm (Dempster et al., 1977). Given an initial value $\boldsymbol{\beta}^{(0)}$, this algorithm proceeds by iterating the following until $||\boldsymbol{\beta}^{(i+1)}-\boldsymbol{\beta}^{(i)}||^{2}_{2}$ converges:

• •

using MCMC to simulate $M$ draws $\boldsymbol{s}^{(1)},\dots,\boldsymbol{s}^{(M)}$ from the full conditional distribution of $\boldsymbol{s}$ given $\boldsymbol{\beta}^{(i)}$ and $\boldsymbol{\theta}$, set $\widehat{\mathbb{E}}[(\frac{\boldsymbol{1}}{\boldsymbol{s}})(\frac{\boldsymbol{1}}{\boldsymbol{s}})^{\prime}|\boldsymbol{\beta}^{(i)},\boldsymbol{\theta}]=\frac{1}{M}\sum_{j=1}^{M}(\frac{\boldsymbol{1}}{\boldsymbol{s}^{(j)}})(\frac{\boldsymbol{1}}{\boldsymbol{s}^{(j)}})^{\prime}$;

• •

set $\boldsymbol{\beta}^{(i+1)}=\text{argmin}_{\boldsymbol{\beta}}h(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\beta},\boldsymbol{\phi})+\frac{1}{2}\boldsymbol{\beta}^{\prime}(\boldsymbol{\Omega}^{-1}\circ\widehat{\mathbb{E}}[(\frac{\boldsymbol{1}}{\boldsymbol{s}})(\frac{\boldsymbol{1}}{\boldsymbol{s}})^{\prime}|\boldsymbol{\beta}^{(i)},\boldsymbol{\theta}])\boldsymbol{\beta}$.

Alternative posterior summaries can be obtained by simulating $M$ values $\boldsymbol{\beta}^{(1)},\dots,\boldsymbol{\beta}^{(M)}$ from the joint posterior distribution of $(\boldsymbol{\beta},\boldsymbol{s})$ given $\boldsymbol{y}$, $\boldsymbol{X}$, $\boldsymbol{\phi}$ and $\boldsymbol{\theta}$, starting with an initial value $\boldsymbol{\beta}^{(0)}$ and iteratively simulating $\boldsymbol{s}^{(j)}$ from the full conditional distribution of $\boldsymbol{s}$ given $\boldsymbol{\beta}^{(j-1)}$ and $\boldsymbol{\theta}$ and simulating $\boldsymbol{\beta}^{(j)}$ from the full conditional distribution of $\boldsymbol{\beta}$ given $\boldsymbol{y}$, $\boldsymbol{X}$, $\boldsymbol{s}^{(j)}$ and $\boldsymbol{\phi}$. When the log-likelihood is quadratic or conditionally quadratic in $\boldsymbol{\beta}$, the full conditional distribution $p(\boldsymbol{\beta}|\boldsymbol{s},\boldsymbol{y},\boldsymbol{X},\boldsymbol{\phi})$ is a multivariate normal distribution.

The full conditional distribution $p(\boldsymbol{s}|\boldsymbol{\beta},\boldsymbol{\Omega},\boldsymbol{\theta})$ can be written as proportional to

[TABLE]

where ‘$/$’ is applied elementwise. The choices $p(\boldsymbol{s}|\boldsymbol{\theta})$ that yield the SPN, SNG and SPB models do not yield standard distributions when $\boldsymbol{\Omega}$ is not a diagonal matrix.

We simulate from the full conditional distribution (3) using generalized elliptical slice sampling (Nishihara et al., 2014). First, we unify all three models by recognizing that simulating $\boldsymbol{s}$ according to (3) is equivalent to simulating unconstrained $\boldsymbol{s}^{\prime}$ according to

[TABLE]

and setting $\boldsymbol{s}=f\left(\boldsymbol{s}^{\prime}\right)$, where $f\left(\boldsymbol{s}^{\prime}\right)=\boldsymbol{s}^{\prime}$ when the SPN prior is used and $f\left(\boldsymbol{s}^{\prime}\right)=\left|\boldsymbol{s}^{\prime}\right|$ when the SNG or SPB priors are used. Such use of the absolute value transformation is common in the generalized elliptical slice sampling literature (Nishihara et al., 2014). Then

we perform a change of variables from \boldsymbol{s}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0} to \boldsymbol{s}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\boldsymbol{V}^{-1/2}\left(\boldsymbol{s}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-\boldsymbol{m}\right), where $\boldsymbol{m}$ and $\boldsymbol{V}$ may depend on the current values of $\boldsymbol{\beta}$, $\boldsymbol{\Omega}$, and $\boldsymbol{\theta}$ but not $\boldsymbol{s}$ and $\boldsymbol{V}^{-1/2}$ refers to the symmetric square root of the pseudoinverse $\boldsymbol{V}^{\dagger}$ of $\boldsymbol{V}$. Applying a change of variables when simulating from conditional distributions is a standard practice, as seen in (Murray et al., 2010).

Simulating $\boldsymbol{s}$ according to (3) is equivalent to simulating \boldsymbol{s}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0} according to

[TABLE]

and setting $\boldsymbol{s}=f\left(\boldsymbol{V}^{1/2}\boldsymbol{s}^{\prime\prime}+\boldsymbol{m}\right)$,

where ‘$/$’ and $\left|\cdot\right|$ are applied elementwise. In turn, simulating \boldsymbol{s}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0} according to (4) is equivalent to simulating four $p\times 1$ vectors $\boldsymbol{r}$, $\boldsymbol{u}$, and $\boldsymbol{\pi}$ according to

[TABLE]

For fixed $\boldsymbol{m}$ and $\boldsymbol{V}$, iteration $i$ of a Gibbs sampler simulating from (5) is as follows:

Simulate $t^{(i)}_{j}\sim\text{normal}(0,1)$ for $j=1,\dots,p$. 2. 2.

Set \boldsymbol{u}^{(i)}=(\boldsymbol{s}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{\prime\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{(i-1)})\text{sin}(\boldsymbol{\pi}^{(i-1)})+\boldsymbol{t}^{(i)}\text{cos}(\boldsymbol{\pi}^{(i-1)}). 3. 3.

Set \boldsymbol{r}^{(i)}=(\boldsymbol{s}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{\prime\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{(i-1)})\text{cos}(\boldsymbol{\pi}^{(i-1)})-\boldsymbol{t}^{(i)}\text{sin}(\boldsymbol{\pi}^{(i-1)}). 4. 4.

Simulate $\pi^{(i)}_{j}$ according to (5) for $j=1,\dots,p$ using univariate slice sampling. 5. 5.

Set \boldsymbol{s}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{\prime\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{(i)}=\boldsymbol{u}^{(i)}\text{sin}(\boldsymbol{\pi}^{(i)})+\boldsymbol{r}^{(i)}\text{cos}(\boldsymbol{\pi}^{(i)}). 6. 6.

Set \boldsymbol{s}^{(i)}=\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f\left(\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{V}^{1/2}\boldsymbol{s}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{(i)}+\boldsymbol{m}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\right)\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}.

Univariate slice sampling is described in the appendix. We can think of steps 1-4 as generating a proposal for \boldsymbol{s}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{(i)} and step 5 as generating an adjustment for the proposal that reflects the relationship between the proposal distribution and the target full conditional distribution. Validity of the proposed algorithm follows from Murray et al. (2010). The vector $\boldsymbol{m}$ and matrix $\boldsymbol{V}$ can be chosen to improve mixing and decrease computational burden. Choices of $\boldsymbol{m}$ and $\boldsymbol{V}$ that provide better approximations to the full conditional distribution $p(\boldsymbol{s}|\boldsymbol{\beta},\boldsymbol{\Omega},\boldsymbol{\theta})$ are likely to result in better mixing, whereas choices of $\boldsymbol{m}$ and $\boldsymbol{V}$ that are easier to compute will improve the speed of the sampling process. Here, $\boldsymbol{m}$ and $\boldsymbol{V}$ are chosen to approximate the mode of p(\boldsymbol{s}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}|\boldsymbol{\beta},\boldsymbol{\Omega},\boldsymbol{\theta}) and corresponding Hessian. We compute $\boldsymbol{m}$, an approximate mode of p(\boldsymbol{s}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{\prime}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}|\boldsymbol{\beta},\boldsymbol{\Omega},\boldsymbol{\theta}), by performing coordinate descent with a large convergence threshold and a small number of maximum iterations. We describe how $\boldsymbol{m}$ can be obtained via coordinate descent under the SHP priors in the appendix.

4.1.1 Simulating from the Posterior Distribution under the SPN Prior

As noted in Section 2, simulation from the joint posterior distribution under the SPN prior is simple. In the linear regression setting with known unit variance where $-h(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\beta})\propto_{\boldsymbol{\beta}}\frac{1}{2}(\boldsymbol{\beta}^{\prime}\boldsymbol{X}^{\prime}\boldsymbol{X}\boldsymbol{\beta}-2\boldsymbol{\beta}^{\prime}\boldsymbol{X}^{\prime}\boldsymbol{y})$, then

[TABLE]

Both full conditional distributions are multivariate normal distributions. As a result, a straightforward Gibbs sampler can be used to simulate from the joint posterior distribution.

4.2 Hyperparameter Estimation

In the previous subsection, we presented a general approach for simulating from the posterior distribution of $\boldsymbol{\beta}$ under the SPN, SNG and SPB priors given hyperparameters $\boldsymbol{\Omega}$ and, in the case of the SPN prior, $\boldsymbol{\Psi}$. When the hyperparameters are unknown, one approach is to assume prior distributions for $\boldsymbol{\phi}$, $\boldsymbol{\Omega}$, and/or $\boldsymbol{\Psi}$. For all three priors, a conjugate inverse-Wishart prior for $\boldsymbol{\Omega}$ is a natural choice. For the SPN prior, a conjugate inverse-Wishart prior for $\boldsymbol{\Psi}$ is likewise natural.

There are situations where conjugate inverse-Wishart priors may have undesirable properties (Schuurman et al., 2016), in which case alternative prior distributions for $\boldsymbol{\Omega}$ and, when the SPN prior is used, $\boldsymbol{\Psi}$, could be used.

Alternatively, hyperparameter estimates can be obtained via maximum marginal likelihood estimation (MMLE) or the method of moments. We can compute the MMLE of the unknown variance components $\boldsymbol{\Omega}$ and, in the case of the SPN prior, $\boldsymbol{\Psi}$, using a Gibbs-within-EM algorithm as described in the appendix. However, maximum marginal likelihood estimation of hyperparameters can converge prohibitively slowly in practice (Roy and Chakraborty, 2016). Furthermore, the Gibbs step can be prohibitively computationally demanding when the data are high dimensional. Fortunately, method of moments type estimates of the unknown variance components can be obtained under the SNG and SPB priors for fixed $c$ and $q$ respectively, and under the symmetric sSPN prior, so long as $\boldsymbol{y}$ is linearly related to $\boldsymbol{X}\boldsymbol{\beta}$. As noted in the Introduction, the prior moments are easy to compute under all three priors

[TABLE]

Furthermore, under the sSPN, SNG and SPB priors, the hyperparameters are correspond to second order moments of $\boldsymbol{\beta}$. When $y$ is linearly related to $\boldsymbol{X}\boldsymbol{\beta}$, a positive semidefinite estimate of $\boldsymbol{\Sigma}$ can be obtained using methods from Perry (2017). Under the sSPN prior, estimates of $\boldsymbol{\Omega}$ and $\boldsymbol{\Psi}$ can be obtained from an estimate $\hat{\boldsymbol{\Sigma}}$ by projecting $\sqrt{|\hat{\boldsymbol{\Sigma}}|}$ and $\text{sign}\{\hat{\boldsymbol{\Sigma}}\}\sqrt{|\hat{\boldsymbol{\Sigma}}|}$ onto the positive semi-definite cone, where $\sqrt{\cdot}$, $|\cdot|$ and $\text{sign}\{\cdot\}$ are applied elementwise. Under the SNG and SPB priors, an estimate of $\boldsymbol{\Omega}$ can be obtained from an estimate $\hat{\boldsymbol{\Sigma}}$ by projecting $\hat{\boldsymbol{\Sigma}}/((1-\mathbb{E}[s_{j}]^{2})\boldsymbol{I}_{p}+\mathbb{E}[s_{j}]^{2}\boldsymbol{1}_{p}\boldsymbol{1}_{p}^{\prime})$ onto the positive semi-definite cone, where ‘$/$’ is applied elementwise, $\mathbb{E}[s_{j}]^{2}=c^{-1}(\Gamma(c+1/2)/\Gamma(c))^{2}$ under the SNG prior and $\mathbb{E}[s_{j}]^{2}=(\pi/2)(\Gamma(2/q)/\sqrt{\Gamma(1/q)\Gamma(3/q)})^{2}$ under the SPB prior.

4.3 Prior Selection

Our focus is on understanding the properties of each novel prior and providing a method for using them in practice. However, we would be remiss to omit guidance on when SNG, SPB, and SPN priors should be used as opposed to more common existing alternatives, how to choose between the SNG, SPB, and SPN priors when their use is warranted, and, in the case of the SNG and SPB priors, how to choose $c$ or $q$.

Ideally, our understanding of the scientific problem would inform our choice of prior. Alternatively, we recommend first choosing a criterion for selecting a prior that is appropriate to the specific application and goals of analysis, e.g. Deviance Information Criterion (DIC), the widely applicable information criterion (WAIC), log marginal likelihood, mean log conditional predictive ordinate (CPO)/leave-one-out cross validation error (LOO) (Roos and Held, 2011; Piironen and Vehtari, 2017; Vehtari et al., 2017). We then recommend performing an initial analysis that compares the chosen criterion under an independent normal, multivariate normal, and independent sparsity inducing prior of interest, e.g. a Laplace prior. If both the multivariate normal and independent sparsity inducing prior provide improved fit compared to an independent normal prior, we recommend considering the SPN, SNG and/or SPB priors for a range of values of $c$ and/or $q$ and selecting the prior that is best according to the chosen criterion.

5 Application

In this section we return to the subset of BCI data discussed in Section 1. We assume the following model for the indictors of whether or not the subject was viewing the target letter $\boldsymbol{y}$ and EEG measurements at time point $j$ and channel $k$, $\boldsymbol{x}_{jk}$ during the first 20 trials,

[TABLE]

where each $\gamma_{jk}$ is a time point and channel specific intercept and each $\beta_{jk}$ is a regression coefficient that describes the relationship between EEG measurements from time point $j$ and channel $k$ and whether or not the subject was viewing the target letter, for $j=1,\dots,p_{1}$ and $k=1,\dots,p_{2}$. Because simulating from the posterior distribution of $\boldsymbol{\beta}$ becomes increasingly computationally demanding as the number of unknown regression coefficients $\beta_{jk}$ increases and because we want to explore the results of using several different SNG and SPB priors, we consider a reduced set of time points, specifically every eighth time point. This reduces $p_{1}$ from $208$ to $26$. We continue to consider all channels, $p_{2}=8$.

Letting $\boldsymbol{\beta}=\text{vec}\left(\boldsymbol{B}\right)$ refer to a vectorized $p_{1}\times p_{2}$ matrix $\boldsymbol{B}$ with elements $b_{jk}=\beta_{jk}$, we assume that $\mathbb{E}\left[\boldsymbol{\beta}\right]=\boldsymbol{0}$ and $\mathbb{V}\left[\boldsymbol{\beta}\right]=\boldsymbol{\Sigma}$, where $\boldsymbol{\Sigma}$ is separable, i.e. $\boldsymbol{\Sigma}=\boldsymbol{\Sigma}_{2}\otimes\boldsymbol{\Sigma}_{1}$. The covariance matrices $\boldsymbol{\Sigma}_{1}$ and $\boldsymbol{\Sigma}_{2}$ characterize relationship of regression coefficients $\boldsymbol{B}$ over time and across channels, respectively. Because the scales of $\boldsymbol{\Sigma}_{1}$ and $\boldsymbol{\Sigma}_{2}$ are not separately identifiable, we assume that the diagonal elements of $\boldsymbol{\Sigma}_{1}$ are exactly equal to 1. We also assume that $\boldsymbol{\Omega}_{1}$ and and $\boldsymbol{\Psi}_{1}$ have autoregressive structures of order one, with $\omega_{1,ij}=\rho^{\left|i-j\right|}_{\boldsymbol{\Omega}}$ and $\psi_{1,ij}=\rho^{\left|i-j\right|}_{\boldsymbol{\Psi}}$. We note that the corresponding marginal variance $\boldsymbol{\Sigma}_{1}$ is autoregressive of order one under the SPN prior but not the SNG and SPB priors. If $\boldsymbol{\Omega}_{1}$ and $\boldsymbol{\Psi}_{1}$ have autoregressive structure of order one with parameters $\rho_{\boldsymbol{\Omega}}$ and $\rho_{\boldsymbol{\Psi}}$, then $\boldsymbol{\Sigma}_{1}$ has autoregressive structure of order one with parameter $\rho=\rho_{\boldsymbol{\Omega}}\rho_{\boldsymbol{\Psi}}$. In contrast, the matrix $\boldsymbol{\Sigma}_{1}=\boldsymbol{\Omega}_{1}\circ\mathbb{E}\left[\boldsymbol{s}\boldsymbol{s}^{\prime}\right]$ does not have autoregressive structure of order one.

We assume $\rho_{\boldsymbol{\Omega}}\sim\text{beta}\left(\left(p_{1}+1\right)/2,\left(p_{1}+1\right)/2\right)$ and $\boldsymbol{\Omega}^{-1}_{2}\sim\text{Wishart}(p_{2}+2,\kappa^{-1}\boldsymbol{I}_{p_{2}})$, where $\kappa=\frac{1}{p_{1}p_{2}-1}\sum_{j=1}^{p_{1}}\sum_{k=1}^{p_{2}}\left(\hat{\beta}_{logit,jk}-\bar{\hat{\beta}}_{logit}\right)^{2}$, $\bar{\hat{\beta}}_{logit}=\frac{1}{p_{1}p_{2}}\sum_{j=1}^{p_{1}}\sum_{k=1}^{p_{2}}\hat{\beta}_{logit,jk}$, and $\hat{\beta}_{logit,jk}$ are the logistic regression estimates of $\beta_{jk}$ depicted in Figure 1. When using the SPN prior, we also assume $\boldsymbol{\Psi}^{-1}_{2}\sim\text{Wishart}(p_{2}+2,\boldsymbol{I}_{p_{2}})$ and $\rho_{\boldsymbol{\Psi}}\sim\text{beta}\left(\left(p_{1}+1\right)/2,\left(p_{1}+1\right)/2\right)$. We assume conjugate inverse-Wishart priors for convenience, as our primary goal is to compare estimates under the SNG, SPB, and SPN priors for $\boldsymbol{\beta}$. We use the latent variable representation of the logistic regression model introduced in Polson et al. (2013), simulate from the full conditional distribution of $\boldsymbol{\beta}$ jointly, use the elliptical slice sampling procedure described in Section 4 to simulate from the full conditional distribution of $\boldsymbol{s}$, and use univariate slice sampling to simulate from the full conditional distributions of $\rho_{\boldsymbol{\Omega}}$ and $\rho_{\boldsymbol{\Psi}}$. We simulate 20 chains of 13,500 samples from the posterior distribution, discard the initial $1,000$ from each chain as burn-in, and thin by a factor of $5$. The remaining 50,000 samples have minimum effective sample sizes shown in the appendix.

Figure 8 shows posterior mean estimates of elements of $\boldsymbol{B}$ under a multivariate normal prior (SNO) and nine SHP priors. Because shrinkage priors are often used to select nonzero elements of $\boldsymbol{B}$, approximate 90% intervals that do not include zero are also indicated. In cases where a sparse point estimate of $\boldsymbol{B}$ is desired, the methods of Hahn and Carvalho (2015), Li and Pati (2017), Piironen et al. (2020), or Woody et al. (2021) can alternatively be used. The estimated posterior means and approximate 90% intervals that do not include zero display similar structure across channels and over time. They also show increasing shrinkage of individual elements of $\boldsymbol{\beta}$ as $c$ or $q\rightarrow 0$. Even the SNG and SPB priors with $c=0.3$ and $q=0.65$, which encourage sparsity very aggressively, produce estimates are structured, i.e. estimates of elements of $\boldsymbol{B}$ corresponding to the same channel or time point tend to be similar regardless of whether they are nearly equal to zero or large in magnitude. In comparison to the normal prior, the more aggressive sparsity inducing priors facilitate identification of especially strong signals; for instance, they suggest that whether or not the subject is viewing the target letter is strongly correlated with EEG measurements at the fifth and twelfth time points on nearly all of the channels.

Figure 9 shows posterior mean estimates of channel-by-channel correlation matrices corresponding to $\boldsymbol{\Sigma}_{2}$ and indications of approximate 90% intervals that do not include zero as well as posterior mean estimates of correlations across consecutive time points $\sigma_{1,i\left(i+1\right)}$ and approximate 90% intervals. There is evidence of dependence, especially across channels. The amount of dependence across channels and over time under the SNG and SPB priors is decreasing as $c$ or q\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\rightarrow 0. Surprisingly, the SPN prior estimates weak correlations despite being able to accommodate arbitrary $\boldsymbol{\Sigma}$. This may be due to especially strong shrinkage of elements of $\boldsymbol{\Sigma}$ under the assumed prior distributions for the hyperparameters.

The available prior information, which suggests that $\boldsymbol{B}$ is sparse with strong signals observed across several channels at a subset of time points that roughly correspond to the expected timing of the P300 wave, can be used to assess the relative merits of the different priors. The SNG and SPB priors with $c=10$ and $q=1.75$, respectively, produce estimates that are not sparse and so similar to the estimates produced under a normal prior that there is no added benefit relative to a multivariate normal prior. Meanwhile, the SNG prior with $c=0.3$ and the SPN prior produce correlations that are not distinguishable from zero and thus provide no added benefit relative to an independent sparsity inducing prior and are implausibly sparse. This leaves the SNG and SPB priors with $c=1$ and $q=1$, which are equivalent and generalize the independent Laplace prior, the SNG prior with $c=0.5$ and the SPB priors with $q=0.78$, and $q=0.65$. The SNG prior with $c=0.5$ misses signals that occur later in each trial that are likely to correspond to the known timing of the P300 wave. The SPB priors with $q=0.78$ and $q=0.65$ provide estimates that are qualitatively very similar to the estimates obtained under the SPB prior with $q=1$. Accordingly, we recommend the multivariate Laplace prior, which corresponds to the SNG prior with $c=1$ and the SPB prior with $q=1$; it is the simplest SHP prior that produces estimates that align with the available scientific context.

For the purposes of demonstration, we also consider three criteria for prior selection; WAIC, DIC, and LOO. Approximations to these criteria are shown in Figure 10. Regardless of the criteria considered, the sparsity inducing independent product normal prior and the structured multivariate normal prior both outperform the the independent normal prior and therefore suggests that SHP priors are worth considering. All structured priors outperform their independent counterparts. A detailed comparison of the estimates of $\boldsymbol{B}$ obtained under the structured priors we consider and their independent counterparts is included in the appendix. Among structured priors, less sparsity inducing priors perform better and the multivariate normal prior performs best in this specific setting according to these criteria.

To conclude, we compare approximate 90% intervals that do not include zero to “ground truth” 90% confidence intervals that do not include zero based on simple logistic regression estimates computed from all 240 trials. Table 2 shows that normal and nearly normal priors provide the best true positive rates, while the multivariate Laplace prior and SPB priors with $q<1$ provide the best true positive rates relative to false positive rates.

6 Discussion

We introduce the SHP class of novel structured shrinkage priors for regression coefficients and show that they can encourage both sparsity and structure, which can be difficult to simultaneously model using existing prior distributions. We provide a parsimonious and general approach to simulation from the posterior distributions under SHP priors based on elliptical slice sampling and demonstrate how SHP priors can improve interpretability of estimated regression coefficients relative to multivariate normal or independent priors.

We have focused on the development of structured shrinkage prior distributions for regression coefficients, however the same distributions could be used to model errors as in Finegold and Drton (2011) or as alternatives to Gaussian processes. We could also construct shrinkage priors for covariance matrices by extending Daniels and Pourahmadi (2002) or building on the matrix-$F$ prior distribution for covariance matrices introduced in Mulder and Pericchi (2018), which assumes that $\boldsymbol{\Sigma}\stackrel{{\scriptstyle d}}{{=}}\boldsymbol{S}\boldsymbol{\Omega}\boldsymbol{S}^{\prime}$, where $\boldsymbol{\Omega}^{-1}$ has an inverse-Wishart distribution with scale $\boldsymbol{I}_{p}$ and degrees of freedom $\delta+k-1$ and $\boldsymbol{S}\boldsymbol{S}^{\prime}$ has a Wishart distribution with scale $\boldsymbol{\Xi}$ and degrees of freedom $\nu$. We could define distributions for covariance matrices according to $\boldsymbol{\Sigma}\stackrel{{\scriptstyle d}}{{=}}\text{diag}\left\{\boldsymbol{s}\right\}\boldsymbol{\Omega}\text{diag}\left\{\boldsymbol{s}\right\}$ for $\boldsymbol{s}\sim\text{normal}\left(\boldsymbol{0},\boldsymbol{\Psi}\right)$, $s^{2}_{j}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\text{gamma}(c,c)$, or $s^{2}_{j}$ independently distributed according to a polynomially tilted positive $\alpha$-stable distribution with index of stability $\alpha=q/2$ and $\mathbb{E}[s^{2}_{j}]=1$. Also, the elliptical slice sampling procedure we use to construct a tractable Gibbs sampler for simulating from $p(\boldsymbol{s}|\boldsymbol{\beta},\boldsymbol{\theta})$ under the SPN, SNG and SPB priors can be to perform posterior inference under other novel structured generalizations of other shrinkage priors that can be represented using Hadamard products involving a normal random vector, including the horseshoe and Dirichlet-Laplace priors (Polson and Scott, 2010; Bhattacharya et al., 2015). Popularity of the horseshoe prior suggests that a structured generalization may be especially valuable. Additional extensions might be to assess whether or not conditions for the convergence of elliptical slice sampling algorithms provided in Natarovskii et al. (2021) hold for the proposed elliptical slice sampling procedure and to improve scalability of the proposed elliptical slice sampling procedure with respect to $p$, the total number of penalized covariates. As shown in the appendix, seconds per posterior draw increases and minimum effective sample size decreases under these structured priors as more penalized covariates are used. The time costs of using structured priors as opposed to their independent counterparts increases with the number of penalized covariates and is greater for priors that encourage sparsity more aggressively.

Another extension to this work would treat $c$ and $q$ as unknown under the SNG and SPB priors, respectively, and perform maximum marginal likelihood estimation of $c$ or $q$ or assume prior distributions for $c$ and $q$.

Last, an extension might consider modifications of the proposed priors that continue to allow for structure in the presence of extreme elementwise sparsity, e.g. extensions of SNG and SPB priors that allow for correlated scales and extensions of the SPN prior that can encourage shrinkage more aggressively by considering elementwise products of more than two normal vectors.

Acknowledgements

This work was partially supported by NSF grants DGE-1256082 and DMS-1505136.

Supplementary Materials

Supplementary material available online includes proofs of all the propositions and additional numerical results. A stand-alone package for implementing the methods described in this paper can be downloaded from https://github.com/maryclare/sspcomp.

Appendix A Exploratory Data

Appendix B Relationship Between $\boldsymbol{\Omega}$, $\boldsymbol{\Psi}$ and Fourth-Order Prior Cross Moments

Let $\boldsymbol{C}_{\boldsymbol{\Omega}}$ be the correlation matrix corresponding to $\boldsymbol{\Omega}$, and recall that parametrize $\boldsymbol{\Psi}$ as having diagonal elements equal to $1$. Off-diagonal elements of $\boldsymbol{C}_{\boldsymbol{\Omega}}$ and $\boldsymbol{\Psi}$ determine fourth-order cross moments of elements of $\boldsymbol{\beta}$. Recalling that the marginal covariance matrix $\boldsymbol{\Sigma}=\boldsymbol{\Omega}\circ\boldsymbol{\Psi}$ and letting $\boldsymbol{C}_{\Sigma}$ be the correlation matrix corresponding to $\boldsymbol{\Sigma}$, we have

[TABLE]

We can see that these fourth-order cross moments depend on the values of both $c_{\Omega,jk}$ and $c_{\Psi,jk}$, in addition to their product $c_{\Sigma,jk}$.

Appendix C Generalized Gamma Rate Mixture Representation of Polynomially Tilted Positive $\alpha$-Stable Random Variables

Following Devroye (2009), we can write polynomially titled positive $\alpha$-stable $s^{2}_{j}$ as:

[TABLE]

where $f(\delta_{j}|\alpha)=\text{sin}(\alpha\delta_{j})^{\frac{\alpha}{1-\alpha}}\text{sin}((1-\alpha)\delta_{j})/\text{sin}(\delta_{j})^{\frac{1}{1-\alpha}}$ and $\delta_{j}\in(0,\pi)$. The density $f(\delta_{j}|\alpha)$ is non-standard, however Devroye (2009) provides a following method for simulating from $\delta_{j}$ from this distribution.

Appendix D Univariate Marginal Distributions

Intuitively, it is clear from the stochastic representation that the marginal distributions are the same as the corresponding univariate shrinkage prior. We can show this directly as follows. The joint marginal prior distribution of $\boldsymbol{\beta}$ is

[TABLE]

Then $p(\beta_{1})$ is given by

[TABLE]

The term $(*)$ is equal to $\int p(\boldsymbol{\beta}_{-1}|\beta_{1}/s_{1})d\beta_{2}\dots d\beta_{p}$. This is the integral of a density, and accordingly $(*)=1$ and $p(\beta_{1})=\int p(\beta_{1}/s_{1})p(s_{1})/|s_{1}|ds_{1}$.

Appendix E Proofs of Propositions

E.1 Propositions 2.1 and 2.3

First, we prove the following lemma:

Lemma E.1

For $\alpha>0$ and $\gamma\in\mathbb{R}$, $\int_{-\infty}^{\infty}\frac{1}{|s|}\text{exp}\{-\alpha(s^{2}-\gamma s)\}ds=+\infty$.

First let’s consider this integral when $\gamma>0$:

[TABLE]

Because the integrand is nonnegative for all $s$, if any of these terms evaluate to $+\infty$, the entire integral evaluates to $+\infty$. Let’s examine the middle integral. Note that over this range, $s-\gamma\leq 0$, $\text{exp}\{-\alpha s\}\leq 1$ and accordingly, $\text{exp}\{-\alpha s\}^{s-\gamma}\geq 1$.

[TABLE]

Now let’s consider the same integral when $\gamma<0$:

[TABLE]

Again, if any of these terms evaluate to $+\infty$, the entire integral evaluates to $+\infty$. Again, let’s consider the middle term. Over this interval, $\text{exp}\{-\alpha s\}\geq 1$ and $s-\gamma\geq 0$. It follows that $\text{exp}\{-\alpha s\}^{s-\gamma}\geq 1$ and:

[TABLE]

Now we’ll consider one last case where $\gamma=0$:

[TABLE]

As in the previous cases, any of these terms evaluate to $+\infty$, the entire integral evaluates to $+\infty$. Examining the middle term one last time, we have:

[TABLE]

E.1.1 Proof of Proposition 2.1:

Now we can evaluate the marginal density $p(\boldsymbol{\beta}=\boldsymbol{b}|\boldsymbol{\Omega},\boldsymbol{\Psi})$ when $b_{j}=0$ for some $j\in\{1,\dots,p\}$. Without loss of generality, set $\beta_{1}=0$. Letting $\boldsymbol{b}_{-1}$ and $\boldsymbol{s}_{-1}$ refer to the vectors $\boldsymbol{b}$ and $\boldsymbol{s}$ each with the first element removed, $(\boldsymbol{\Omega}^{-1})_{-1,-1}$ be the matrix $\boldsymbol{\Omega}^{-1}$ with the first row and column removed, $(\boldsymbol{\Psi}^{-1})_{-1,-1}$ be the matrix $\boldsymbol{\Psi}^{-1}$ with the first row and column removed and $(\boldsymbol{\Psi}^{-1})_{-1,1}$ be the first column of $\boldsymbol{\Psi}^{-1}$ excluding the first element. For $\boldsymbol{b}=(0,\boldsymbol{b}_{-1})$,

[TABLE]

Applying Lemma E.1, the term denoted by $(*)$ evaluates to $+\infty$ for every value of $\boldsymbol{s}_{-1}$.

E.1.2 Proof of Proposition 2.3:

Proposition 2.3 follows from Proposition 2.1 and Lemma E.1. Again, letting $\beta_{1}=0$, for $\boldsymbol{\beta}=(0,\boldsymbol{\beta}_{-1})$,

[TABLE]

Again, $(*)$ evaluates to $+\infty$ for any value of $\boldsymbol{s}_{-1}$ and does not depend at all on $\boldsymbol{\beta}_{-1}$.

E.2 Proofs of Propositions 2.2 and 2.4

First, we prove the following lemma:

Lemma E.2

For $0<c<1/2$, $\int_{0}^{\infty}(s^{2})^{c-1/2-1}\text{exp}\{-cs^{2}\}ds^{2}=+\infty$.

We can break the integral into two nonnegative components, as follows:

[TABLE]

Now let’s examine the first component. When $s^{2}<1/c$, $\text{exp}\{-cs^{2}\}\geq 1-cs^{2}$ and

[TABLE]

E.2.1 Proof of Proposition 2.2

Now we can evaluate the marginal density $p(\boldsymbol{\beta}=\boldsymbol{b}|c,\boldsymbol{\Omega})$ when $b_{j}=0$ for some $j\in\{1,\dots,p\}$. Without loss of generality, set $\beta_{1}=0$. Letting $\boldsymbol{b}_{-1}$ and $\boldsymbol{s}_{-1}$ refer to the vectors $\boldsymbol{b}$ and $\boldsymbol{s}$ each with the first element removed and $(\boldsymbol{\Omega}^{-1})_{-1,-1}$ be the matrix $\boldsymbol{\Omega}^{-1}$ with the first row and column removed. For $\boldsymbol{b}=(0,\boldsymbol{b}_{-1})$,

[TABLE]

Applying Lemma E.2, the term denoted by $(*)$ evaluates to $+\infty$ for every value of $\boldsymbol{s}_{-1}$.

E.2.2 Proof of Proposition 2.4

Proposition 2.4 follows from Proposition 2.2 and Lemma E.2. For $\boldsymbol{\beta}=(0,\boldsymbol{\beta}_{-1})$,

[TABLE]

Again, $(*)$ evaluates to $+\infty$ for any value of $\boldsymbol{s}_{-1}$ and does not depend at all on $\boldsymbol{\beta}_{-1}$.

Appendix F Kurtosis of SHP $\beta_{j}$

F.1 SNG $\beta_{j}$

When $s^{2}_{j}\sim\text{gamma}(c,c)$, we have:

[TABLE]

A standard normal random variable has $\mathbb{E}[z^{4}_{j}]=3$. It follows that

[TABLE]

F.2 SPB $\beta_{j}$

The kurtosis of a random variable is not a function of the overall scale, so without loss of generality, let $\mathbb{V}[\beta_{j}]=1$. When $\boldsymbol{\beta}$ is distributed according to an $SPB$ prior, elements $\beta_{j}$ each have a power/bridge distribution with density

[TABLE]

It follows that

[TABLE]

It follows that the kurtosis of $\beta_{j}$ is

[TABLE]

Appendix G Expectation of $s_{j}$

G.1 SNG $\beta_{j}$ with $\mathbb{V}[\beta_{j}]=1$

When $s^{2}_{j}\sim\text{gamma}(c,c)$, we have:

[TABLE]

G.2 SPB $\beta_{j}$ with $\mathbb{V}[\beta_{j}]=1$

This is a little challenging because working with the stable distribution directly is difficult. However, given knowledge of normal moments and the marginal distribution of $\beta$, we can “back out” $\mathbb{E}[s]$. Let $\beta=sz$, where $z$ is a standard normal random variable and $1/s^{2}$ has an $\alpha$-stable distribution on the positive real line. When the stable prior is parametrized to yield $\mathbb{E}[\beta^{2}]=1$, we have:

[TABLE]

Now, to determine $\mathbb{E}[s]$, note that $\mathbb{E}[|\beta|]=\mathbb{E}[s]\mathbb{E}[|z|]$. We have

[TABLE]

According to Wikipedia, $\mathbb{E}[|z|]=\sqrt{2/\pi}$. It follows that

[TABLE]

Appendix H Prior Conditional Distributions of $s^{2}_{j}$ and $s_{j}$

H.1 SPB $\beta_{j}$

When using the SPB prior, we have $s^{2}_{j}\stackrel{{\scriptstyle d}}{{=}}\Gamma(1/(2\alpha))\xi_{j}^{\frac{1-\alpha}{\alpha}}/(2\Gamma(3/(2\alpha)))$, where $\xi_{j}\sim\text{gamma}((1+\alpha)/(2\alpha),f(\delta_{j}|\alpha))$ and $p(\delta_{j}|\alpha)\propto f(\delta_{j}|\alpha)^{\frac{\alpha-1}{2\alpha}}$ on $(0,\pi)$. We want to compute the prior conditional distribution $\pi(s^{2}_{j}|\delta_{j},\alpha)$ for use in the elliptical slice sampler for $\boldsymbol{s}$. We have:

[TABLE]

We can find the density of $s^{2}_{j}$ via change of variables. We have

[TABLE]

Then the prior conditional distribution is

[TABLE]

In practice, we end up working with $s_{j}$, so once more

[TABLE]

H.2 SNG Prior $\beta_{j}$

Under the SNG prior, we have

[TABLE]

In practice, we end up working with $s_{j}$, which has the density

[TABLE]

Appendix I Maximum correlation for bivariate SHP $\boldsymbol{\beta}$

I.1 SNG $\boldsymbol{\beta}$

Consider $\boldsymbol{\Omega}=(1-\omega)\boldsymbol{I}_{2}+\omega\boldsymbol{1}_{2}\boldsymbol{1}_{2}^{\prime}$. When $s^{2}_{j}\sim\text{gamma}(c,c)$, we have:

[TABLE]

where $|\omega|\leq 1$ by positive semidefiniteness of $\boldsymbol{\Omega}$. It follows that for fixed $c$, the largest possible value of $|\rho|$ is $(\Gamma(c+1/2)/\Gamma(c))^{2}/c$.

I.2 SPB $\boldsymbol{\beta}$

Consider $\boldsymbol{\Omega}=(1-\omega)\boldsymbol{I}_{2}+\omega\boldsymbol{1}_{2}\boldsymbol{1}_{2}^{\prime}$. We have:

[TABLE]

where $|\omega|\leq 1$ by positive semidefiniteness of $\boldsymbol{\Omega}$. It follows that for fixed $q$, the largest possible value of $|\rho|$ is $(\frac{\pi}{2})(\frac{\Gamma(2/q)}{\sqrt{\Gamma(1/q)\Gamma(3/q)}})^{2}$.

Appendix J Posterior Simulation Details

J.1 Univarate Slice Sampling

In several places in this paper, we use a univariate slice sampling algorithm described in Neal (2003) to simulate values of a random variable $x\in[a,b]$ with density proportional to $\text{exp}\{g(x)\}$. Given a previous value $\tilde{x}$, we can simulate a new value $\tilde{x}^{\prime}$ from a full conditional using univariate slice sampling as follows:

Draw $e\sim\text{exp}(1)$. 2. 2.

Draw $d\sim\text{uniform}(a,b)$.

• •

If $g(\tilde{x})-e\leq g(d)$, set $\tilde{x}^{\prime}=d$.

• •

If $g(\tilde{x})-e>g(d)$:

1. (a)

If $d<\tilde{x}$, set $a=d$, else if $d\geq\tilde{x}$ set $b=d$. 2. (b)

Return to 3.

J.2 Simulation from Full Conditional Distribution for $\delta_{j}$

When using the SPB prior, we have $s^{2}_{j}\stackrel{{\scriptstyle d}}{{=}}\Gamma(1/(2\alpha))\xi_{j}^{\frac{1-\alpha}{\alpha}}/(2\Gamma(3/(2\alpha)))$, where $\xi_{j}\sim\text{gamma}((1+\alpha)/(2\alpha),f(\delta_{j}|\alpha))$ and $p(\delta_{j}|\alpha)\propto f(\delta_{j}|\alpha)^{\frac{\alpha-1}{2\alpha}}$ on $(0,\pi)$. Suppose $s^{2}_{j}$ is fixed. Then $\xi_{j}=(2\Gamma(3/(2\alpha))s^{2}_{j}/\Gamma(1/(2\alpha)))^{\frac{\alpha}{1-\alpha}}$. Then we can write the full conditional distribution of $\delta_{j}$ as

[TABLE]

Apply the univariate slice sampling algorithm described in Section J.1 using $g(x)=\text{log}(f(x|\alpha))-f(x|\alpha)\xi_{j}$ and initial values $a=0$ and $b=\pi$.

J.3 Simulation from Full Conditional Distribution for $\rho$

In the real data application, I consider the setting where $\boldsymbol{B}$ and $\boldsymbol{S}$ are $p_{1}\times p_{2}$ matrices and $\boldsymbol{\beta}=\text{vec}(\boldsymbol{B}/\boldsymbol{S})$ has covariance matrix $\mathbb{V}[\boldsymbol{\beta}/\boldsymbol{s}]=\boldsymbol{\Omega}_{2}\otimes\boldsymbol{\Omega}_{1}$, where $\boldsymbol{\Omega}_{1}$ is depends on a single autoregressive parameter, $\rho$, with $\omega_{1,ij}=\rho^{|i-j|}$. The full conditional distribution of $\rho$ is:

[TABLE]

where elements of $\boldsymbol{\Omega}_{1}^{-1}$ are given by

[TABLE]

and $\pi(\rho)$ is the density of an assumed prior distribution for $\rho$. We assume a uniform prior on $(-1,1)$, i.e. $\text{beta}(1,1)$ prior for $(\rho+1)/2$.

J.4 Coordinate Descent for $\boldsymbol{s}$

For any element of $\boldsymbol{s}$, the

[TABLE]

The optimal value of $s$ will satisfy

[TABLE]

When $\kappa_{6}$ is an integer (as is the case for the SNG and SPN priors) we can find all values of $s$ that satisfy this quickly and exactly using fast polynomial root finding functions from the polynom package for R.

When $\kappa_{6}$ is not an integer, find values of $s_{F}$ and $s_{C}$ that maximize

[TABLE]

respectively. This gives us an interval for the solution $s$. This gives us a (hopefully small) interval which contains the maximizing $s$ for the original problem, and we can compute the maximizing $s$ for the original problem using bisection.

[TABLE]

Appendix K Hyperparameter Maximum Marginal Likelihood

The marginal likelihood is given by:

[TABLE]

Given an initial value $\boldsymbol{\Omega}^{(0)}$, Gibbs-within-EM estimates of $\boldsymbol{\Omega}$ can be obtained by iterating the following until $||\boldsymbol{\Omega}^{(i+1)}-\boldsymbol{\Omega}^{(i)}||_{F}$ converges:

• •

using MCMC to simulate $M$ draws $(\boldsymbol{z}^{(1)},s^{(1)}),\dots,(\boldsymbol{z}^{(M)},\boldsymbol{s}^{(M)})$ from the joint posterior distribution of $(\boldsymbol{z},\boldsymbol{s})$ given $\boldsymbol{\Omega}^{(i)}$, $\boldsymbol{\phi}$ and $\boldsymbol{\theta}$, set $\widehat{\mathbb{E}}[\boldsymbol{z}\boldsymbol{z}^{\prime}|\boldsymbol{\Omega}^{(i)},\boldsymbol{X},\boldsymbol{y},\boldsymbol{\phi},\boldsymbol{\theta}]=\frac{1}{M}\sum_{j=1}^{M}\boldsymbol{z}^{(j)}(\boldsymbol{z}^{(j)})^{\prime}$;

• •

set $\boldsymbol{\Omega}^{(i+1)}=\text{argmin}_{\boldsymbol{\Omega}>0}\text{log}(|\boldsymbol{\Omega}|)+\frac{1}{2}\text{tr}(\boldsymbol{\Omega}^{-1}\widehat{\mathbb{E}}[\boldsymbol{z}\boldsymbol{z}^{\prime}|\boldsymbol{\Omega}^{(i)},\boldsymbol{X},\boldsymbol{y},\boldsymbol{\phi},\boldsymbol{\theta}])$.

In practice, we will rarely assume $\boldsymbol{\Omega}$ is unstructured, as doing so requires estimating $p+p(p-1)/2$ unknown parameters from a single observation of a $p\times 1$ vector. Instead, we will parametrize $\boldsymbol{\Omega}$ as a function of lower-dimensional parameters, e.g. as a function of a single autoregressive parameter $\rho$ or as a function of several separable covariance matrices $\boldsymbol{\Omega}=\boldsymbol{\Omega}_{K}\otimes\dots\otimes\boldsymbol{\Omega}_{1}$ which correspond to variation along different modes of the matrix or tensor of regression coefficients $\boldsymbol{B}$. However, the general approach to maximum marginal likelihood estimation of $\boldsymbol{\Omega}$ remains the same, with minimization over $\boldsymbol{\Omega}$ in the second step replaced by minimization over the lower-dimensional components.

Appendix L Seconds per Posterior Draw

Appendix M Minimum Effective Sample Sizes

Appendix N Comparison to Independent Priors

Appendix O Effective Sample Sizes for Application