Spike-and-Slab Group Lassos for Grouped Regression and Sparse Generalized Additive Models

TL;DR

This paper introduces the spike-and-slab group lasso (SSGL), a Bayesian method for variable selection and estimation in grouped linear regression and sparse generalized additive models, with efficient algorithms and uncertainty quantification.

Contribution

It extends the spike-and-slab lasso to nonparametric models, develops a fast MAP estimation algorithm, and provides theoretical guarantees for high-dimensional settings.

Findings

Efficient block coordinate ascent algorithm for MAP estimation.

Theoretical posterior concentration rates in high dimensions.

Successful application to simulations and real data.

Abstract

We introduce the spike-and-slab group lasso (SSGL) for Bayesian estimation and variable selection in linear regression with grouped variables. We further extend the SSGL to sparse generalized additive models (GAMs), thereby introducing the first nonparametric variant of the spike-and-slab lasso methodology. Our model simultaneously performs group selection and estimation, while our fully Bayes treatment of the mixture proportion allows for model complexity control and automatic self-adaptivity to different levels of sparsity. We develop theory to uniquely characterize the global posterior mode under the SSGL and introduce a highly efficient block coordinate ascent algorithm for maximum a posteriori (MAP) estimation. We further employ de-biasing methods to provide uncertainty quantification of our estimates. Thus, implementation of our model avoids the computational intensiveness of Markov chain Monte Carlo (MCMC) in high dimensions. We derive posterior concentration rates for both grouped linear regression and sparse GAMs when the number of covariates grows at nearly exponential rate with sample size. Finally, we illustrate our methodology through extensive simulations and data analysis.

Tables

Table 1. Table 1 : Coverage probabilities for de-biasing simulation.

	$\rho$	Important covariates	Null covariates
$n=100,G=100$	0.0	0.83	0.93
	0.7	0.85	0.94
$n=300,G=300$	0.0	0.93	0.95
	0.7	0.92	0.95

Table 2. Table 2 : Results for SSGL and Group Lasso on the Bardet-Biedl syndrome gene expression data set. In parentheses, we report the standard errors for the CV prediction error.

	SSGL	Group Lasso
# groups selected	12	83
10-fold CV error	0.012 (0.003)	0.017 (0.008)

Table 3. Table 3 : Standardized out-of-sample root mean squared prediction error averaged across 1000 replications for the data sets in Section C.1 . An RMSE of 1 indicates the best performance within a data set.

Data	SSGL	GroupLasso	RandomForest	SoftBart	SuperLearner	BSGS	GroupSpike
Tecator 1	1.41	1.57	2.75	1.93	1.00	5.16	1.67
Tecator 2	1.25	1.58	2.91	1.97	1.00	6.77	1.41
Tecator 3	1.14	1.38	1.94	1.81	1.10	3.31	1.00
BloodBrain	1.10	1.04	1.00	1.01	1.00	1.24	1.13
Wipp	1.44	1.30	1.46	1.00	1.17	4.68	1.30

Table 4. Table 4 : Probes found by SSGL on the Bardet-Biedl syndrome gene expression data set. The probes which were also found by the Group Lasso have nonzero group norm values.

Probe ID	Gene Symbol	SSGL Norm	Group Lasso Norm
1374131_at		0.034
1383749_at	Phospho1	0.067	0.088
1393735_at		0.033	0.002
1379029_at	Zfp62	0.074
1383110_at	Klhl24	0.246
1384470_at	Maneal	0.087	0.005
1395284_at		0.014
1383673_at	Nap1l2	0.045
1379971_at	Zc3h6	0.162
1384860_at	Zfp84	0.008
1376747_at		0.489	0.002
1379094_at		0.220

Table 5. Table 5 : Table displays the terms for which SSGL was found to be enriched in a gene ontology enrichment analysis, ordered by statistical significance.

	SSGL: enriched terms in gene ontology enrichment analysis
1	alpha-mannosidase activity
2	RNA polymerase II intronic transcription regulatory region sequence-specific DNA binding
3	mannosidase activity
4	intronic transcription regulatory region sequence-specific DNA binding
5	intronic transcription regulatory region DNA binding

Equations

$$\mathbf{Y}=\sum _{g=1}^G{\mathbf{X}}_g{\mathbf{\beta }}_g+\mathbf{\epsilon },$$

$$\underset{\mathbf{\beta }}{\mathrm{argmin}}\frac{1}{2}\Vert \mathbf{Y}-\sum _{g=1}^G{\mathbf{X}}_g{\mathbf{\beta }}_g{\Vert }_2^2+\lambda \sum _{g=1}^G\sqrt{m_g}\Vert {\mathbf{\beta }}_g{\Vert }_2,$$

$$\displaystyle\begin{array}[]{l}\pi(\bm{\beta}|\bm{\gamma})=\displaystyle\prod_{g=1}^{G}[(1-\gamma_{g})\delta_{0}(\bm{\beta}_{g})+\gamma_{g}\pi(\bm{\beta}_{g})],\\ \pi(\bm{\gamma}|\theta)=\displaystyle\prod_{g=1}^{G}\theta^{\gamma_{g}}(1-\theta)^{1-\gamma_{g}},\\ \theta\sim\pi(\theta),\end{array}$$

$$\pi (\mathbf{\beta }|\theta )=\prod _{j=1}^p[(1-\theta )\psi ({\beta }_j|{\lambda }_0)+\theta \psi ({\beta }_j|{\lambda }_1)],$$

$$\mathbf{\Psi }({\mathbf{\beta }}_g|\lambda )=C_g{\lambda }^{m_g}\exp \left(-\lambda \Vert {\mathbf{\beta }}_g{\Vert }_2\right),$$

$$\pi (\mathbf{\beta }|\theta )=\prod _{g=1}^G\left[(1-\theta )\mathbf{\Psi }({\mathbf{\beta }}_g|{\lambda }_0)+\theta \mathbf{\Psi }({\mathbf{\beta }}_g|{\lambda }_1)\right],$$

$$\pi ({\sigma }^2)\propto {\sigma }^{-2}.$$

$$pe{n}_S(\mathbf{\beta }|\theta )$$

$$p_{\theta }^{\ast }({\mathbf{\beta }}_g)=\frac{\theta {\mathbf{\Psi }}_1({\mathbf{\beta }}_g)}{\theta {\mathbf{\Psi }}_1({\mathbf{\beta }}_g)+(1-\theta ){\mathbf{\Psi }}_0({\mathbf{\beta }}_g)}.$$

$$\frac{\partial pe{n}_S(\mathbf{\beta }|\theta )}{\partial \Vert {\mathbf{\beta }}_g{\Vert }_2}=-{\lambda }_{\theta }^{\ast }({\mathbf{\beta }}_g)$$

$${\lambda }_{\theta }^{\ast }({\mathbf{\beta }}_g)={\lambda }_1{p}_{\theta }^{\ast }({\mathbf{\beta }}_g)+{\lambda }_0[1-p_{\theta }^{\ast }({\mathbf{\beta }}_g)].$$

$${\mathbf{X}}_g^T(\mathbf{Y}-\mathbf{X}\hat{\mathbf{\beta }})={\sigma }^2{\lambda }_{\theta }^{\ast }({\hat{\mathbf{\beta }}}_g)\frac{{\hat{\mathbf{\beta }}}_g}{\Vert {\mathbf{\beta }}_g{\Vert }_2}$$

$$\Vert {\mathbf{X}}_g^T(\mathbf{Y}-\mathbf{X}\hat{\mathbf{\beta }}){\Vert }_2\le {\sigma }^2{\lambda }_{\theta }^{\ast }({\hat{\mathbf{\beta }}}_g)$$

$${\hat{\mathbf{\beta }}}_g=\frac{1}{n}{\left(1-\frac{{\sigma }^2{\lambda }_{\theta }^{\ast }({\hat{\mathbf{\beta }}}_g)}{\Vert {\mathbf{z}}_g{\Vert }_2}\right)}_{+}{\mathbf{z}}_g$$

$$\Delta =\underset{{\mathbf{\beta }}_g}{\inf }\left\{\frac{n\Vert {\mathbf{\beta }}_g{\Vert }_2}{2}-\frac{{\sigma }^{2}pe{n}_S(\mathbf{\beta }|\theta )}{\Vert {\mathbf{\beta }}_g{\Vert }_2}\right\}.$$

$$h({\mathbf{\beta }}_g)=[{\lambda }_{\theta }^{\ast }({\mathbf{\beta }}_g)-{\lambda }_1{]}^2+\frac{2n}{{\sigma }^2}\log p_{\theta }^{\ast }({\mathbf{\beta }}_g).$$

$${\Delta }^L<\Delta <{\Delta }^U$$

$${\Delta }^L$$

$${\Delta }^U$$

$$0<d<\frac{2n}{{\sigma }^2}-{\left(\frac{n}{{\sigma }^2({\lambda }_0-{\lambda }_1)}-\frac{\sqrt{2n}}{\sigma }\right)}^2$$

$$\pi (\mathbf{\beta })$$

$$=\left(\prod _{g=1}^G{C}_g{\lambda }_1^{m_g}\right)e^{-{\lambda }_1{\sum }_{g=1}^G\Vert {\mathbf{\beta }}_g{\Vert }_2}{\int }_0^1\frac{{\theta }^G}{{\prod }_{g=1}^G{p}_{\theta }^{\ast }({\mathbf{\beta }}_g)}d\pi (\theta ),$$

$$pe{n}_{NS}(\mathbf{\beta })$$

$$\frac{\partial pe{n}_{NS}(\mathbf{\beta })}{\partial \Vert {\mathbf{\beta }}_g{\Vert }_2}\equiv {\lambda }^{\ast }({\mathbf{\beta }}_g;{\mathbf{\beta }}_{\backslash g}),$$

$${\lambda }^{\ast }({\mathbf{\beta }}_g;{\mathbf{\beta }}_{\backslash g})=p^{\ast }({\mathbf{\beta }}_g;{\mathbf{\beta }}_{\backslash g}){\lambda }_1+[1-p^{\ast }({\mathbf{\beta }}_g;{\mathbf{\beta }}_{\backslash g})]{\lambda }_0$$

$$p^{\ast }({\mathbf{\beta }}_g;{\mathbf{\beta }}_{\backslash g})\equiv p_{{\theta }_g}^{\ast }({\mathbf{\beta }}_g),\text{with}{\theta }_g=\mathbb{E}[\theta |{\mathbf{\beta }}_{\backslash g}].$$

$$\pi (\theta |\hat{\mathbf{\beta }})\propto {\theta }^{a-1}(1-\theta {)}^{b-1}(1-\theta z{)}^{G-\hat{q}}\prod _{g=1}^{\hat{q}}(1-\theta x_g),$$

$$\mathbb{E}[\theta |\hat{\mathbf{\beta }}]=\frac{{\int }_0^1{\theta }^a(1-\theta {)}^{b-1}(1-\theta z{)}^{G-\hat{q}}{\prod }_{g=1}^{\hat{q}}(1-\theta x_g)d\theta }{{\int }_0^1{\theta }^{a-1}(1-\theta {)}^{b-1}(1-\theta z{)}^{G-\hat{q}}{\prod }_{g=1}^{\hat{q}}(1-\theta x_g)d\theta }.$$

$$\mathbb{E}[\theta |\hat{\mathbf{\beta }}]=\frac{a+\hat{q}}{a+b+G}.$$

$${\hat{\mathbf{\beta }}}_g\leftarrow \frac{1}{n}{\left(1-\frac{{\sigma }^2{\lambda }_{\hat{\theta }}^{\ast }({\hat{\mathbf{\beta }}}_g)}{\Vert {\mathbf{z}}_g{\Vert }_2}\right)}_{+}{\mathbf{z}}_g,$$

$${\mathbf{\beta }}_g^{(k)}\leftarrow \frac{1}{n}{\left(1-\frac{{\sigma }^{2(k)}{\lambda }^{\ast }({\mathbf{\beta }}_g^{(k-1)};{\theta }^{(k)})}{\Vert {\mathbf{z}}_g{\Vert }_2}\right)}_{+}{\mathbf{z}}_g\mathbb{I}(\Vert {\mathbf{z}}_g{\Vert }_2>{\Delta }^U)$$

$${\sigma }^{2(k)}\leftarrow \frac{\Vert \mathbf{Y}-\mathbf{X}{\mathbf{\beta }}^{(k-1)}{\Vert }_2^2}{n+2}.$$

Full text

1 Introduction

1.1 Regression with Grouped Variables

Group structure arises in many statistical applications. For example, in multifactor analysis of variance, multi-level categorical predictors are each represented by a group of dummy variables. In genomics, genes within the same pathway may form a group at the pathway or gene set level and act in tandem to regulate a biological system. In each of these scenarios, the response $\bm{Y}_{n\times 1}$ can be modeled as a linear regression problem with $G$ groups:

[TABLE]

where $\bm{\varepsilon}\sim\mathcal{N}_{n}(\bm{0},\sigma^{2}\bm{I}_{n})$, $\bm{\beta}_{g}$ is a coefficients vector of length $m_{g}$, and $\bm{X}_{g}$ is an $n\times m_{g}$ covariate matrix corresponding to group $g=1,\ldots G$. Even in the absence of grouping information about the covariates, the model (1.1) subsumes a wide class of important nonparametric regression models called generalized additive models (GAMs). In GAMs, continuous covariates may be represented by groups of basis functions which have a nonlinear relationship with the response. We defer further discussion of GAMs to Section 5.

It is often of practical interest to select groups of variables that are most significantly associated with the response. To facilitate this group-level selection, Yuan and Lin [54] introduced the group lasso, which solves the optimization problem,

[TABLE]

where $||\cdot||_{2}$ is the $\ell_{2}$ norm. In the frequentist literature, many variants of model (1.2) have been introduced, which use some combination of $\ell_{1}$ and $\ell_{2}$ penalties on the coefficients of interest (e.g., [15, 19, 39]).

In the Bayesian framework, selection of relevant groups under model (1.1) is often done by placing spike-and-slab priors on each of the groups $\bm{\beta}_{g}$ (e.g., [49, 23, 50, 29]). These priors typically take the form,

[TABLE]

where $\bm{\gamma}$ is a binary vector that indexes the $2^{G}$ possible models, $\theta\in(0,1)$ is the mixing proportion, $\delta_{0}$ is a point mass at $\bm{0}_{m_{g}}\in\mathbb{R}^{m_{g}}$ (the “spike”), and $\pi(\bm{\beta}_{g})$ is an appropriate “slab” density (typically a multivariate normal distribution or a scale-mixture multivariate normal density). With a well-chosen prior on $\theta$, this model will favor parsimonious models in very high dimensions, thus avoiding the curse of dimensionality.

1.2 The Spike-and-Slab Lasso

For Bayesian variable selection, point mass spike-and-slab priors (1.6) are interpretable, but they are computationally intractable in high dimensions, due in large part to the combinatorial complexity of updating the discrete indicators $\bm{\gamma}$. As an alternative, fully continuous variants of spike-and-slab models have been developed. For continuous spike-and-slab models, the point mass spike $\delta_{0}$ is replaced by a continuous density heavily concentrated around $\bm{0}_{m_{g}}$. This not only mimics the point mass but it also facilitates more efficient computation, as we describe later.

In the context of sparse normal means estimation and univariate linear regression, Ročková [34] and Ročková and George [36] introduced the univariate spike-and-slab lasso (SSL). The SSL places a mixture prior of two Laplace densities on the individual coordinates $\beta_{j}$, i.e.

[TABLE]

where $\theta\in(0,1)$ is the mixing proportion and $\psi(\cdot|\lambda)$ denotes a univariate Laplace density indexed by hyperparameter $\lambda$, i.e. $\psi(\beta|\lambda)=\frac{\lambda}{2}e^{-\lambda|\beta|}$. Typically, we set $\lambda_{0}\gg\lambda_{1}$ so that the spike is heavily concentrated about zero. Unlike (1.6), the SSL model (1.7) does not place any mass on exactly sparse vectors. Nevertheless, the global posterior mode under the SSL prior may be exactly sparse. Meanwhile, the slab stabilizes posterior estimates of the larger coefficients so they are not downward biased. Thus, the SSL posterior mode can be used to perform variable selection and estimation simultaneously.

The spike-and-slab lasso methodology has now been adopted for a wide number of statistical problems. Apart from univariate linear regression, it has been used for factor analysis [35, 27], multivariate regression [8], covariance/precision matrix estimation [8, 10, 20], causal inference [2], generalized linear models (GLMs) [44, 42], and Cox proportional hazards models [43].

While the SSL (1.7) induces sparsity on individual coefficients (through the posterior mode), it does not account for group structure of covariates. For inference with structured data in GLMs, Tang et al. [42] utilized the univariate spike-and-slab lasso prior (1.7) for grouped data where each group had a group-specific sparsity-inducing parameter, $\theta_{g}$, instead of a single $\theta$ for all coefficients. However, this univariate SSL prior does not feature the “all in, all out” selection property of the original group lasso of Yuan and Lin [54] or the grouped and multivariate SSL prior, which we develop in this work.

In this paper, we introduce the spike-and-slab group lasso (SSGL) for Bayesian grouped regression and variable selection. Under the SSGL prior, the global posterior mode is exactly sparse, thereby allowing the mode to automatically threshold out insignificant groups of coefficients. To widen the use of spike-and-slab lasso methodology for situations where the linear model is too inflexible, we extend the SSGL to sparse generalized additive models by introducing the nonparametric spike-and-slab lasso (NPSSL). To our knowledge, our work is the first to apply the spike-and-slab lasso methodology outside of a parametric setting. Our contributions can be summarized as follows:

We propose a new group spike-and-slab prior for estimation and variable selection in both parametric and nonparametric settings. Unlike frequentist methods which rely on separable penalties, our model has a non-separable and self-adaptive penalty which allows us to automatically adapt to ensemble information about sparsity. 2. 2.

We introduce a highly efficient block coordinate ascent algorithm for global posterior mode estimation. This allows us to rapidly identify significant groups of coefficients, while thresholding out insignificant ones. 3. 3.

We show that de-biasing techniques that have been used for the original lasso [45] can be extended to our SSGL model to provide valid inference on the estimated regression coefficients. 4. 4.

For both grouped regression and sparse additive models, we derive near-optimal posterior contraction rates for both the regression coefficients $\bm{\beta}$ and the unknown variance $\sigma^{2}$ under the SSGL prior.

The rest of the paper is structured as follows. In Section 2, we introduce the spike-and-slab group lasso (SSGL). In Section 3, we characterize the global posterior mode and introduce efficient algorithms for fast maximum a posteriori (MAP) estimation and variable selection. In Section 4, we utilize ideas from the de-biased lasso to perform inference on the SSGL model. In Section 5, we extend the SSGL to nonparametric settings by proposing the nonparametric spike-and-slab lasso (NPSSL). In Section 6, we present asymptotic theory for the SSGL and the NPSSL. Finally, in Sections 7 and 8, we provide extensive simulation studies and use our models to analyze real data sets.

1.3 Notation

We use the following notations. For two nonnegative sequences $\{a_{n}\}$ and $\{b_{n}\}$, we write $a_{n}\asymp b_{n}$ to denote $0<\lim\inf_{n\rightarrow\infty}a_{n}/b_{n}\leq\lim\sup_{n\rightarrow\infty}a_{n}/b_{n}<\infty$. If $\lim_{n\rightarrow\infty}a_{n}/b_{n}=0$, we write $a_{n}=o(b_{n})$ or $a_{n}\prec b_{n}$. We use $a_{n}\lesssim b_{n}$ or $a_{n}=O(b_{n})$ to denote that for sufficiently large $n$, there exists a constant $C>0$ independent of $n$ such that $a_{n}\leq Cb_{n}$. For a vector $\bm{v}\in\mathbb{R}^{p}$, we let $\lVert\bm{v}\rVert_{1}:=\sum_{i=1}^{p}|v_{i}|$, $\lVert\bm{v}\rVert_{2}:=\sqrt{\sum_{i=1}^{p}v_{i}^{2}}$, and $\lVert\bm{v}\rVert_{\infty}:=\max_{1\leq i\leq p}|v_{i}|$ denote the $\ell_{1}$, $\ell_{2}$, and $\ell_{\infty}$ norms respectively. For a symmetric matrix $\bm{A}$, we let $\lambda_{\min}(\bm{A})$ and $\lambda_{\max}(\bm{A})$ denote its minimum and maximum eigenvalues.

2 The Spike-and-Slab Group Lasso

Let $\bm{\beta}_{g}$ denote a real-valued vector of length $m_{g}$. We define the group lasso density as

[TABLE]

where $C_{g}=2^{-m_{g}}\pi^{-(m_{g}-1)/2}\left[\Gamma\left((m_{g}+1)/2\right)\right]^{-1}$. This prior has been previously considered by [18, 49] for Bayesian inference in the grouped regression model (1.1). Kyung et al. [18] considered a single prior (2.1) on each of the $\bm{\beta}_{g}$’s, while Xu and Ghosh [49] employed (2.1) as the slab in the point-mass mixture (1.6). These authors implemented their models using MCMC.

In this manuscript, we introduce a continuous spike-and-slab prior with the group lasso density (2.1) for both the spike and the slab. The continuous nature of our prior is critical in facilitating efficient coordinate ascent algorithms for MAP estimation that allow us to bypass the use of MCMC. Letting $\bm{\beta}=(\bm{\beta}_{1}^{T},\ldots,\bm{\beta}_{G}^{T})^{T}$ under model (1.1), the spike-and-slab group lasso (SSGL) is defined as:

[TABLE]

where $\bm{\Psi}(\cdot|\lambda)$ denotes the group lasso density (2.1) indexed by hyperparameter $\lambda$, and $\theta\in(0,1)$ is a mixing proportion. $\lambda_{0}$ corresponds to the spike which shrinks the entire vector $\bm{\beta}_{g}$ towards $\bm{0}_{m_{g}}$, while $\lambda_{1}$ corresponds to the slab. For shorthand notation, we denote $\bm{\Psi}(\bm{\beta}_{g}|\lambda_{0})$ as $\bm{\Psi}_{0}(\bm{\beta}_{g})$ and $\bm{\Psi}(\bm{\beta}_{g}|\lambda_{1})$ as $\bm{\Psi}_{1}(\bm{\beta}_{g})$ going forward.

Under the grouped regression model (1.1), we place the SSGL prior (2.2) on $\bm{\beta}$. In accordance with the recommendations of [28], we do not scale our prior by the unknown $\sigma$. Instead, we place an independent Jeffreys prior on $\sigma^{2}$, i.e.

[TABLE]

The mixing proportion $\theta$ in (2.2) can either be fixed deterministically or endowed with a prior $\theta\sim\pi(\theta)$. We will discuss this in detail in Section 3.

3 Characterization and Computation of the Global Posterior Mode

Throughout this section, we let $p$ denote the total number of covariates, i.e. $p=\sum_{g=1}^{G}m_{g}$. Our goal is to find the maximum a posteriori estimates of the regression coefficients $\bm{\beta}\in\mathbb{R}^{p}$. This optimization problem is equivalent to a penalized likelihood method in which the logarithm of the prior (2.2) may be reinterpreted as a penalty on the regression coefficients. Similarly to Ročková and George [36], we will leverage this connection between the Bayesian and frequentist paradigms and introduce the SSGL penalty. This strategy combines the adaptivity of the Bayesian approach with the computational efficiency of existing algorithms in the frequentist literature.

A key component of the SSGL model is $\theta$, the prior expected proportion of groups with large coefficients. Ultimately, we will pursue a fully Bayes approach and place a prior on $\theta$, allowing the SSGL to adapt to the underlying sparsity of the data and perform an automatic multiplicity adjustment [38]. For ease of exposition, however, we will first consider the case where $\theta$ is fixed, echoing the development of Ročková and George [36]. In this situation, the regression coefficients $\bm{\beta}_{g}$ are conditionally independent a priori, resulting in a separable SSGL penalty. Later we will consider the fully Bayes approach, which will yield the non-separable SSGL penalty.

Definition 1**.**

Given $\theta\in(0,1)$, the separable SSGL penalty is defined as

[TABLE]

where

[TABLE]

The separable SSGL penalty is almost the logarithm of the original prior (2.2); the only modification is an additive constant to ensure that $pen_{S}(\mathbf{0}_{p}|\theta)=0$. The connection between the SSGL and penalized likelihood methods is made clearer when considering the derivative of the separable SSGL penalty, given in the following lemma.

Lemma 1**.**

The derivative of the separable SSGL penalty satisfies

[TABLE]

where

[TABLE]

Similarly to the SSL, the SSGL penalty is a weighted average of the two regularization parameters, $\lambda_{1}$ and $\lambda_{0}$. The weight $p^{*}_{\theta}(\bm{\beta}_{g})$ is the conditional probability that $\bm{\beta}_{g}$ was drawn from the slab distribution rather than the spike. Hence, the SSGL features an adaptive regularization parameter which applies different amounts of shrinkage to each group, unlike the group lasso which applies the same shrinkage to each group.

3.1 The Global Posterior Mode

Similarly to the group lasso [54], the separable nature of the penalty (3.1) lends itself naturally to a block coordinate ascent algorithm which cycles through the groups. In this section, we first outline the group updates resulting from the Karush-Kuhn-Tucker (KKT) conditions. The KKT conditions provide necessary conditions for the global posterior mode. We then derive a more refined condition for the global mode to aid in optimization for multimodal posteriors.

Following Huang et al. [14], we assume that within each group, covariates are orthonormal, i.e. $\bm{X}_{g}^{T}\bm{X}_{g}=n\mathbf{I}_{m_{g}}$ for $g=1,\dots,G$. If this assumption does not hold, then the $\bm{X}_{g}$ matrices can be orthonormalized before fitting the model. As noted by Breheny and Huang [4], orthonormalization can be done without loss of generality since the resulting solution can be transformed back to the original scale.

Proposition 1**.**

The necessary conditions for $\widehat{\bm{\beta}}=(\widehat{\bm{\beta}}_{1}^{T},\dots,\widehat{\bm{\beta}}_{G}^{T})^{T}$ to be a global mode are:

[TABLE]

Equivalently,

[TABLE]

where $\bm{z}_{g}=\bm{X}_{g}^{T}\left[\bm{Y}-\sum_{l\neq g}\bm{X}_{l}\widehat{\bm{\beta}}_{l}\right].$

Proof.

Follows immediately from Lemma 1 and subdifferential Calculus. ∎

The above characterization for the global mode is necessary, but not sufficient. A more refined characterization may be obtained by considering the group-wise optimization problem, noting that the global mode is also a maximizer of the $g$th group, keeping all other groups fixed.

Proposition 2**.**

The global mode $\widehat{\bm{\beta}}_{g}=\mathbf{0}_{m_{g}}$ if and only if $\lVert\bm{z}_{g}\rVert_{2}\leq\Delta$, where

[TABLE]

The proof for Proposition 2 can be found in Appendix D.2. Unfortunately, the threshold $\Delta$ is difficult to compute. We instead find an approximation to this threshold. An upper bound is simply that of the soft-threshold solution (3.7), with $\Delta\leq\sigma^{2}\lambda^{*}(\bm{\beta}_{g})$. However, when $\lambda_{0}$ is large, this bound may be improved. Similarly to Ročková and George [36], we provide improved bounds on the threshold in Theorem 1. This result requires the function $h:\mathbb{R}^{m_{g}}\to\mathbb{R}$, defined as:

[TABLE]

Theorem 1**.**

When $(\lambda_{0}-\lambda_{1})>2\sqrt{n}/\sigma$ and $h(\mathbf{0}_{m_{g}})>0$, the threshold $\Delta$ is bounded by:

[TABLE]

where

[TABLE]

and

[TABLE]

When $\lambda_{0}$ is large, $d\to 0$ and the lower bound on the threshold approaches the upper bound, yielding the approximation $\Delta=\Delta^{U}$. We will ultimately use this approximation in our block coordinate ascent algorithm.

3.2 The Non-Separable SSGL penalty

As discussed earlier, a key reason for adopting a Bayesian strategy is that it allows the model to borrow information across groups and self-adapt to the true underlying sparsity in the data. This is achieved by placing a prior on $\theta$, the proportion of groups with non-zero coefficients. We now outline this fully Bayes strategy and the resulting non-separable SSGL penalty. With the inclusion of the prior $\theta\sim\pi(\theta)$, the marginal prior for the regression coefficients has the following form:

[TABLE]

The non-separable SSGL penalty is then defined similarly to the separable penalty, where again we have centered the penalty to ensure $pen_{NS}(\mathbf{0}_{p})=0$.

Definition 2**.**

The non-separable SSGL (NS-SSGL) penalty with $\theta\sim\pi(\theta)$ is defined as

[TABLE]

Although the penalty (3.14) appears intractable, intuition is again obtained by considering the derivative. Following the same line of argument as Ročková and George [36], the derivative of (3.14) is given in the following lemma.

Lemma 2**.**

[TABLE]

where

[TABLE]

and

[TABLE]

That is, the marginal prior from (3.14) is rendered tractable by considering each group of regression coefficients separately, conditional on the remaining coefficients. Such a conditional strategy is motivated by the group-wise updates for the separable penalty considered in the previous section. Thus, our optimization strategy for the non-separable penalty will be very similar to the separable case, except instead of a fixed value for $\theta$, we will impute the mean of $\theta$ conditioned on the remaining regression coefficients.

We now consider the form of the conditional mean, $\mathbb{E}[\theta|\widehat{\bm{\beta}}_{\backslash g}]$. As noted by Ročková and George [36], when the number of groups is large, this conditional mean can be replaced by $\mathbb{E}[\theta|\widehat{\bm{\beta}}]$; we will proceed with the same approximation. For the prior on $\theta$, we will use the standard beta prior $\theta\sim\mathcal{B}(a,b)$. With the choices $a=1$ and $b=G$ for these hyperparameters, this prior results in an automatic multiplicity adjustment for the regression coefficients [38].

We now examine the conditional distribution $\pi(\theta|\widehat{\bm{\beta}})$. Suppose that the number of groups with non-zero coefficients is $\widehat{q}$, and assume without loss of generality that the first $\widehat{q}$ groups have non-zero coefficients. Then,

[TABLE]

with $z=1-\frac{\lambda_{1}}{\lambda_{0}}$ and $x_{g}=(1-\frac{\lambda_{1}}{\lambda_{0}}e^{\lVert\widehat{\bm{\beta}}_{g}\rVert_{2}(\lambda_{0}-\lambda_{1})})$. Similarly to Ročková and George [36], this distribution is a generalization of the Gauss hypergeometric distribution. Consequently, the expectation may be written as

[TABLE]

While the above expression (3.20) appears laborious to compute, it admits a much simpler form when $\lambda_{0}$ is very large. Using a slight modification to the arguments of [33], we obtain this simpler form in Lemma 3.

Lemma 3**.**

Assume $\pi(\theta|\widehat{\bm{\beta}})$ is distributed according to (3.19). Let $\widehat{q}$ be the number of groups with non-zero coefficients. Then as $\lambda_{0}\to\infty$,

[TABLE]

The proof for Lemma 3 is in Appendix D.2. We note that the expression (3.21) is essentially the usual posterior mean of $\theta$ under a beta prior. Intuitively, as $\lambda_{0}$ diverges, the weights $p_{\theta}^{*}(\bm{\beta}_{g})$ concentrate at zero and one, yielding the familiar form for $\mathbb{E}[\theta|\widehat{\bm{\beta}}]$. With this in hand, we are now in a position to outline the block coordinate ascent algorithm for the non-separable SSGL.

3.3 Optimization

The KKT conditions for the non-separable SSGL penalty yield the following necessary condition for the global mode:

[TABLE]

where $\bm{z}_{g}=\bm{X}_{g}^{T}\left[\bm{Y}-\sum_{l\neq g}\bm{X}_{l}\widehat{\bm{\beta}}_{l}\right]$ and $\widehat{\theta}$ is the mean (3.21), conditioned on the previous value of $\bm{\beta}$. As before, (3.22) is sufficient for a local mode, but not the global mode. When $p\gg n$ and $\lambda_{0}$ is large, the posterior will be highly multimodal. As in the separable case, we require a refined thresholding scheme that will eliminate some of these suboptimal local modes from consideration. In approximating the group-wise conditional mean $\mathbb{E}[\theta|\widehat{\bm{\beta}}_{\backslash g}]$ with $\mathbb{E}[\theta|\widehat{\bm{\beta}}]$, we do not require group-specific thresholds. Instead, we can use the threshold given in Proposition 2 and Theorem 1 where $\theta$ is replaced with the current update (3.21). In particular, we shall use the upper bound $\Delta^{U}$ in our block coordinate ascent algorithm.

Similarly to Ročková and George [36], we combine the refined threshold, $\Delta^{U}$ with the soft thresholding operation (3.22), to yield the following update for $\widehat{\bm{\beta}}_{g}$ at iteration $k$:

[TABLE]

where $\theta^{(k)}=\mathbb{E}[\theta|\bm{\beta}^{(k-1)}]$. Technically, $\theta$ should be updated after each group $\bm{\beta}_{g}$ is updated. In practice, however, there will be little change after one group is updated and so we will update both $\theta$ and $\Delta^{U}$ after every $M$ iterations with a default value of $M=10$.

With the Jeffreys prior $\pi(\sigma^{2})\propto\sigma^{-2}$, the error variance $\sigma^{2}$ also has a closed form update:

[TABLE]

The complete optimization algorithm is given in Algorithm 1 of Appendix A.1. The computational complexity of this algorithm is $\mathcal{O}(np)$ per iteration, where $p=\sum_{g=1}^{G}m_{g}$. It takes $\mathcal{O}(nm_{g})$ operations to compute the partial residual $\bm{z}_{g}$ for the $g$th group, for a total cost of $\mathcal{O}(n\sum_{g=1}^{G}m_{g})=\mathcal{O}(np)$. Similarly, it takes $\mathcal{O}(np)$ cost to compute the sum of squared residuals $\lVert\bm{Y}-\bm{X}\widehat{\bm{\beta}}\rVert_{2}^{2}$ to update the variance parameter $\sigma^{2}$. The computational complexity of our algorithm matches that of the usual gradient descent algorithms for lasso and group lasso [9].

As a non-convex method, it is not guaranteed that SSGL will find the global posterior mode, only a local mode. However, the refined thresholding scheme (Theorem 1) and a warm start initialization strategy (described in detail in Appendix A.2) enable SSGL to eliminate a number sub-optimal local modes from consideration in a similar manner to Ročková and George [36]. To briefly summarize the initialization strategy, we tune $\lambda_{0}$ from an increasing sequence of values, and we further scale $\lambda_{0}$ by $\sqrt{m_{g}}$ for each $g$th group to ensure that the amount of penalization is on the same scale for groups of potentially different sizes [14]. Meanwhile, we keep $\lambda_{1}$ fixed at a small value so that selected groups have minimal shrinkage. See Appendix A.2 for detailed discussion of choosing $(\lambda_{0},\lambda_{1})$.

4 Approaches to Inference

While the above procedure allows us to find the posterior mode of $\bm{\beta}$, providing a measure of uncertainty around our estimate is a challenging task. One possible solution is to run MCMC where the algorithm is initialized at the posterior mode. By starting the MCMC chain at the mode, the algorithm should converge faster. However, this is still not ideal, as it can be computationally burdensome in high dimensions. Instead, we will adopt ideas from a recent line of research ([46, 16]) based on de-biasing estimates from high-dimensional regression. These ideas were derived in the context of lasso regression, and we will explore the extent to which they work for the SSGL penalty. Define $\widehat{\bm{\Sigma}}=\bm{X}^{T}\bm{X}/n$ and let $\widehat{\bm{\Theta}}$ be an approximate inverse of $\widehat{\bm{\Sigma}}$. We define

[TABLE]

where $\widehat{\bm{\beta}}$ is the MAP estimator of $\bm{\beta}$ under the SSGL model. By [46], this quantity $\widehat{\bm{\beta}}_{d}$ has the following asymptotic distribution:

[TABLE]

For our inference procedure, we replace the population variance $\sigma^{2}$ in (4.2) with the modal estimate $\widehat{\sigma}^{2}$ from the SSGL model. To estimate $\widehat{\bm{\Theta}}$, we utilize the nodewise regression approach developed in [24, 46]. We describe this estimation procedure for $\widehat{\bm{\Theta}}$ in Appendix A.3.

Let $\widehat{\beta}_{dj}$ denote the $j$th coordinate of $\widehat{\bm{\beta}}_{d}$. We have from (4.2) that the $100(1-\alpha)\%$ asymptotic pointwise confidence intervals for $\beta_{j},j=1,\ldots,p$, are

[TABLE]

where $c(\alpha,n,\widehat{\sigma}^{2}):=\Phi^{-1}(1-\alpha/2)\sqrt{\widehat{\sigma}^{2}(\widehat{\bm{\Theta}}\widehat{\bm{\Sigma}}\widehat{\bm{\Theta}}^{T})_{jj}/n}$ and $\Phi(\cdot)$ denotes the cdf of $\mathcal{N}(0,1)$. It should be noted that our posterior mode estimates should have less bias than existing estimates such as the group lasso. Therefore, the goal of the de-biasing procedure is less about de-biasing the posterior mode estimates, and more about providing an estimator with an asymptotic normal distribution from which we can perform inference.

To assess the ability of this procedure to obtain accurate confidence intervals (4.3) with $\alpha=0.05$, we run a small simulation study with $n=100$, $G=100$ or $n=300,G=300$, and each of the $G$ groups having $m=2$ covariates. We generate the covariates from a multivariate normal distribution with mean $\bm{0}$ and an AR(1) covariance structure with correlation $\rho$. The two covariates from each group are the linear and squared term from the original covariates. We set the first seven elements of $\bm{\beta}$ equal to $(0,0.5,0.25,0.1,0,0,0.7)$ and the remaining elements equal to zero. Lastly, we try $\rho=0$ and $\rho=0.7$. Table 1 shows the coverage probabilities across 1000 simulations for all scenarios looked at. We see that important covariates, i.e. covariates with a nonzero corresponding $\beta_{j}$, have coverage near 0.85 when $n=100$ under either correlation structure, though this increases to nearly the nominal rate when $n=300$. The remaining covariates (null covariates) achieve the nominal level regardless of the sample size or correlation present.

5 Nonparametric Spike-and-Slab Lasso

We now introduce the nonparametric spike-and-slab lasso (NPSSL). The NPSSL allows for flexible modeling of a response surface with minimal assumptions regarding its functional form. We consider two cases for the NPSSL: (i) a main effects only model, and (ii) a model with both main and interaction effects.

5.1 Main Effects

We first consider the main effects NPSSL model. Here, we assume that the response surface may be decomposed into the sum of univariate functions of each of the $p$ covariates. That is, we have the following model:

[TABLE]

Following Ravikumar et al. [32], we assume that each $f_{j}$, $j=1,\dots,p$, may be approximated by a linear combination of basis functions $\mathcal{B}_{j}=\{g_{j1},\dots,g_{jd}\}$, i.e.,

[TABLE]

where $\bm{\beta}_{j}=(\beta_{j1},\dots,\beta_{jd})^{T}$ are the unknown weights. Let $\widetilde{\bm{X}}_{j}$ denote the $n\times d$ matrix with the $(i,k)$th entry $\widetilde{\bm{X}}_{j}(i,k)=g_{jk}(X_{ij})$. Then, (5.1) may be represented in matrix form as

[TABLE]

where $\bm{\delta}$ is a vector of the lower-order truncation bias. Note that we assume the response $\bm{Y}$ has been centered and so we do not include a grand mean $\bm{\mu}$ in (5.3). Thus, we do not require the main effects to integrate to zero as in [48]. We do, however, require the matrices $\widetilde{\bm{X}}_{j},j=1,\ldots,p$, to be orthogonal, as discussed in Section 3. Note that the entire design matrix does not need to be orthogonal; only the group-specific matrices need to be. We can enforce this in practice by either using orthonormal basis functions or by orthornormalizing the $\widetilde{\bm{X}}_{j}$ matrices before fitting the model.

We assume that $\bm{Y}$ depends on only a small number of the $p$ covariates so that many of the $f_{j}$’s have a negligible contribution to (5.1). This is equivalent to assuming that most of the weight vectors $\bm{\beta}_{j}$ have all zero elements. If the $j$th covariate is determined to be predictive of $\bm{Y}$, then $f_{j}$ has a non-negligible contribution to (5.1). In this case, we want to include the entire basis function approximation to $f_{j}$ in the model.

The above situation is a natural fit for the SSGL. We have $p$ groups where each group is either included as a whole or not included in the model. The design matrices for each group are exactly the matrices of basis functions, $\widetilde{\bm{X}}_{j},j=1,\ldots,p$. We will utilize the non-separable SSGL penalty developed in Section 3.2 to enforce this group-sparsity behavior in the model (5.3). More specifically, we seek to maximize the objective function with respect to $\bm{\beta}=(\bm{\beta}_{1}^{T},\dots,\bm{\beta}_{p}^{T})^{T}\in\mathbb{R}^{pd}$ and $\sigma^{2}$:

[TABLE]

To find the estimators of $\bm{\beta}$ and $\sigma^{2}$, we use Algorithm 1 in Appendix A.1. Similar additive models have been proposed by a number of authors including Ravikumar et al. [32] and Wei et al. [48]. However, our proposed NPSSL method has a number of advantages. First, we allow the noise variance $\sigma^{2}$ to be unknown, unlike Ravikumar et al. [32]. Accurate estimates of $\sigma^{2}$ are important to avoid overfitting the noise beyond the signal. Secondly, we use a block-descent algorithm to quickly target the modes of the posterior, whereas Wei et al. [48] utilize MCMC. Finally, our SSGL algorithm automatically thresholds negligible groups to zero, negating the need for a post-processing thresholding step.

5.2 Main and Interaction Effects

The main effects model (5.1) allows for each covariate to have a nonlinear contribution to the model, but assumes a linear relationship between the covariates. In some applications, this assumption may be too restrictive. For example, in the environmental exposures data which we analyze in Section 8.2, we may expect high levels of two toxins to have an even more adverse effect on a person’s health than high levels of either of the two toxins. Such an effect may be modeled by including interaction effects between the covariates.

Here, we extend the NPSSL to include interaction effects. We consider only second-order interactions between the covariates, but our model can easily be extended to include even higher-order interactions. We assume that the interaction effects may be decomposed into the sum of bivariate functions of each pair of covariates, yielding the model:

[TABLE]

For the interaction terms, we follow Wei et al. [48] and approximate $f_{kl}$ using the outer product of the basis functions of the interacting covariates:

[TABLE]

where $\bm{\beta}_{kl}=(\beta_{kl11},\dots,\beta_{kl1d^{*}},\beta_{kl21},\dots,\beta_{kld^{*}d^{*}})^{T}\in\mathbb{R}^{d^{*2}}$ is the vector of unknown weights. We let $\widetilde{\bm{X}}_{kl}$ denote the $n\times d^{*2}$ matrix with rows

[TABLE]

where $\bm{g}_{k}(X_{ik})=(g_{k1}(X_{ik}),\dots,g_{kd^{*}}(X_{ik}))^{T}$. Then, (5.5) may be represented in matrix form as

[TABLE]

where $\bm{\delta}$ is a vector of the lower-order truncation bias. We again assume $\bm{Y}$ has been centered and so do not include a grand mean in (5.7). We do not constrain $f_{kl}$ to integrate to zero as in Wei et al. [48]. However, we do ensure that the main effects are not in the linear span of the interaction functions. That is, we require the “main effect” matrices $\widetilde{\bm{X}}_{l}$ and $\widetilde{\bm{X}}_{k}$ to be orthogonal to the “interaction” matrix $\widetilde{\bm{X}}_{kl}$. This condition is needed to maintain identifiability for both the main and interaction effects in the model. In practice, we enforce this condition by setting the interaction design matrix to be the residuals of the regression of $\widetilde{\bm{X}}_{k}\circ\widetilde{\bm{X}}_{l}$ on $\widetilde{\bm{X}}_{k}$ and $\widetilde{\bm{X}}_{l}$.

Note that the current representation does not enforce strong hierarchy. That is, interaction terms can be included even if their corresponding main effects are removed from the model. However, the NPSSL model can be easily modified to accommodate strong hierarchy. If hierarchy is desired, the “interaction” matrices can be augmented to contain both main and interaction effects, as in Lim and Hastie [21], i.e. the “interaction” matrices in (5.7) would be $\widetilde{\bm{X}}_{kl}^{\textrm{aug}}=[\widetilde{\bm{X}}_{k},\widetilde{\bm{X}}_{l},\widetilde{\bm{X}}_{kl}]$, instead of simply $\widetilde{\bm{X}}_{kl}$. This augmented model is overparameterized since the main effects still have their own separate design matrices as well (to ensure that main effects can still be selected even if $\bm{\beta}_{kl}^{\textrm{aug}}=\bm{0}$). However, this ensures that interaction effects are only selected if the corresponding main effects are also in the model.

In the interaction model, we either include $\bm{\beta}_{kl}$ in the model (5.7) if there is a non-negligible interaction between the $k$th and $l$th covariates, or we estimate $\widehat{\bm{\beta}}_{kl}=\mathbf{0}_{d^{*2}}$ if such an interaction is negligible. With the non-separable SSGL penalty, the objective function is:

[TABLE]

where $\bm{\beta}=(\bm{\beta}_{1}^{T},\dots,\bm{\beta}_{p}^{T},\bm{\beta}_{12}^{T},\dots\bm{\beta}_{(p-1)p}^{T})^{T}\in\mathbb{R}^{pd+p(p-1)d^{*2}/2}.$ We can again use Algorithm 1 in Appendix A.1 to find the modal estimates of $\bm{\beta}$ and $\sigma^{2}$.

6 Asymptotic Theory for the SSGL and NPSSL

In this section, we derive asymptotic properties for the separable SSGL and NPSSL models. We first note some differences between our theory and the theory in Ročková and George [36]. First, we prove joint consistency in estimation of both the unknown $\bm{\beta}$ and the unknown $\sigma^{2}$, whereas [36] proved their result only for $\bm{\beta}$, assuming known variance $\sigma^{2}=1$. Secondly, Ročková and George [36] established convergence rates for the global posterior mode and the full posterior separately, whereas we establish a contraction rate $\epsilon_{n}$ for the full posterior only. Our rate $\epsilon_{n}$ satisfies $\epsilon_{n}\rightarrow 0$ as $n\rightarrow\infty$ (i.e. the full posterior collapses to the true $(\bm{\beta},\sigma^{2})$ almost surely as $n\rightarrow\infty$), and hence, it automatically follows that the posterior mode is a consistent estimator of $(\bm{\beta},\sigma^{2})$. Finally, we also derive a posterior contraction rate for nonparametric additive regression, not just linear regression. All proofs for the theorems in this section can be found in Appendix D.3.

6.1 Grouped Linear Regression

We work under the frequentist assumption that there is a true model,

[TABLE]

where $\bm{\beta}_{0}=(\bm{\beta}_{01}^{T},\ldots,\bm{\beta}_{0G}^{T})^{T}$ and $\sigma_{0}^{2}\in(0,\infty)$. Denote $\bm{X}=[\bm{X}_{1},\ldots,\bm{X}_{G}]$ and $\bm{\beta}=(\bm{\beta}_{1}^{T},\ldots,\bm{\beta}_{G}^{T})^{T}$. Suppose we endow $(\bm{\beta},\sigma^{2})$ under model (6.1) with the following prior:

[TABLE]

where $c_{0}>0$ and $d_{0}>0$ are fixed constants and the hyperparameters $(a,b)$ in the prior on $\theta$ are to be chosen later.

Remark 1**.**

In our implementation of the SSGL model, we endowed $\sigma^{2}$ with an improper prior, $\pi(\sigma^{2})\propto\sigma^{-2}$. This can be viewed as a limiting case of the $\mathcal{IG}(c_{0},d_{0})$ prior with $c_{0}\rightarrow 0,d_{0}\rightarrow 0$. This improper prior is fine for implementation since it leads to a proper posterior, but for our theoretical investigation, we require the priors on $(\bm{\beta},\sigma^{2})$ to be proper.

6.1.1 Posterior Contraction Rates

Let $m_{\max}=\max_{1\leq j\leq G}m_{g}$ and let $p=\sum_{g=1}^{G}m_{g}$. Let $S_{0}$ be the set containing the indices of the true nonzero groups, where $S_{0}\subseteq\{1,\ldots,G\}$ with cardinality $s_{0}=\lvert S_{0}\rvert$. We make the following assumptions:

1. (A1)

Assume that $G\gg n$, $\log(G)=o(n)$, and $m_{\max}=O(\log G/\log n)$. 2. (A2)

The true number of nonzero groups satisfies $s_{0}=o(n/\log G)$. 3. (A3)

There exists a constant $k>0$ so that $\lambda_{\max}(\bm{X}^{T}\bm{X})\leq kn^{\alpha}$, for some $\alpha\in[1,\infty)$. 4. (A4)

Let $\xi\subset\{1,\ldots,G\}$, and let $\bm{X}_{\xi}$ denote the submatrix of $\bm{X}$ that contains the submatrices with groups indexed by $\xi$. There exist constants $\nu_{1}>0$, $\nu_{2}>0$, and an integer $\bar{p}$ satisfying $s_{0}=o(\bar{p})$ and $\bar{p}=o(s_{0}\log n)$, so that $n\nu_{1}\leq\lambda_{\min}(\bm{X}_{\xi}^{T}\bm{X}_{\xi})\leq\lambda_{\max}(\bm{X}_{\xi}^{T}\bm{X}_{\xi})\leq n\nu_{2}$ for any model of size $\lvert\xi\rvert\leq\bar{p}$. 5. (A5)

$\lVert\bm{\beta}_{0}\rVert_{\infty}=O(\log G).$

Assumption (A1) allows the number of groups $G$ and total number of covariates $p$ to grow at nearly exponential rate with sample size $n$. The size of each individual group may also grow as $n$ grows, but should grow at a slower rate than $n/\log n$. Assumption (A2) specifies the growth rate for the true model size $s_{0}$. Assumption (A3) bounds the eigenvalues of $\bm{X}^{T}\bm{X}$ from above and is less stringent than requiring all the eigenvalues of the Gram matrix ($\bm{X}^{T}\bm{X}/n$) to be bounded away from infinity. Assumption (A4) ensures that $\bm{X}^{T}\bm{X}$ is locally invertible over sparse sets. In general, conditions (A3)-(A4) are difficult to verify, but they can be shown to hold with high probability for certain classes of matrices where the rows of $\bm{X}$ are independent and sub-Gaussian [25, 30]. Finally, Assumption (A5) places a restriction on the growth rate of the maximum signal size for the true $\bm{\beta}_{0}$.

We now state our main theorem on the posterior contraction rates for the SSGL prior (6.2) under model (6.1). Let $\mathbb{P}_{0}$ denote the probability measure underlying the truth (6.1) and $\Pi(\cdot|\bm{Y})$ denote the posterior distribution under the prior (6.2) for $(\bm{\beta},\sigma^{2})$.

Theorem 2** (posterior contraction rates).**

Let $\epsilon_{n}=\sqrt{s_{0}\log G/n}$, and suppose that Assumptions (A1)-(A5) hold. Under model (6.1), suppose that we endow $(\bm{\beta},\sigma^{2})$ with the prior (6.2). For the hyperparameters in the $\mathcal{B}(a,b)$ prior on $\theta$, we choose $a=1,b=G^{c}$, $c>2$. Further, we set $\lambda_{0}=(1-\theta)/\theta$ and $\lambda_{1}\asymp 1/n$ in the SSGL prior. Then

[TABLE]

[TABLE]

[TABLE]

for some $M_{1}>0,M_{2}>0$.

Remark 2**.**

In the case where $G=p$ and $m_{1}=\ldots=m_{G}=1$, the $\ell_{2}$ and prediction error rates in (6.3)-(6.4) reduce to the familiar optimal rates of $\sqrt{s_{0}\log p/n}$ and $\sqrt{s_{0}\log p}$ respectively.

Remark 3**.**

Eq. (6.5) demonstrates that our model also consistently estimates the unknown variance $\sigma^{2}$, therefore providing further theoretical justification for placing an independent prior on $\sigma^{2}$, as advocated by Moran et al. [28].

6.1.2 Dimensionality Recovery

Although the posterior mode is exactly sparse, the SSGL prior is absolutely continuous so it assigns zero mass to exactly sparse vectors. To approximate the model size under the SSGL model, we use the following generalized notion of sparsity [3]. For $\omega_{g}>0$, we define the generalized inclusion indicator and generalized dimensionality, respectively, as

[TABLE]

In contrast to [3, 36], we allow the threshold $\omega_{g}$ to be different for each group, owing to the fact that the group sizes $m_{g}$ may not necessarily all be the same. However, the $\omega_{g}$’s, $g=1,\ldots,G$, should still tend towards zero as $n$ increases, so that $|\bm{\gamma}(\bm{\beta})|$ provides a good approximation to $\#\{g:\bm{\beta}_{g}\neq\mathbf{0}_{m_{g}}\}$.

Consider as the threshold,

[TABLE]

Note that for large $\lambda_{0}$, this threshold rapidly approaches zero. Analogous to [34, 36], any vectors $\bm{\beta}_{g}$ that satisfy $\lVert\bm{\beta}_{g}\rVert_{2}=\omega_{g}$ correspond to the intersection points between the two group lasso densities in the separable SSGL prior (2.2), or when the second derivative $\partial^{2}pen_{S}(\bm{\beta}|\theta)/\partial\lVert\bm{\beta}_{g}\rVert_{2}^{2}=0.5$. The value $\omega_{g}$ represents the turning point where the slab has dominated the spike, and thus, the sharper the spike (when $\lambda_{0}$ is large), the smaller the threshold.

Using the notion of generalized dimensionality (6.6) with (6.7) as the threshold, we have the following theorem.

Theorem 3** (dimensionality).**

Suppose that the same conditions as those in Theorem 2 hold. Then under (6.1), for sufficiently large $M_{3}>0$,

[TABLE]

Theorem 3 shows that the expected posterior probability that the generalized dimension is a constant multiple larger than the true model size $s_{0}$ is asymptotically vanishing. In other words, the SSGL posterior concentrates on sparse sets.

6.2 Sparse Generalized Additive Models (GAMs)

Assume there is a true model,

[TABLE]

where $\sigma_{0}^{2}\in(0,\infty)$. Throughout this section, we assume that all the covariates $\bm{X}_{i}=(X_{i1},\ldots,X_{ip})^{T}$ have been standardized to lie in $[0,1]^{p}$ and that $f_{0j}\in\mathcal{C}^{\kappa}[0,1],j=1,\ldots,p$. That is, the true functions are all at least $\kappa$-times continuously differentiable over $[0,1]$, for some $\kappa\in\mathbb{N}$. Suppose that each $f_{0j}$ can be approximated by a linear combination of basis functions, $\{g_{j1},\ldots,g_{jd}\}$. In matrix notation, (6.9) can then be written as

[TABLE]

where $\widetilde{\bm{X}}_{j}$ denotes an $n\times d$ matrix where the $(i,k)$th entry is $\widetilde{\bm{X}}_{j}(i,k)=g_{jk}(X_{ij})$, the $\bm{\beta}_{0j}$’s are $d\times 1$ vectors of basis coefficients, and $\bm{\delta}$ denotes an $n\times 1$ vector of lower-order bias.

Denote $\widetilde{\bm{X}}=[\widetilde{\bm{X}}_{1},\ldots,\widetilde{\bm{X}}_{p}]$ and $\bm{\beta}=(\bm{\beta}_{1}^{T},\ldots,\bm{\beta}_{p}^{T})^{T}$. Under (6.9), suppose that we endow $(\bm{\beta},\sigma^{2})$ in (6.10) with the prior (6.2). We have the following assumptions:

1. (B1)

Assume that $p\gg n$, $\log p=o(n)$, and $d\asymp n^{1/(2\kappa+1)}$. 2. (B2)

The number of true nonzero functions satisfies

[TABLE] 3. (B3)

There exists a constant $k_{1}>0$ so that for all $n$, $\lambda_{\max}(\widetilde{\bm{X}}^{T}\widetilde{\bm{X}})\leq k_{1}n$. 4. (B4)

Let $\xi\subset\{1,\ldots,p\}$, and let $\widetilde{\bm{X}}_{\xi}$ denote the submatrix of $\widetilde{\bm{X}}$ that contains the submatrices indexed by $\xi$. There exists a constant $\nu_{1}>0$ and an integer $\bar{p}$ satisfying $s_{0}=o(\bar{p})$ and $\bar{p}=o(s_{0}\log n)$, so that $\lambda_{\min}(\widetilde{\bm{X}}_{\xi}^{T}\widetilde{\bm{X}}_{\xi})\geq n\nu_{1}$ for any model of size $\lvert\xi\rvert\leq\bar{p}$. 5. (B5)

$\lVert\bm{\beta}_{0}\rVert_{\infty}=O(\log p).$ 6. (B6)

The bias $\bm{\delta}$ satisfies $\lVert\bm{\delta}\rVert_{2}\lesssim\sqrt{s_{0}n}d^{-\kappa}$.

Assumptions (B1)-(B5) are analogous to assumptions (A1)-(A5). Assumptions (B3)-(B4) are difficult to verify but can be shown to hold if appropriate basis functions for the $g_{jk}$’s are used, e.g. cubic B-splines [52, 48]. Finally, Assumption (B6) bounds the approximation error incurred by truncating the basis expansions to be of size $d$. This assumption is satisfied, for example, by B-spline basis expansions [56, 48].

Let $\widetilde{\mathbb{P}}_{0}$ denote the probability measure underlying the truth (6.9) and $\Pi(\cdot|\bm{Y})$ denote the posterior distribution under NPSSL model with the prior (6.2) for $(\bm{\beta},\sigma^{2})$ in (6.10). Further, let $f(\bm{X}_{i})=\sum_{j=1}^{p}f_{j}(X_{ij})$ and $f_{0}(\bm{X}_{i})=\sum_{j=1}^{p}f_{0j}(X_{ij})$, and define the empirical norm $\lVert\cdot\rVert_{n}$ as

[TABLE]

Let $\mathcal{F}$ denote the infinite-dimensional set of all possible additive functions $f=\sum_{j=1}^{p}f_{j}$, where each $f_{j}$ can be represented by a $d$-dimensional basis expansion. In Raskutti et al. [31], it was shown that the minimax estimation rate for $f_{0}=\sum_{j=1}^{p}f_{0j}$ under squared $\ell_{2}$ error loss is $\epsilon_{n}^{2}\asymp s_{0}\log p/n+s_{0}n^{-2\kappa/(2\kappa+1)}$. The next theorem establishes that the NPSSL model achieves this minimax posterior contraction rate.

Theorem 4** (posterior contraction rates).**

Let $\epsilon_{n}^{2}=s_{0}\log p/n+s_{0}n^{-2\kappa/(2\kappa+1)}$. Suppose that Assumptions (B1)-(B6) hold. Under model (6.10), suppose that we endow $(\bm{\beta},\sigma^{2})$ with the prior (6.2) (replacing $G$ with $p$). For the hyperparameters in the $\mathcal{B}(a,b)$ prior on $\theta$, we choose $a=1,b=p^{c}$, $c>2$. Further, we set $\lambda_{0}=(1-\theta)/\theta$ and $\lambda_{1}\asymp 1/n$ in the SSGL prior. Then

[TABLE]

[TABLE]

for some $\widetilde{M}_{1}>0$.

Let the generalized dimensionality $|\bm{\gamma}(\bm{\beta})|$ be defined as before in (6.6) (replacing $G$ with $p$), with $\omega_{g}$ from (6.7) as the threshold (replacing $m_{g}$ with $d$). The next theorem shows that under the NPSSL, the expected posterior probability that the generalized dimension size is a constant multiple larger than the true model size $s_{0}$ asymptotically vanishes.

Theorem 5** (dimensionality).**

Suppose that the same conditions as those in Theorem 4 hold. Then under (6.10), for sufficiently large $\widetilde{M}_{2}>0$,

[TABLE]

7 Simulation Studies

In this section, we will evaluate our method in a number of settings. For the SSGL approach, we fix $\lambda_{1}=1$ and use cross-validation to choose from $\lambda_{0}\in\{1,2,\ldots,100\}$. For the prior $\theta\sim\mathcal{B}(a,b)$, we set $a=1,b=G$ so that $\theta$ is small with high probability. We will compare our SSGL approach with the following methods:

GroupLasso: the group lasso [54] 2. 2.

BSGS: Bayesian sparse group selection [6] 3. 3.

SoftBart: soft Bayesian additive regression tree (BART) [22] 4. 4.

RandomForest: random forests [5] 5. 5.

SuperLearner: super learner [47] 6. 6.

GroupSpike: point-mass spike-and-slab priors (1.6) placed on groups of coefficients555Code to implement GroupSpike is included in the Supplementary data. Due to the discontinuous prior, GroupSpike is not amenable to a MAP finding algorithm and has to be implemented using MCMC.

In our simulations, we will look at the mean squared error (MSE) for estimating $f(\bm{X}_{\textrm{new}})$ averaged over a new sample of data $\bm{X}_{\textrm{new}}$. We will also evaluate the variable selection properties of the different methods using precision and recall, where $\text{precision}=\text{TP}/(\text{TP}+\text{FP})$, $\text{recall}=\text{TP}/(\text{TP}+\text{FN})$, and TP, FP, and FN denote the number of true positives, false positives, and false negatives respectively. Note that we will not show precision or recall for the SuperLearner, which averages over different models and different variable selection procedures and therefore does not have one set of variables that are deemed significant.

7.1 Sparse Semiparametric Regression

Here, we will evaluate the use of our proposed SSGL procedure in sparse semiparametric regression with $p$ continuous covariates. Namely, we implement the NPSSL main effects model described in Section 5.1. In Appendix B, we include more simulation studies of the SSGL approach under both sparse and dense settings, as well as a simulation study showing that we are accurately estimating the residual variance $\sigma^{2}$.

We let $n=100,p=300$. We generate independent covariates from a standard uniform distribution, and we let the true regression surface take the following form:

[TABLE]

with variance $\sigma^{2}=1$.

To implement the SSGL approach, we estimate the mean response as

[TABLE]

where $\widetilde{\bm{X}}_{j}$ is a design matrix of basis functions used to capture the possibly nonlinear effect of $X_{j}$ on $Y$. For the basis functions in $\widetilde{\bm{X}}_{j},j=1,\ldots,p$, we use natural splines with degrees of freedom $d$ chosen from $d\in\{2,3,4\}$ using cross-validation. Thus, we are estimating a total of between 600 and 1200 unknown basis coefficients.

We run 1000 simulations and average all of the metrics considered over each simulated data set. Figure 1 shows the results from this simulation study. The GroupSpike approach has the best performance in terms of MSE, followed closely by SSGL, with the next best approach being SoftBart. In terms of recall, the SSGL and GroupLasso approaches perform the best, indicating the highest power in detecting the significant groups. This comes with a loss of precision as the GroupSpike and SoftBart approaches have the best precision among all methods.

Although the GroupSpike method performed best in this scenario, the SSGL method was much faster. As we show in Appendix B.5, when $p=4000$, fitting the SSGL model with a sufficiently large $\lambda_{0}$ takes around three seconds to run. This is almost 50 times faster than running 100 MCMC iterations of the GroupSpike method (never mind the total time it takes for the GroupSpike model to converge). Our experiments demonstrate that the SSGL model gives comparable performance to the “theoretically ideal” point mass spike-and-slab in a fraction of the computational time.

7.2 Interaction Detection

We now explore the ability of the SSGL approach to identify important interaction terms in a nonparametric regression model. To this end, we implement the NPSSL model with interactions from Section 5.2. We generate 25 independent covariates from a standard uniform distribution with a sample size of 300. Data is generated from the model:

[TABLE]

with variance $\sigma^{2}=1$. While this may not seem like a high-dimensional problem, we will consider all two-way interactions, and there are 300 such interactions. The important two-way interactions are between $X_{1}$ and $X_{2}$ and between $X_{3}$ and $X_{5}$. We evaluate the performance of each method and examine the ability of SSGL to identify important interactions while excluding all of the remaining interactions. Figure 2 shows the results for this simulation setting. The SSGL, GL, GroupSpike, and SoftBart approaches all perform well in terms of out-of-sample mean squared error, with GroupSpike slightly outperforming the competitors. The SSGL also does a very good job at identifying the two important interactions. The $(X_{1},X_{2})$ interaction is included in 97% of simulations, while the $(X_{3},X_{5})$ interaction is included 100% of the time. All other interactions are included in only a small fraction of simulated data sets.

8 Real Data Analysis

Here, we will illustrate the SSGL procedure in two distinct settings: 1) evaluating the SSGL’s performance on a data set where $n=120$ and $p=15,000$, and 2) identifying important (nonlinear) main effects and interactions of environmental exposures. In Appendix C, we evaluate the predictive performance of our approach on benchmark data sets where $p<n$, compared to several other state-of-the-other methods. Our results show that in both the $p\gg n$ and $p<n$ settings, the SSGL maintains good predictive accuracy.

8.1 Bardet-Biedl Syndrome Gene Expression Study

We now analyze a microarray data set consisting of gene expression measurements from the eye tissue of 120 laboratory rats666Data accessed from the Gene Expression Omnibus www.ncbi.nlm.nih.gov/geo (accession no. GSE5680).. The data was originally studied by Scheetz et al. [37] to investigate mammalian eye disease, and later analyzed by Breheny and Huang [4] to demonstrate the performance of their group variable selection algorithm. In this data, the goal is to identify genes which are associated with the gene TRIM32. TRIM32 has previously been shown to cause Bardet-Biedl syndrome [7], a disease affecting multiple organs including the retina.

The original data consists of 31,099 probe sets. Following Breheny and Huang [4], we included only the 5,000 probe sets with the largest variances in expression (on the log scale). For these probe sets, we considered a three-term natural cubic spline basis expansion, resulting in a grouped regression problem with $n=120$ and $p=15,000$. We implemented SSGL with regularization parameter values $\lambda_{1}=1$ and $\lambda_{0}$ ranging on an equally spaced grid from 1 to 500. We compared SSGL with the group lasso [54], implemented using the R package gglasso [51].

As shown in Table 2, SSGL selected much fewer groups than the group lasso. Namely, SSGL selected 12 probe sets, while the group lasso selected 83 probe sets. Moreover, SSGL achieved a smaller 10-fold cross-validation error than the group lasso, albeit within range of random variability (Table 2). These results demonstrate that the SSGL achieves strong predictive accuracy, while also achieving the most parsimony. The groups selected by both SSGL and the group lasso are displayed in Table 4 of Appendix C. Interestingly, only four of the 12 probes selected by SSGL were also selected by the group lasso.

We next conducted gene ontology enrichment analysis on the group of genes found by each of the methods using the R package clusterProfiler [53]. This software determines whether subsets of genes known to act in a biological process are overrepresented in a group of genes, relative to chance. If such a subset is significant, the group of genes is said to be “enriched” for that biological process. With a false discovery rate of 0.01, SSGL had five enriched terms, while the group lasso had none. The terms for which SSGL was enriched included RNA binding, a biological process with which the response gene TRIM32 is associated.777https://www.genecards.org/cgi-bin/carddisp.pl?gene=TRIM32 (accessed 03/01/20) These findings show the ability of SSGL to find biologically meaningful signal in the data. Additional details for our gene ontology enrichment analysis can be found in Appendix C.

8.2 Environmental Exposures in the NHANES Data

Here, we analyze data from the 2001-2002 cycle of the National Health and Nutrition Examination Survey (NHANES), which was previously analyzed by Antonelli et al. [1]. We aim to identify which organic pollutants are associated with changes in leukocyte telomere length (LTL) levels. Telomeres are segments of DNA that help to protect chromosomes, and LTL levels are commonly used as a proxy for overall telomere length. LTL levels have previously been shown to be associated with adverse health effects [13], and recent studies within the NHANES data have found that organic pollutants can be associated with telomere length [26].

We use the SSGL approach to evaluate whether any of 18 organic pollutants are associated with LTL length and whether there are any significant interactions among the pollutants also associated with LTL length. In addition to the 18 exposures, there are 18 additional demographic variables which we adjust for in our model. We model the effects of the 18 exposures on LTL length using spline basis functions with two degrees of freedom. For the interaction terms, this leads to four terms for each pair of interactions, and we orthogonalize these terms with respect to the main effects. In total, this leads to a data set with $n=1003$ and $p=666$.

Our model selects four significant main effects and six significant interaction terms. In particular, PCB 3, PCB 11, Furan 1, and Furan 4 are identified as the important main effects in the model. Figure 3 plots the exposure response curves for these exposures. We see that each of these four exposures has a positive association with LTL length, which agrees with results seen in [26] that saw positive relationships between persistent organic pollutants and telomere length. Further, our model identifies more main effects and more interactions than previous analyses of these data, e.g. [1], which could lead to more targeted future research in understanding how these pollutants affect telomere length. Additional discussion and analysis of the NHANES data set can be found in Appendix C.

9 Discussion

We have introduced the spike-and-slab group lasso (SSGL) model for variable selection and linear regression with grouped variables. We also extended the SSGL model to generalized additive models with the nonparametric spike-and-slab lasso (NPSSL). The NPSSL can efficiently identify both nonlinear main effects and higher-order nonlinear interaction terms. Moreover, our prior performs an automatic multiplicity adjustment and self-adapts to the true sparsity pattern of the data through a non-separable penalty. For computation, we introduced highly efficient coordinate ascent algorithms for MAP estimation and employed de-biasing methods for uncertainty quantification. An R package implementing the SSGL model can be found at https://github.com/jantonelli111/SSGL.

Although our model performs group selection, it does so in an “all-in-all-out” manner, similar to the original group lasso [54]. Future work will be to extend our model to perform both group selection and within-group selection of individual coordinates. We are currently working to extend the SSGL to perform bilevel selection.

We are also working to extend the nonparametric spike-and-slab lasso so it can adapt to even more flexible regression surfaces than the generalized additive model. Under the NPSSL model, we used cross-validation to tune a single value for the degrees of freedom. In reality, different functions can have vastly differing degrees of smoothness, and it will be desirable to model anisotropic regression surfaces while avoiding the computational burden of tuning the individual degrees of freedom over a $p$-dimensional grid.

Acknowledgments

Dr. Ray Bai, Dr. Gemma Moran, and Dr. Joseph Antonelli contributed equally and wrote this manuscript together, with input and suggestions from all other listed co-authors. The bulk of this work was done when the first listed author was a postdoc at the Perelman School of Medicine, University of Pennsylvania, under the mentorship of the last two authors. The authors are grateful to three anonymous reviewers, the Associate Editor, and the Editor whose thoughtful comments and suggestions helped to improve this manuscript. The authors would also like to thank Ruoyang Zhang, Peter Bühlmann, and Edward George for helpful discussions.

Funding

Dr. Ray Bai and Dr. Mary Boland were funded in part by generous funding from the Perelman School of Medicine, University of Pennsylvania. Dr. Ray Bai and Dr. Yong Chen were funded by NIH grants 1R01AI130460 and 1R01LM012607.

Appendix A Additional Computational Details

A.1 SSGL Block-Coordinate Ascent Algorithm

A.2 Tuning Hyperparameters, Initializing Values, and Updating the Variance in Algorithm 1

We keep the slab hyperparameter $\lambda_{1}$ fixed at a small value. We have found that our results are not very sensitive to the choice of $\lambda_{1}$. This parameter controls the variance of the slab component of the prior, and the variance must simply be large enough to avoid overshrinkage of important covariates. For the default implementation, we recommend fixing $\lambda_{1}=1$. This applies minimal shrinkage to the significant groups of coefficients and affords these groups the ability to escape the pull of the spike.

Meanwhile, we choose the spike parameter $\lambda_{0}$ from an increasing ladder of values. We recommend selecting $\lambda_{0}\in\{1,2,...,100\}$, which represents a range from hardly any penalization to very strong penalization. Below, we describe precisely how to tune $\lambda_{0}$. To account for potentially different group sizes, we use the same $\lambda_{0}$ for all groups but multiply $\lambda_{0}$ by $\sqrt{m_{g}}$ for each $g$th group, $g=1,\ldots,G$. As discussed in [14], further scaling of the penalty by group size is necessary in order to ensure that the same degree of penalization is applied to potentially different sized groups. Otherwise, larger groups may be erroneously selected simply because they are larger (and thus have larger $\ell_{2}$ norm), not because they contain significant entries.

When the spike parameter $\lambda_{0}$ is very large, the continuous spike density approximates the point-mass spike. Consequently, we face the computational challenge of navigating a highly multimodal posterior. To ameliorate this problem for the spike-and-slab lasso, Ročková and George [36] recommend a “dynamic posterior exploration” strategy in which the slab parameter $\lambda_{1}$ is held fixed at a small value and $\lambda_{0}$ is gradually increased along a grid of values. Using the solution from a previous $\lambda_{0}$ as a “warm start” allows the procedure to more easily find optimal modes. In particular, when $(\lambda_{1}-\lambda_{0})^{2}\leq 4$, the posterior is convex.

Moran et al. [28] modify this strategy for the unknown $\sigma^{2}$ case. This is because the posterior is always non-convex when $\sigma^{2}$ is unknown. Namely, when $p\gg n$ and $\lambda_{0}\approx\lambda_{1}$, the model can become saturated, causing the residual variance to go to zero. To avoid this suboptimal mode at $\sigma^{2}=0$, Moran et al. [28] recommend fixing $\sigma^{2}$ until the $\lambda_{0}$ value at which the algorithm starts to converge in less than 100 iterations. Then, $\bm{\beta}$ and $\sigma^{2}$ are simultaneously updated for the next largest $\lambda_{0}$ in the sequence. The intuition behind this strategy is we first find a solution to the convex problem (in which $\sigma^{2}$ is fixed) and then use this solution as a warm start for the non-convex problem (in which $\sigma^{2}$ can vary).

We pursue a similar “dynamic posterior exploration” strategy with the modification for the unknown variance case for the SSGL in Algorithm 1 of Section A.1. A key aspect of this algorithm is how to choose the maximum value of $\lambda_{0}$. Ročková and George [36] recommend this maximum to be the $\lambda_{0}$ value at which the estimated coefficients stabilize. An alternative approach is to choose the maximum $\lambda_{0}$ using cross-validation, a strategy which is made computationally feasible by the speed of our block coordinate ascent algorithm. In our experience, the dynamic posterior exploration strategy favors more parsimonious models than cross-validation. In the simulation studies in Section 7, we utilize cross-validation to choose $\lambda_{0}$, as there, our primary goal is predictive accuracy rather than parsimony.

Following [28], we initialize $\bm{\beta}^{*}=\bm{0}_{p}$ and $\theta^{*}=0.5$. We also initialize $\sigma^{*2}$ to be the mode of a scaled inverse chi-squared distribution with degrees of freedom $\nu=3$ and scale parameter chosen such that the sample variance of $\bm{Y}$ corresponds to the 90th quantile of the prior. We have found this initialization to be quite effective in practice at ensuring that Algorithm 1 converges in less than 100 iterations for sufficiently large $\lambda_{0}$.

A.3 Additional Details for the Inference Procedure

Here, we describe the nodewise regression procedure for estimating $\widehat{\bm{\Theta}}$ in Section 4. This approach for estimating the inverse of the covariance matrix $\widehat{\bm{\Sigma}}=\bm{X}^{T}\bm{X}/n$ was originally proposed and studied theoretically in [24] and [46].

For each $j=1,\dots,p$, let $\bm{X}_{j}$ denote the $j$th column of $\bm{X}$ and $\bm{X}_{-j}$ denote the submatrix of $\bm{X}$ with the $j$th column removed. Define $\widehat{\bm{\gamma}}_{j}$ as

[TABLE]

Now we can define the components of $\widehat{\bm{\gamma}}_{j}$ as $\widehat{\bm{\gamma}}_{j,k}$ for $k=1,\dots,p$ and $k\neq p$, and create the following matrix:

[TABLE]

Lastly, let $\widehat{\bm{T}}^{2}=\text{diag}(\widehat{\tau}_{1}^{2},\widehat{\tau}_{2}^{2},\dots,\widehat{\tau}_{p}^{2})$, where

[TABLE]

We can proceed with $\widehat{\bm{\Theta}}=\widehat{\bm{T}}^{-2}\widehat{\bm{C}}$. This choice is used because it puts an upper bound on $||\widehat{\bm{\Sigma}}\widehat{\bm{\Theta}}_{j}^{T}-\bm{e}_{j}||_{\infty}$. Other regression models such as the original spike-and-slab lasso [36] could be used instead of the lasso [45] regressions for each covariate. However, we will proceed with this choice, as it has already been studied theoretically and shown to have the required properties to be able to perform inference for $\bm{\beta}$.

Appendix B Additional Simulation Results

Here, we present additional results which include different sample sizes than those seen in the manuscript, assessment of the SSGL procedure under dense settings, estimates of $\sigma^{2}$, timing comparisons, and additional figures.

B.1 Increased Sample Size for Sparse Simulation

Here, we present the same sparse simulation setup as that seen in Section 7.1, though we will increase $n$ from 100 to 300. Figure 4 shows the results and we see that they are very similar to those from the manuscript, except that the mean squared error (MSE) for the SSGL approach is now nearly as low as the MSE for the GroupSpike approach, and the precision score has improved substantially.

B.2 Dense Model

Here, we generate independent covariates from a standard normal distribution, and we let the true regression surface take the following form

[TABLE]

with variance $\sigma^{2}=1$. In this model, there are no strong predictors of the outcome, but rather a large number of predictors which have small impacts on the outcome. Here, we display results for both $n=100$ and $p=300$, as well as $n=300$ and $p=300$, as the qualitative results change across the different sample sizes. Our simulation results can be seen in Figure 5. When the sample size is 100, the SSGL procedure performs the best in terms of both MSE and recall score, while all approaches do poorly with the precision score. When the sample size increases to 300, the SSGL approach still performs quite well in terms of MSE and recall, though the GroupLasso and GroupSpike approaches are slightly better in terms of MSE. The SSGL approach still maintains a low precision score while the GroupSpike approach has a very high precision once the sample size is large enough.

B.3 Estimation of $\sigma^{2}$

To evaluate our ability to estimate $\sigma^{2}$ and confirm our theoretical results that the posterior of $\sigma^{2}$ contracts around the true parameter, we ran a simulation study using the following data generating model:

[TABLE]

with $\sigma^{2}=1$. We vary $n\in\{50,100,500,1000,2000\}$ and we set $G=n$ to confirm that the estimates are centering around the truth as both the sample size and covariate dimension grows. We use groups of size two that contain both the linear and quadratic term for each covariate. Note that in this setting, the total number of regression coefficients actually exceeds the sample size since each group has two terms, leading to a total of $p=2G$ coefficients in the model.

Figure 6 shows box plots of the estimates for $\sigma^{2}$ across all simulations for each sample size and covariate dimension. We see that for small sample sizes there are some estimates well above 1 or far smaller than 1. This is because either some important variables are excluded (so the sum of squared residuals gets inflated), or too many variables are included and the model is overfitted (leading to small $\widehat{\sigma}^{2}$). These problems disappear as the sample size grows to 500 or larger, where we observe that the estimates are closely centering around the true $\sigma^{2}=1$. Figure 6 confirms our theoretical results in Theorem 2 and Theorem 4, which state that as $n,G\rightarrow\infty$, the posterior $\pi(\sigma^{2}|\bm{Y})$ contracts around the true $\sigma^{2}$.

B.4 Large Number of Groups

We now generate data with $n=200$ and $G=2000$, where each group contains three predictors. We generate data from the following model:

[TABLE]

where we set $\bm{\beta}_{g}=\bm{0}$ for $g=1,\dots 1996$. For the final four groups, we draw individual coefficient values from independent normal distributions with mean 0 and standard deviation 0.4. These coefficients are redrawn for each data set in the simulation study, and therefore, the results are averaging over many possible combinations of magnitudes for the true nonzero coefficients. We see that the best performing approach in this scenario is the GroupSpike approach, followed by the SSGL approach. The SSGL approach outperforms group lasso in terms of MSE and precision, while group lasso has a slightly higher recall score.

B.5 Computation Time

In this study, we evaluate the computational speed of the SSGL procedure in comparison with the fully Bayesian GroupSpike approach that places point-mass spike-and-slab priors on groups of coefficients. We fix $n=300$ and vary the number of groups $G\in\{100,200,\dots,2000\}$, with two elements per group. For the SSGL approach, we keep track of the computation time for estimating the model for $\lambda_{0}=20$. For large values of $\lambda_{0}$, it typically takes 100 or fewer iterations for the SSGL method to converge. For the GroupSpike approach, we keep track of the computation time required to run 100 MCMC iterations. Both SSGL and GroupSpike models were run on an Intel E5-2698v3 processor.

In any given data set, the computation time will be higher than the numbers presented here because the SSGL procedure typically requires fitting the model for multiple values of $\lambda_{0}$, while the GroupSpike approach will likely take far more than 100 MCMC iterations to converge, especially in higher dimensions. Nonetheless, this should provide a comparison of the relative computation speeds for each approach.

The average CPU time in seconds can be found in Figure 8. We see that the SSGL approach is much faster as it is able to estimate all the model parameters for a chosen $\lambda_{0}$ in just a couple of seconds, even for $G=2000$ (or $p=4000$). When $p=2000$, the SSGL returned a final solution in roughly three seconds on average, whereas GroupSpike required over two minutes to run 100 iterations (and would most likely require many more iterations to converge). This is to be expected as the GroupSpike approach relies on MCMC. Figure 8 shows the large computational gains that can be achieved using our MAP finding algorithm.

Appendix C Additional Results and Discussion for Real Data Examples

In this section, we perform additional data analysis of the SSGL method on benchmark datasets where $p<n$ to demonstrate that the SSGL model also works well in low-dimensional settings. We also provide additional analyses and discussion of the two real data examples analyzed in Section 8.

C.1 Testing Predictive Performance of the SSGL on Datasets Where $p<n$

We first look at three data sets which have been analyzed in a number of manuscripts, most recently in [22]. The tecator data set is available in the caret package in R [17] and has three different outcomes $\bm{Y}$ to analyze. Specifically, this data set looks at using 100 infrared absorbance spectra to predict three different features of chopped meat with a sample size of 215. The Blood-Brain data is also available in the caret package and aims to predict levels of a particular compound in the brain given 134 molecular descriptors with a sample size of 208. Lastly, the Wipp data set contains 300 samples with 30 features from a computer model used to estimate two-phase fluid flow [41]. For each of these data sets, we hold out 20 of the subjects in the data as a validation sample and see how well the model predicts the outcome in the held-out data. We repeat this 1000 times and compute the root mean squared error (RMSE) for prediction.

Table 3 shows the results for each of the methods considered in the simulation study. The results are standardized so that for each data set, the RMSE is divided by the minimum RMSE for that data set. This means that the model with an RMSE of 1 had the best predictive performance, and all others should be greater than 1, with the magnitude indicating how poor the performance was. We see that the top performer across the data sets was SuperLearner, which is not surprising given that SuperLearner is quite flexible and averages over many different prediction models. Our simulation studies showed that SuperLearner may not work as well when $p>n$. However, the data sets considered here all have $p<n$, which could explain its improved performance here. Among the other approaches, SSGL performs quite well as it has RMSE’s relatively close to 1 for all the data sets considered.

C.2 Additional Details for Bardet-Biedl Analysis

Here we present additional results for the Bardet-Biedl Syndrome gene expression analysis conducted in Section 8.1. Table 4 displays the 12 probes found by SSGL. Table 5 displays the terms for which SSGL was enriched in a gene ontology enrichment analysis.

C.3 Additional Details for Analysis of NHANES Data

Here we will present additional results from the NHANES data analysis in Section 8.2. Here, the aim is to identify environmental exposures that are associated with leukocyte telomere length. In the NHANES data, we have measurements from 18 persistent organic pollutants. Persistent organic pollutants are toxic chemicals that have potential to adversely affect health. They are known to remain in the environment for long periods of time and can travel through wind, water, or even the food chain. Our data set consists of 11 polychlorinated biphenyls (PCBs), three Dioxins, and four Furans. We want to understand the impact that these can have on telomere length, and to understand if any of these pollutants interact in their effect on humans.

The data also contains additional covariates that we will adjust for such as age, a squared term for age, gender, BMI, education status, race, lymphocyte count, monocyte count, cotinine level, basophil count, eosinophil count, and neutrophil count. To better understand the data set, we have shown the correlation matrix between all organic pollutants and covariates in Figure 9. We can see that the environmental exposures are all fairly positively correlated with each other. In particular, the PCBs are highly correlated among themselves. The correlation across chemical types, such as the correlation between PCBs and Dioxins or Furans are lower, though still positively correlated. The correlation between the covariates that we place into our model and the exposures is generally extremely low, and the correlation among the individual covariates is also low, with the exception of a few blood cell types as seen in the upper right of Figure 9.

As discussed in Section 8.2, when we fit the SSGL model to this data set, we identified four main effects (plotted in Figure 3). Our model also identified six interactions as having nonzero parameter estimates. The identified interactions are PCB 10 - PCB 7, Dioxin 1 - PCB 11, Dioxin 2 - PCB 2, Dioxin 2 - Dioxin 1, Furan 1 - PCB 10, and Furan 4 - Furan 3. We see that there are interactions both within a certain type of pollutant (Dioxin and Dioxin, etc.) and across pollutant types (Furan and PCB).

Lastly, looking at Figure 3, we can see that the exposure response curves for the four identified main effects are relatively linear, particularly for PCB 11 and Furan 1. With this in mind, we also ran our SSGL model with one degree of freedom splines for each main effect. Note that this does not require a model that handles grouped covariate structures as the main effects and interactions in this case are both single parameters. Cross-validated error from the model with one degree of freedom is nearly identical to the model with two degrees of freedom, though the linear model selects far more terms. The linear model selects six main effect terms and 20 interaction terms. As the two models provide similar predictive performance but the model with two degrees of freedom is far more parsimonious, we elect to focus on the model with two degrees of freedom.

Appendix D Proofs of Main Results

D.1 Preliminary Lemmas

Before proving the main results in the paper, we first prove the following lemmas.

Lemma 4**.**

Suppose that $\bm{\beta}_{g}\in\mathbb{R}^{m_{g}}$ follows a group lasso density indexed by $\lambda$, i.e. $\bm{\beta}_{g}\sim\bm{\Psi}(\bm{\beta}_{g}|\lambda)$. Then

[TABLE]

Proof.

The group lasso density, $\bm{\Psi}(\bm{\beta}_{g}|\lambda)$, is the marginal density of a scale mixture,

[TABLE]

Therefore, using iterated expectations, we have

[TABLE]

∎

Lemma 5**.**

Suppose $\sigma^{2}>0,\sigma_{0}^{2}>0$. Then for any $\epsilon_{n}\in(0,1)$ such that $\epsilon_{n}\rightarrow 0$ as $n\rightarrow\infty$, we have for sufficiently large $n$,

[TABLE]

Proof.

For large $n$, $\epsilon_{n}<1/2$, so $2\epsilon_{n}/(1-\epsilon_{n})<4\epsilon_{n},-2\epsilon_{n}/(1+\epsilon_{n})>-4\epsilon_{n}$, and thus,

[TABLE]

and hence,

[TABLE]

. ∎

Lemma 6**.**

Suppose that a vector $\bm{z}\in\mathbb{R}^{m}$ can be decomposed into subvectors, $\bm{z}=[\bm{z}_{1}^{\prime},\ldots,\bm{z}_{d}^{\prime}]$, where $\sum_{i=1}^{d}\lvert\bm{z}_{i}\rvert=m$ and $\lvert\bm{z}_{i}\rvert$ denotes the length of $\bm{z}_{i}$. Then $\lVert\bm{z}\rVert_{2}\leq\sum_{i=1}^{d}\lVert\bm{z}_{i}\rVert_{2}$.

Proof.

We have

[TABLE]

∎

D.2 Proofs for Section 3

Proof of Proposition 2.

This result follows from an adaptation of the arguments of Zhang and Zhang [55]. The group-specific optimization problem is:

[TABLE]

We first note that the optimization problem (D.1) is equivalent to maximizing the objective

[TABLE]

where $\varphi$ is the angle between $\bm{z}_{g}$ and $\bm{\beta}_{g}$. Then, when $\lVert\bm{z}_{g}\rVert_{2}<\Delta$, the second factorized term of (D.4) is always less than zero, and so $\widehat{\bm{\beta}}_{g}=\mathbf{0}_{m_{g}}$ must be the global maximizer of $L$. On the other hand, when the global maximizer $\widehat{\bm{\beta}}_{g}=\mathbf{0}_{m_{g}}$, then the second factorized term must always be less than zero, otherwise $\widehat{\bm{\beta}}_{g}=\mathbf{0}_{m_{g}}$ would no longer be the global maximizer and so $\lVert\bm{z}_{g}\rVert_{2}<\Delta$. ∎

Proof of Lemma 3.

We have

[TABLE]

When $\lambda_{0}\to\infty$, we have $z\to 1$ and $x_{g}\to-\infty$ for all $g=1,\dots,\widehat{q}$. Hence,

[TABLE]

∎

D.3 Proofs for Section 6

In this section, we use proof techniques from [29, 40, 48] rather than the ones in [36]. However, none of these other papers considers both continuous spike-and-slab priors for groups of regression coefficients and an independent prior on the unknown variance.

Proof of Theorem 2.

Our proof is based on first principles of verifying Kullback-Leibler (KL) and testing conditions (see e.g., [11]). We first prove (6.3) and (6.5).

Part I: Kullback-Leibler conditions. Let $f\sim\mathcal{N}_{n}(\bm{X}\bm{\beta},\sigma^{2}\bm{I}_{n})$ and $f_{0}\sim\mathcal{N}_{n}(\bm{X}\bm{\beta}_{0},\sigma_{0}^{2}\bm{I}_{n})$, and let $\Pi(\cdot)$ denote the prior (6.2). We first show that for our choice of $\epsilon_{n}=\sqrt{s_{0}\log G/n}$,

[TABLE]

for some constant $C_{1}>0$, where $K(\cdot,\cdot)$ denotes the KL divergence and $V(\cdot,\cdot)$ denotes the KL variation. The KL divergence between $f_{0}$ and $f$ is

[TABLE]

and the KL variation between $f_{0}$ and $f$ is

[TABLE]

Define the two events $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ as follows:

[TABLE]

and

[TABLE]

Following from (D.9)-(D.13), we may write $\Pi(K(f_{0},f)\leq\epsilon_{n}^{2},V(f_{0},f)\leq\epsilon_{n}^{2})=\Pi(\mathcal{A}_{2}|\mathcal{A}_{1})\Pi(\mathcal{A}_{1})$. We derive lower bounds for $\Pi(\mathcal{A}_{1})$ and $\Pi(\mathcal{A}_{2}|\mathcal{A}_{1})$ separately. Noting that we may rewrite $\mathcal{A}_{1}$ as

[TABLE]

and expanding $\log(\sigma_{0}^{2}/\sigma^{2})$ in the powers of $1-\sigma_{0}^{2}/\sigma^{2}$ to get $\sigma_{0}^{2}/\sigma^{2}-1-\log(\sigma_{0}^{2}/\sigma^{2})\sim(1-\sigma_{0}^{2}/\sigma^{2})^{2}/2$, it is clear that $\mathcal{A}_{1}\supset\mathcal{A}_{1}^{\star}$, where $\mathcal{A}_{1}^{\star}=\{\sigma^{2}:\sigma_{0}^{2}/(\epsilon_{n}+1)\leq\sigma^{2}\leq\sigma_{0}^{2}\}$. Thus, since $\sigma^{2}\sim\mathcal{IG}(c_{0},d_{0})$, we have for sufficiently large $n$,

[TABLE]

Thus, from (D.3), we have

[TABLE]

since $n\epsilon_{n}^{2}\rightarrow\infty$. Next, we consider $\Pi(\mathcal{A}_{2}|\mathcal{A}_{1})$. We have

[TABLE]

and conditional on $\mathcal{A}_{1}$, we have that the previous display is bounded above by

[TABLE]

for large $n$ (since $\epsilon_{n}<1$ when $n$ is large). Since $\mathcal{A}_{1}\supset\mathcal{A}_{1}^{\star}$, where $\mathcal{A}_{1}^{\star}$ was defined earlier, the left-hand side of both expressions in (D.13) can be bounded above by a constant multiple of $\lVert\bm{X}(\bm{\beta}-\bm{\beta}_{0})\rVert_{2}^{2}$, conditional on $\mathcal{A}_{1}$. Therefore, for some constant $b_{1}>0$, $\Pi(\mathcal{A}_{2}|\mathcal{A}_{1})$ is bounded below by

[TABLE]

where we used Assumption (A3) in the second inequality, and in the third inequality, we used the fact that conditional on $\theta$, the SSGL prior is separable, so $\pi(\bm{\beta}|\theta)=\pi_{S_{0}}(\bm{\beta}|\theta)\pi_{S_{0}^{c}}(\bm{\beta}|\theta)$. We proceed to lower-bound each bracketed integrand term in (D.3) separately. Changing the variable $\bm{\beta}-\bm{\beta}_{0}\rightarrow\bm{b}$ and using the fact that $\pi_{S_{0}}(\bm{\beta}|\theta)>\theta^{s_{0}}\prod_{g\in S_{0}}\left[C_{g}\lambda_{1}^{m_{g}}e^{-\lambda_{1}\lVert\bm{\beta}_{g}\rVert_{2}}\right]$ and $\lVert\bm{z}\rVert_{2}\leq\lVert\bm{z}\rVert_{1}$ for any vector $\bm{z}$, we have as a lower bound for the first term in (D.3),

[TABLE]

Each of the integral terms in (D.17) is the probability of the first $m_{g}$ events of a Poisson process happening before time $b_{1}\epsilon_{n}/2s_{0}\sqrt{kn^{\alpha}}$. Using similar arguments as those in the proof of Lemma 5.1 of [29], we obtain as a lower bound for the product of integrals in (D.17),

[TABLE]

Combining (D.17)-(D.3), we have the following lower bound for the first bracketed term in (D.3):

[TABLE]

Now, noting that $\pi_{S_{0}^{c}}(\bm{\beta}|\theta)>(1-\theta)^{G-s_{0}}\prod_{g\in S_{0}^{c}}\left[C_{g}\lambda_{0}^{m_{g}}e^{-\lambda_{0}\lVert\bm{\beta}_{g}\rVert_{2}}\right]$, we further bound the second bracketed term in (D.3) from below. Let $\check{\pi}(\cdot)$ denote the density, $\check{\pi}(\bm{\beta}_{g})=C_{g}\lambda_{0}^{m_{g}}e^{-\lambda_{0}\lVert\bm{\beta}_{g}\rVert_{2}}$. We have

[TABLE]

where we used an application of the Markov inequality and Lemma 4 in the second line. Combining (D.19)-(D.3) gives as a lower-bound for (D.3),

[TABLE]

Let us consider the second bracketed term in (D.3) first. By assumption, $\lambda_{0}=(1-\theta)/\theta$. Further, $\lambda_{0}^{2}=(1-\theta)^{2}/\theta^{2}$ is monotonically decreasing in $\theta$ for $\theta\in(0,1)$. Hence, for constant $c>2$ in the $\mathcal{B}(1,G^{c})$ prior on $\theta$, a lower bound for the second bracketed term in (D.3) is

[TABLE]

where in the third line, we used our assumptions about the growth rates for $m_{\max}$, $G$, and $s_{0}$ in Assumptions (A1)-(A2) and the fact that $c>2$. In the fourth line, we used the fact that $(1-1/x)^{x}\rightarrow e^{-1}$ as $x\rightarrow\infty$ and $\theta\sim\mathcal{B}(1,G^{c})$. In the sixth line, we used the fact that the bracketed term in the fifth line can be bounded below by $e^{-1}-e^{-2}$ for sufficiently large $n$.

Combining (D.3)-(D.3), we obtain for sufficiently large $n$,

[TABLE]

We examine each of the terms in (D.3) separately. By Assumptions (A1) and (A5) and the fact that $\lambda_{1}\asymp 1/n$, we have

[TABLE]

and

[TABLE]

Next, using the facts that $x!\leq x^{x}$ for $x\in\mathbb{N}$ and Assumption (A1), we have

[TABLE]

Using the fact that the normalizing constant, $C_{g}=2^{-m_{g}}\pi^{-(m_{g}-1)/2}[\Gamma((m_{g}+1)/2)]^{-1}$, we also have

[TABLE]

where we used the fact that $\Gamma((m_{g}+1)/2)\leq\Gamma(m_{g}+1)=m_{g}!$. Finally, since $\lambda_{1}\asymp 1/n$ and using Assumption (A1) that $m_{\max}=O(\log G/\log n)$, we have

[TABLE]

Finally, it is clear by the definition of $n\epsilon_{n}^{2}$ and the fact that $c>2$ is a constant that

[TABLE]

Combining all of the above, together with (D.3), we have

[TABLE]

By (D.15) and (D.24), we may choose a large constant $C_{1}>0$, so that

[TABLE]

so the Kullback-Leibler condition (D.9) holds.

Part II: Testing conditions. To complete the proof, we show the existence of a sieve $\mathcal{F}_{n}$ such that

[TABLE]

for positive constant $C_{2}>C_{1}+2$, where $C_{1}$ is the constant from (D.9), and a sequence of test functions $\phi_{n}\in[0,1]$ such that

[TABLE]

and

[TABLE]

for some $C_{4}>0$, where $\nu$ is from Assumption (A4). Recall that $\omega_{g}\equiv\omega_{g}(\lambda_{0},\lambda_{1},\theta)=\frac{1}{\lambda_{0}-\lambda_{1}}\log\left[\frac{1-\theta}{\theta}\frac{\lambda_{0}^{m_{g}}}{\lambda_{1}^{m_{g}}}\right]$. Choose $C_{3}\geq C_{1}+2+\log 3$, and consider the sieve,

[TABLE]

where $c_{0}$ is from $\mathcal{IG}(c_{0},d_{0})$ prior on $\sigma^{2}$ and $\lvert\bm{\gamma}(\bm{\beta})\rvert$ denotes the generalized dimensionality (6.6).

We first verify (D.25). We have

[TABLE]

We focus on bounding each of the terms in (D.29) separately. First, let $\theta_{0}=C_{3}s_{0}\log G/G^{c}$, where $c>2$ is the constant in the $\mathcal{B}(1,G^{c})$ prior on $\theta$. Similarly as in the proof of Theorem 6.3 in Ročková and George [36], we have $\pi(\bm{\beta}_{g}|\theta)<2\theta C_{g}\lambda_{1}^{m_{g}}e^{-\lambda_{1}\lVert\bm{\beta}_{g}\rVert_{2}}$ for all $\lVert\bm{\beta}_{g}\rVert_{2}>\omega_{g}$. We have for any $\theta\leq\theta_{0}$ that

[TABLE]

where we used the assumption that $\lambda_{1}\asymp 1/n$, the definition of $\omega_{g}$, and the fact that $\theta\leq\theta_{0}$ to bound the first integral term from above by $\prod_{g\in S}(1/n)^{m_{g}}\leq n^{-|S|}$, and we bounded the second integral term above by one. We then have

[TABLE]

Note that since $s_{0}=o(n/\log G)$ by Assumption (A1), $G\gg n$, and $c>2$, we have $\theta_{0}\leq C_{3}n/G^{c}<1/G^{2}$ for sufficiently large $n$. Following from (D.3), we thus have that for sufficiently large $n$,

[TABLE]

where we used the inequality ${G\choose k}\leq(eG/k)^{k}$ in the third line of the display.

Next, since $\theta\sim\mathcal{B}(1,G^{c})$, we have

[TABLE]

Combining (D.3)-(D.3), we have that

[TABLE]

Finally, we have as a bound for the second term in (D.29),

[TABLE]

Combining (D.29)-(D.3), we have

[TABLE]

and so given our choice of $C_{3}$, (D.29) is asymptotically bounded from above by $\exp(-C_{2}n\epsilon_{n}^{2})$ for some $C_{2}\geq C_{1}+2$. This proves (D.25).

We now proceed to prove (D.26). Our proof is based on the technique used in Song and Liang [40] with suitable modifications. For $\xi\subset\{1,\ldots,G\}$, let $\bm{X}_{\xi}$ denote the submatrix of $\bm{X}$ with submatrices indexed by $\xi$, where $\lvert\xi\rvert\leq\bar{p}$ and $\bar{p}$ is from Assumption (A4). Let $\widehat{\bm{\beta}}_{\xi}=(\bm{X}_{\xi}^{T}\bm{X}_{\xi})^{-1}\bm{X}_{\xi}^{T}\bm{Y}$ and $\bm{\beta}_{0\xi}$ denote the subvector of $\bm{\beta}_{0}$ with groups indexed by $\xi$. Let $m_{\xi}=\sum_{g\in\xi}m_{g}$, and let $\widehat{\sigma}_{\xi}^{2}=\lVert\bm{Y}-\bm{X}_{\xi}\widehat{\bm{\beta}}_{\xi}\rVert_{2}^{2}/(n-m_{\xi})$. Note that $\widehat{\bm{\beta}}_{\xi}$ and $\widehat{\sigma}_{\xi}^{2}$ both exist and are unique because of Assumptions (A1), (A2), and (A4) (which combined, gives us that $m_{\xi}=o(n)$).

Let $\widetilde{p}$ be an integer satisfying $\widetilde{p}\asymp s_{0}$ and $\widetilde{p}\leq\bar{p}-s_{0}$, where $\bar{p}$ is from Assumption (A4), and the specific choice of $\widetilde{p}$ will be given below. Recall that $S_{0}$ is the set of true nonzero groups with cardinality $s_{0}=\lvert S_{0}\rvert$. Similar to [40], we consider the test function $\phi_{n}=\max\{\phi_{n}^{\prime},\tilde{\phi}_{n}\}$, where

[TABLE]

Because of Assumption (A4), we have $\widetilde{p}\prec n$ and $\widetilde{p}\prec n\epsilon_{n}^{2}$. Additionally, since $\epsilon_{n}=o(1)$, we can use almost identical arguments as those used to establish (A.5)-(A.6) in the proof of Theorem A.1 of [40] to show that for any $\xi$ satisfying $\xi\supset S_{0},|\xi|\leq\widetilde{p}$,

[TABLE]

for some constant $\hat{c}_{4}>0$, and for any $\xi$ satisfying $\xi\supset S_{0},|\xi|\leq\widetilde{p}$,

[TABLE]

for some $\tilde{c}_{4}>0$. Using the proof of Theorem A.1 in [40], we may then choose $\widetilde{p}=\lfloor\min\{c_{4}^{\prime},\tilde{c}_{4}\}n\epsilon_{n}^{2}/(2\log G)\rfloor$, and then

[TABLE]

for some $\check{c}_{4}>0$. Next, define the set,

[TABLE]

By Lemma 5, we have

[TABLE]

Similar to [40], we consider $\mathcal{C}\subset\widehat{\mathcal{C}}\cup\widetilde{\mathcal{C}}$, where

[TABLE]

and so an upper bound for (D.38) is

[TABLE]

Let $\tilde{\xi}=\{g:\lVert\bm{\beta}_{g}\rVert_{2}>\omega_{g}\}\cup S_{0}$, $m_{\tilde{\xi}}=\sum_{g\in\tilde{\xi}}m_{g}$, and $\tilde{\xi}^{c}=\{1,\ldots,G\}\setminus\tilde{\xi}$. For any $f\in\mathcal{F}_{n}$ such that $(\bm{\beta},\sigma^{2})\in\hat{\mathcal{C}}\cup\tilde{\mathcal{C}}$, we must have then that $\lvert\tilde{\xi}\rvert\leq C_{3}s_{0}+s_{0}\leq\bar{p}$, by Assumption (A4). By (D.3), the prior puts exponentially vanishing probability on values of $\theta>\theta_{0}$ where $\theta_{0}=C_{3}s_{0}\log G/G^{c}<1/(G^{2}+1)$ for large $G$. Since $\lambda_{0}=(1-\theta)/\theta$ is monotonic decreasing in $\theta$, we have that with probability greater than $1-e^{-C_{3}n\epsilon_{n}^{2}}$, $\lambda_{0}\geq G^{2}$. Combining this fact with Assumption (A3) and using $\mathcal{F}_{n}$ in (D.28), we have that for any $f\in\mathcal{F}_{n},(\bm{\beta},\sigma^{2})\in\hat{C}\cup\tilde{C}$ and sufficiently large $n$,

[TABLE]

where $\nu$ is from Assumption (A4). In the above display, we used Lemma 6 in the second inequality, while the last inequality follows from our assumptions on $(\theta,\lambda_{0},\lambda_{1})$ and $m_{\max}$, so one can show that the bracketed term in the third line is asymptotically bounded above by $D\sqrt{\nu_{1}}\sqrt{n^{1-\alpha}}\sigma_{0}\epsilon_{n}$ for large $n$ and any constant $D>0$. Thus, using nearly identical arguments as those used to prove Part I of Theorem A.1 in [40], we have

[TABLE]

where the noncentrality parameter $\zeta$ satisfies $\zeta\leq n\epsilon_{n}^{2}\nu_{1}\sigma_{0}^{2}/16\sigma^{2}$, and the last inequality follows from the fact that the noncentral $\chi^{2}$ distribution is subexponential and Bernstein’s inequality (see Lemmas A.1 and A.2 in [40]).

Using the arguments in Part I of the proof of Theorem A.1 in [40], we also have that for large $n$,

[TABLE]

where the second inequality in the above display holds since $\lVert\bm{\beta}_{\tilde{\xi}}-\bm{\beta}_{0\tilde{\xi}}\rVert_{2}\geq\lVert\bm{\beta}-\bm{\beta}_{0}\rVert_{2}-\lVert\bm{\beta}_{\tilde{\xi}^{c}}\rVert_{2}$, and since (D.3) can be further bounded from above by $\sqrt{kn^{\alpha}}\sqrt{\nu_{1}}\sigma_{0}\epsilon_{n}$ and thus $\lVert\bm{\beta}_{\tilde{\xi}^{c}}\rVert\leq\sqrt{\nu_{1}}\sigma_{0}\epsilon_{n}$. Therefore, we have for $f\in\mathcal{F}_{n},(\bm{\beta},\sigma^{2})\in\tilde{C}$,

[TABLE]

while by Assumption (A4) and (D.3), we also have

[TABLE]

and then we used the fact that on the set $\tilde{C}$, $\sigma=\sigma_{0}$. The last three inequalities in (D.3) follow from Assumption (A4), the fact that $\lvert\widetilde{\xi}\rvert\leq\bar{p}\prec n\epsilon_{n}^{2}$, and the fact that for all $m>0$, $\Pr(\chi^{2}_{m}\geq x)\leq\exp(-x/4)$ whenever $x\geq 8m$. Altogether, combining (D.38)-(D.3), we have that

[TABLE]

where $\hat{c}_{4}>0$ and $\tilde{c}_{4}>0$ are the constants from (D.3) and (D.3).

Now set $C_{4}=\min\{\hat{c}_{4},\tilde{c}_{4},\check{c}_{4}\}$, where $\check{c}_{4}$ is the constant from (D.37). By and (D.37) and (D.43), this choice of $C_{4}$ will satisfy both testing conditions (D.26) and (D.27).

Since we have verified (D.9) and (D.25)-(D.27) for $\epsilon_{n}=\sqrt{s_{0}\log G/n}$, we have

[TABLE]

and

[TABLE]

i.e. we have proven (6.3) and (6.5).

Part III. Posterior contraction under prediction error loss. The proof is very similar to the proof of (6.3). The only difference is the testing conditions. We use the same sieve $\mathcal{F}_{n}$ as that in (D.28) so that (D.25) holds, but now, we need to show the existence of a different sequence of test functions $\tau_{n}\in[0,1]$ such that

[TABLE]

and

[TABLE]

Let $\widetilde{p}$ be the same integer from (D.36) and consider the test function $\tau_{n}=\max\{\tau_{n}^{\prime},\tilde{\tau}_{n}\}$, where

[TABLE]

Using Assumption (A4) that for any $\xi\subset\{1,\ldots,G\}$ such that $\lvert\xi\rvert\leq\bar{p}$, $\lambda_{\max}(\bm{X}_{\xi}^{T}\bm{X}_{\xi})\leq n\nu_{2}$ for some $\nu_{2}>0$, we have that

[TABLE]

and so

[TABLE]

Therefore, using similar steps as those in Part II of the proof, we can show that our chosen sequence of tests $\tau_{n}$ satisfies (D.44) and (D.45). We thus arrive at

[TABLE]

i.e. we have proven (6.4). ∎

Proof of Theorem 3.

According to Part I of the proof of Theorem 2, we have that for $\epsilon_{n}=\sqrt{s_{0}\log G/n}$,

[TABLE]

for some $C>0$. Thus, by Lemma 8.10 of [12], there exist positive constants $C_{1}$ and $C_{2}$ such that the event,

[TABLE]

satisfies

[TABLE]

Define the set $\mathcal{T}=\{\bm{\beta}:\lvert\bm{\gamma}(\bm{\beta})\rvert\leq C_{3}s_{0}\}$, where we choose $C_{3}>1+C_{2}$. We must show that $\mathbb{E}_{0}\Pi(\mathcal{T}^{c}|\bm{Y})\rightarrow 0$ as $n\rightarrow\infty.$ The posterior probability $\Pi(\mathcal{T}^{c}|\bm{Y})$ is given by

[TABLE]

By (D.48), the denominator of (D.49) is bounded below by $e^{-(1+C_{2})n\epsilon_{n}^{2}}$. For the numerator of (D.49), we have as an upper bound,

[TABLE]

Using the same arguments as (D.3)-(D.34) in the proof of Theorem 2, we can show that

[TABLE]

Combining (D.47)-(D.50), we have that

[TABLE]

since $C_{3}>1+C_{2}$. This proves (6.8). ∎

Proof of Theorem 4.

Let $f_{0j}(\bm{X}_{j})$ be an $n\times 1$ vector with $i$th entry equal to $f_{0j}(X_{ij})$. Note that proving posterior contraction with respect to the empirical norm (6.11) is equivalent to proving that

[TABLE]

so to prove the theorem, it suffices to prove (D.52). Let $f\sim\mathcal{N}_{n}(\widetilde{\bm{X}}\bm{\beta},\sigma^{2}\bm{I}_{n})$ and $f_{0}\sim\mathcal{N}_{n}(\widetilde{\bm{X}}\bm{\beta}_{0}+\bm{\delta},\sigma_{0}^{2}\bm{I}_{n})$, and let $\Pi(\cdot)$ denote the prior (6.2). Similar to the proof for Theorem 2, we show that for our choice of $\epsilon_{n}^{2}=s_{0}\log p/n+s_{0}n^{-2\kappa/(2\kappa+1)}$ and some constant $C_{1}>0$,

[TABLE]

and the existence of a sieve $\mathcal{F}_{n}$ such that

[TABLE]

for positive constant $C_{2}>C_{1}+2$, and a sequence of test functions $\phi_{n}\in[0,1]$ such that

[TABLE]

and

[TABLE]

for some $C_{4}>0$ and $\tilde{c}_{0}>0$.

We first verify (D.53). The KL divergence between $f_{0}$ and $f$ is

[TABLE]

and the KL variation between $f_{0}$ and $f$ is

[TABLE]

Define the two events $\widetilde{\mathcal{A}}_{1}$ and $\widetilde{\mathcal{A}}_{2}$ as follows:

[TABLE]

and

[TABLE]

Following from (D.57)-(D.60), we have $\Pi(K(f_{0},f)\leq n\epsilon_{n}^{2},V(f_{0},f)\leq n\epsilon_{n}^{2})=\Pi(\widetilde{\mathcal{A}}_{2}|\widetilde{\mathcal{A}}_{1})\Pi(\widetilde{\mathcal{A}}_{1})$. Using the steps we used to prove (D.15) in part I of the proof of Theorem 2, we have

[TABLE]

for some sufficiently large $C_{1}>0$. Following similar reasoning as in the proof of Theorem 2, we also have for some $b_{2}>0$,

[TABLE]

Using Assumptions (B3) and (B6), we then have

[TABLE]

and so (D.62) can be asymptotically lower bounded by

[TABLE]

where we used Assumption (B1) that $d\asymp n^{1/(2\kappa+1)}$. Using very similar arguments as those used to prove (D.24), this term can also be lower bounded by $\exp(-C_{1}n\epsilon_{n}^{2}/2)$. Altogether, we have

[TABLE]

Combining (D.61) and (D.63), we have that (D.53) holds. To verify (D.54), we choose $C_{3}\geq C_{1}+2+\log 3$ and use the same sieve $\mathcal{F}_{n}$ as the one we employed in the proof of Theorem 2 (eq. (D.28)), and then (D.54) holds for our choice of $\mathcal{F}_{n}$.

Finally, we follow the recipe of Wei et al. [48] and Song and Liang [40] to construct our test function $\phi_{n}$ which will satisfy both (D.55) and (D.56). For $\xi\subset\{1,\ldots,p\}$, let $\widetilde{\bm{X}}_{\xi}$ denote the submatrix of $\widetilde{\bm{X}}$ with submatrices indexed by $\xi$, where $\lvert\xi\rvert\leq\bar{p}$ and $\bar{p}$ is from Assumption (B4). Let $\widehat{\bm{\beta}}_{\xi}=(\widetilde{\bm{X}}_{\xi}^{T}\widetilde{\bm{X}}_{\xi})^{-1}\widetilde{\bm{X}}_{\xi}^{T}\bm{Y}$ and $\bm{\beta}_{0\xi}$ denote the subvector of $\bm{\beta}_{0}$ with basis coefficients appearing in $\xi$. Then the total number of elements in $\widehat{\bm{\beta}}_{\xi}$ is $d\lvert\xi\rvert$. Finally, let $\widehat{\sigma}_{\xi}^{2}=\bm{Y}^{T}(\bm{I}_{n}-\bm{H}_{\xi})\bm{Y}/(n-d\lvert\xi\rvert)$, where $\bm{H}_{\xi}=\widetilde{\bm{X}}_{\xi}(\widetilde{\bm{X}}_{\xi}^{T}\widetilde{\bm{X}}_{\xi})^{-1}\widetilde{\bm{X}}_{\xi}^{T}$ is the hat matrix for the subgroup $\xi$.

Let $\widetilde{p}$ be an integer satisfying $\widetilde{p}\asymp s_{0}$ and $\widetilde{p}\leq\bar{p}-s_{0}$, where $\bar{p}$ is from Assumption (B4) and the specific choice for $\widetilde{p}$ will be given later. Recall that $S_{0}$ is the set of true nonzero groups with cardinality $s_{0}=\lvert S_{0}\rvert$. Similar to [48], we consider the test function, $\phi_{n}=\max\{\phi_{n}^{\prime},\tilde{\phi}_{n}\}$, where

[TABLE]

for some positive constants $c_{0}^{\prime}$ and $\tilde{c}_{0}$. Using Assumptions (B1) and (B4), we have that for any $\xi$ in our test $\phi_{n}$, $d\lvert\xi\rvert\leq d(\widetilde{p}+s_{0})\leq d\bar{p}\prec n\epsilon_{n}^{2}$. Using essentially the same arguments as those in the proof for Theorem 4.1 in [48], we have that for any $\xi$ which satisfies $\xi\supset S_{0}$ so that $\lvert\xi\rvert\leq\widetilde{p}+s_{0}$,

[TABLE]

for some $c_{0}^{\prime\prime}>0$. By Assumption (B3), we also have

[TABLE]

and using the fact that $\lVert\bm{\delta}\rVert_{2}\lesssim\sqrt{ns_{0}}d^{-\kappa}\lesssim\tilde{c}_{0}\sigma_{0}\sqrt{n}\epsilon_{n}/2$ (by Assumptions (B1) and (B6)), we have that for any $\xi$ such that $\xi\supset S_{0},\lvert\xi\rvert\leq\widetilde{p}+s_{0}$,

[TABLE]

for some $\tilde{c}_{4}>0$, where we used the proof of Theorem A.1 in [40] to arrive at the final inequality. Again, as in the proof of Theorem A.1 of [40], we choose $\widetilde{p}=\lfloor\min\{c_{4}^{\prime},\tilde{c}_{4}\}n\epsilon_{n}^{2}/(2\log p)\rfloor$, and then

[TABLE]

for some $\check{c}_{4}>0$. Next, we define the set,

[TABLE]

By Lemma 5, we have

[TABLE]

Similar to [40], we consider $\mathcal{C}\subset\widehat{\mathcal{C}}\cup\widetilde{\mathcal{C}}$, where

[TABLE]

and so an upper bound for (D.3) is

[TABLE]

Using very similar arguments as those used to prove (D.43) in Theorem 2 and using Assumptions (B1) and (B6), so that the bias $\lVert\bm{\delta}\rVert_{2}^{2}\lesssim ns_{0}d^{-2\kappa}\lesssim n\epsilon_{n}^{2}$, we can show that (D.3) can be further bounded from above as

[TABLE]

where $\hat{c}_{4}>0$ and $\tilde{c}_{4}>0$ are the constants from (D.65) and (D.66).

Choose $C_{4}=\min\{\check{c}_{4},\hat{c}_{4},\tilde{c}_{4}\}$, and we have from (D.66) and (D.3) that (D.55) and (D.56) both hold.

Since we have verified (D.53) and (D.54)-(D.56) for our choice of $\epsilon_{n}^{2}=s_{0}\log p/n+s_{0}n^{-2\kappa/(2\kappa+1)}$, it follows that

[TABLE]

and

[TABLE]

i.e. we have proven (D.52), or equivalently, (6.11) and (6.12).

∎

Proof of Theorem 5.

The proof is very similar to the proof of Theorem 3 and is thus omitted. ∎