# Informative co-data learning for high-dimensional Horseshoe regression

Claudio Busatto  Department of Statistics, Computer Science, Applications “G.
Parenti”, University of Florence, Florence, Italy  Mark van de Wiel  Dep.
Epidemiology & Data Science, Amsterdam Public Health research institute,
Amsterdam University Medical Centers, Amsterdam, The Netherlands

###### Abstract

High-dimensional data often arise from clinical genomics research to infer
relevant predictors of a particular trait. A way to improve the predictive
performance is to include information on the predictors derived from prior
knowledge or previous studies. Such information is also referred to as “co-
data”. To this aim, we develop a novel Bayesian model for including co-data in
a high-dimensional regression framework, called Informative Horseshoe
regression (infHS). The proposed approach regresses the prior variances of the
regression parameters on the co-data variables, improving variable selection
and prediction. We implement both a Gibbs sampler and a Variational
approximation algorithm. The former is suited for applications of moderate
dimensions which, besides prediction, target posterior inference, whereas the
computational efficiency of the latter allows handling a very large number of
variables. We show the benefits from including co-data with a simulation
study. Eventually, we demonstrate that infHS outperforms competing approaches
for two genomics applications.

Keywords: Bayesian inference, Variational Bayes, informative shrinkage prior,
Horseshoe prior, co-data information.

##  1 Introduction

The analysis of high-dimensional data sets is a main interest in many
scientific fields. In particular, clinical research often deals with huge
amount of data, such as genes expression or genome-wide methylation levels for
relatively few samples, due to budget or practical constraints. We consider
regression models to study a clinical outcome with this type of data. Since
the number of parameters, $p$, overwhelms sample size, $n$, we develop an
approach to improve the overall performance of the model by incorporating
(prior) external knowledge in the estimating process. Such an external source
of information is referred to as co-data (complementary data; Neuenschwander
et al.,, 2016), as it provides additional information about the covariates. We
consider two different types: continuous, such as $p$-values from previous
studies, or categorical, such as membership to a group, e.g. chromosome or
particular genomic regions.

Several methods allow to incorporate one source of auxiliary information in a
regression framework (Tai and Pan,, 2007; Boonstra et al.,, 2013). A popular
method is the (sparse) group LASSO (Yuan and Lin,, 2006; Simon et al.,, 2013),
which penalizes groups of variables using one common hyperparameter for all
groups. Such a solution is attractive when the number of covariate groups is
large, but lacks flexibility and fails to adapt locally in other settings,
leading to sub-optimal results (Münch et al.,, 2019). More recent work focuses
on the estimation of adaptive penalties with prior variances specific for each
group (Van de Wiel et al.,, 2019; Velten and Huber,, 2019, Münch et al.,
(2019)). These methods, however, are restrictive in use as they either deal
with just one (discrete) co-data source or handle one specific type of outcome
only. In Van Nee et al., (2021), instead, a Ridge regression method is
proposed that allows for multiple co-data sources by regressing the local
variances on the co-data. With their approach, the co-data regression
parameters are estimated independently for each co-data source and are
eventually combined using a vector of weights, where each weight is related to
the importance of a single co-data source.

Here, we present a novel Bayesian method for both linear and binary regression
that accounts for multiple co-data information. In particular, we introduce a
generalization of the Horseshoe regression (Carvalho et al.,, 2010), referred
to as the informative Horseshoe regression model (infHS). We regress the local
variances on the co-data variables, following the work of Van Nee et al.,
(2021). In contrast to Kpogbezan et al., (2019), where the Horseshoe prior is
used with only a single two-group co-data source, our model is flexible with
respect to the co-data type, as it allows for both continuous and discrete co-
data predictors. Moreover, it extends to binary outcome via probit regression
(Albert and Chib,, 1993). Unlike Van Nee et al., (2021), it tackles the sparse
setting and co-data regression parameters related to different sources are
estimated jointly, avoiding multiple regressions for each co-data source
separately.

We first propose a Gibbs sampler for iteratively updating posterior
parameters. We introduce a novel rejection sampling method for sampling from
the non-analytical full-conditional distribution of the local variances. When
the number of variables increases, we rely on the computational methods
presented in Bhattacharya et al., (2016) to sample efficiently from a
multivariate normal density. To make the method applicable to particularly
large $p$ settings, we develop a Variational Bayes approximation to the joint
posterior distribution, using techniques from Münch et al., (2019) to
efficiently optimize the target density of the variational distribution,
rendering an algorithm with computation time linear in $p$. With simulations
and two data applications, we show that both prediction and variable selection
benefit from the inclusion of co-data information, the latter being
particularly relevant under the Horseshoe setting.

The paper in organized as follows. Section 2 introduces the hierarchical
structure of the model and discusses the parametrizations. In Sections 3 and 4
we develop the Gibbs sampling algorithm and the Variational Bayes
approximation to the joint posterior distribution, respectively. Section 5
illustrates the benefit of co-data information on variable selection with a
simulation study, whereas Section 6 presents applications of our model to two
data sets, one from genetics and one from cancer genomics. We conclude with
discussions and possible extensions in Section 7.

##  2 The model

Let $\mathbf{y}\in\mathbb{R}^{n}$ be the response vector and
$\mathbf{X}\in\mathbb{R}^{n\times(p+1)}$ the design matrix, with the first
column of ones. We regress $\mathbf{y}$ on $\mathbf{X}$ using a generalized
linear model (GLM) with regression coefficient vector
$\boldsymbol{\beta}=\left[\beta_{0}\;\beta_{1}\dots\beta_{p}\right]^{\intercal}$.

$\displaystyle
Y_{i}\mid\mathbf{x}_{i},\boldsymbol{\beta}\;\stackrel{{\scriptstyle\mathclap{\mbox{{ind}}}}}{{\sim}}$
$\displaystyle\;p\left(Y_{i}\mid\mathbf{x}_{i},\boldsymbol{\beta}\right)$
$\displaystyle
E_{Y_{i}\mid\mathbf{x}_{i},\boldsymbol{\beta}}\left(Y_{i}\right)=$
$\displaystyle\;g^{-1}\left(\mathbf{x}_{i}^{\intercal}\boldsymbol{\beta}\right),\quad
i=1,\dots,n,$

where $\mathbf{x}_{i}=\left[1\;x_{i1}\dots x_{ip}\right]^{\intercal}$ is the
set of covariates related to observation $i$.

Following Carvalho et al., (2010), the Horseshoe prior locally shrinks
regression parameters toward zero and provides a sparse solution for
$\boldsymbol{\beta}$. We assume a normal prior distribution for each
$\beta_{j}$, where the variance is decomposed in a global scale parameter
$\tau$ and a local shrinkage parameter $\lambda_{j}$. Formally,

$\displaystyle\beta_{0}\mid\sigma^{2},\tau,\lambda_{0}\sim$
$\displaystyle\;\mathcal{N}\left(0,\sigma^{2}\tau^{2}\lambda_{0}^{2}\right),$
$\displaystyle\lambda_{0}\sim$ $\displaystyle\;\mathcal{C}^{+}(0,1),$
$\displaystyle\beta_{j}\mid\sigma^{2},\tau,\lambda_{j}\sim$
$\displaystyle\;\mathcal{N}\left(0,\sigma^{2}\tau^{2}\lambda_{j}^{2}\right),\quad
j=1,\dots,p.$

Suppose that $D$ different co-data sources
$\mathbf{Z}_{d}\in\mathbb{R}^{p\times m_{d}}$ are available, where
$\sum_{d}m_{d}=M$ and $d=1,\dots,D$. In order to capture the external
information effect, Van Nee et al., (2021) introduce parameters $\omega_{d}>0$
and $\boldsymbol{\gamma}_{d}\in\mathbb{R}^{m_{d}}$, $d=1,\dots,D$, to model
covariate-specific shrinkage $\lambda_{j}$, $j=1,\dots,p$, as a function of
$\mathbf{Z}=\left[\mathbf{Z}_{1},\dots,\mathbf{Z}_{D}\right]$. Parameters
$\boldsymbol{\gamma}_{d}$ represent the regression coefficient vector related
to matrix $\mathbf{Z}_{d}$ and are estimated separately for each group $d$,
whereas co-data weights $\omega_{d}>0$ model the relative importance of group
$d$ and are introduced to combine the different co-data sources. Here,
parameter
$\boldsymbol{\gamma}=\left[\boldsymbol{\gamma}_{1}^{\intercal}\dots\boldsymbol{\gamma}_{D}^{\intercal}\right]^{\intercal}$
is update jointly by sampling from its full-conditional distribution. This way
the co-data sources are naturally combined and grouping weights $\omega_{d}$
can be excluded from the model. Therefore, the hierarchical set of prior
distributions in our model is

$\displaystyle\lambda_{j}\mid\mathbf{Z},\boldsymbol{\gamma}\sim$
$\displaystyle\;\mathcal{C}\left(\sum_{d=1}^{D}\left(\mathbf{z}_{j}^{d}\right)^{\intercal}\boldsymbol{\gamma}_{d},s_{0}^{2}\right)\cdot\mathbb{I}_{\left(\lambda_{j}>0\right)},\quad
j=1,\dots,p,$ (1)
$\displaystyle\boldsymbol{\gamma}_{d}\mid\boldsymbol{\Sigma}_{\gamma_{d}}\sim$
$\displaystyle\;\mathcal{N}_{m_{d}}\left(\mathbf{0},\boldsymbol{\Sigma}_{\gamma_{d}}\right),\quad
d=1,\dots,D,$ $\displaystyle\tau\sim$
$\displaystyle\;\mathcal{C}^{+}\left(0,1\right),$
$\displaystyle\sigma^{2}\sim$ $\displaystyle\;\mathcal{IG}\left(v,q\right),$

where $\mathcal{C}$ denotes the Cauchy distribution and $\mathcal{C}^{+}$ its
half-positive part. Here we model the location parameter of the prior local
variances mainly because it allows to sample co-data coefficients
$\boldsymbol{\gamma}$ jointly from a multivariate normal distribution (See
Section 3). Note that when the co-data are not informative and
$\boldsymbol{\gamma}_{d}\to 0$ for each $d$, the model reduces to the ordinary
Horseshoe prior from Carvalho et al., (2010) with $s_{0}^{2}=1$.

In Van de Wiel et al., (2019) the authors estimate prior covariance matrices
$\boldsymbol{\Sigma}_{d}$ from the data with an Empirical Bayes estimator
separately for each source. This approach, however, is computationally
burdensome, as it implies the implementation of multiple MCMC chains until
convergence. Here we consider a group-specific scale parameter
$\kappa_{d}^{2}$ and a ridge-like prior
$\boldsymbol{\Sigma}_{d}=\kappa_{d}^{2}\mathbf{I}_{m_{d}}$. A main advantage
of this approach is the computational efficiency, because only prior scale
parameters $\kappa_{d}^{2}$ have to be updated. A conjugated prior
distribution for $\kappa_{d}^{2}$ is

$\displaystyle\kappa_{d}^{2}\sim$
$\displaystyle\;\mathcal{IG}\left(a_{d},b_{d}\right).$ (2)

Parameters $\kappa_{d}^{2}$ act deep in the model and one can argue that they
should have a small impact on the global estimation process. For this reason,
a non-informative choice for $a_{d}$ and $b_{d}$ should suffice.

##  3 Posterior inference

In this section we introduce a Gibbs sampler that iteratively updates the
parameters by sampling from their full-conditional distributions. We show the
details of the algorithm for the linear regression model, under the assumption

$y_{i}=\mathbf{x}_{i}^{\intercal}\boldsymbol{\beta}+\varepsilon_{i},\quad\varepsilon_{i}\sim\mathcal{N}\left(0,\sigma^{2}\right),\quad
i=1,\dots,n.$

However, Pólya-Gamma latent variables (Polson et al.,, 2013) or probit GLM
(Albert and Chib,, 1993) can be introduced to augment the model and reach a
gaussian full-conditional distribution for $\boldsymbol{\beta}$ in a binary
regression model.

1\. Sampling $\boldsymbol{\beta}$ and $\sigma^{2}$. The full-conditional
distributions of $\boldsymbol{\beta}$ and $\sigma^{2}$ are

$\displaystyle\boldsymbol{\beta}\mid\mathbf{y},\mathbf{X},\sigma^{2},\tau^{2},\lambda_{0}^{2},\boldsymbol{\lambda}\sim$
$\displaystyle\;\mathcal{N}_{p+1}\left(\boldsymbol{\Sigma}_{\beta}^{\star}\mathbf{X}^{\intercal}\mathbf{y},\sigma^{2}\boldsymbol{\Sigma}_{\beta}^{\star}\right),\quad\boldsymbol{\Sigma}_{\beta}^{\star}=\left(\mathbf{X}^{\intercal}\mathbf{X}+\tau^{-2}\boldsymbol{\Lambda}^{-2}\right)^{-1},$
(3)
$\displaystyle\sigma^{2}\mid\mathbf{y},\mathbf{X},\boldsymbol{\beta},\tau^{2},\lambda_{0}^{2},\boldsymbol{\lambda}\sim$
$\displaystyle\;\mathcal{IG}\left(v+\frac{n+p+1}{2},q+\frac{1}{2}\|\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\|_{2}^{2}+\frac{\beta_{0}^{2}}{2\tau^{2}\lambda_{0}^{2}}+\frac{1}{2\tau^{2}}\sum_{j=1}^{p}\frac{\beta_{j}^{2}}{\lambda_{j}^{2}}\right),$

where
$\boldsymbol{\Lambda}=\text{diag}\left(\lambda_{0},\lambda_{1},\dots,\lambda_{p}\right)$.
Sampling from (3) requires the inversion of the $p\times p$ covariance matrix
$\boldsymbol{\Sigma}_{\beta}^{\star}$, which becomes computationally
infeasible when the number of covariates $p$ increases, as naive inversion is
of order $\mathcal{O}\left(p^{3}\right)$. For high-dimensional problems we
rely on the strategies introduced in Bhattacharya et al., (2016) and Johndrow
et al., (2020) to reduce the computational costs of sampling regression
parameters $\boldsymbol{\beta}$ to $\mathcal{O}\left(n^{2}p\right)$
operations.

2\. Sampling $\boldsymbol{\lambda}$, $\boldsymbol{\gamma}$ and
$\boldsymbol{\kappa}^{2}$. The prior distributions for $\lambda_{j}$, $j\neq
0$, and $\boldsymbol{\gamma}$ in (1) are not conjugated. To this aim, we rely
on the data-augmentation step proposed in Geweke, (1993) to reach a conjugated
framework and jointly update parameter $\boldsymbol{\gamma}$ by sampling from
a multivariate normal distribution.

######  Proposition 1.

Let $x$ and $y$ be random variables such that

$x\sim\mathcal{C}\left(m,s^{2}\right)\quad\text{and}\quad
y\sim\mathsf{t}_{v}\left(m,s^{2}\right),$

where $t_{v}$ denotes the Student-$t$ distribution. If $v=1$, then $X\sim Y$.

######  Proposition 2.

Let $x$ and $\theta$ be random variables such that

$x\mid
m,s^{2},\varphi^{2}\sim\mathcal{N}\left(m,s^{2}\varphi^{2}\right)\quad\text{and}\quad\varphi^{2}\sim\mathcal{IG}\left(\frac{v}{2},\frac{v}{2}\right),$

then $x\mid m,s^{2},v\sim\mathsf{t}_{v}\left(m,s^{2}\right)$ (Geweke,, 1993).

Following Propositions 1 and 2, after the introduction of a latent factor
$\varphi_{j}^{2}\sim\mathcal{IG}(1/2,1/2)$, the prior distribution of
$\lambda_{j}$ can be conveniently re-written as a normal distribution
truncated at $0$. This way the normal prior for co-data regression
coefficients $\boldsymbol{\gamma}$ is conjugated and the parameters can easily
be updated by sampling from a multivariate normal distribution. The set of
prior distributions for parameters $\lambda_{j}$, $\varphi_{j}^{2}$,
$\boldsymbol{\gamma}_{d}$ can be re-written as

$\displaystyle\lambda_{j}\mid\mathbf{Z},\boldsymbol{\gamma},\varphi_{j}^{2}\sim$
$\displaystyle\;\mathcal{N}\left(\sum_{d=1}^{D}\left(\mathbf{z}_{j}^{d}\right)^{\intercal}\boldsymbol{\gamma}_{d},s_{0}^{2}\varphi_{j}^{2}\right)\cdot\mathbb{I}_{\left(\lambda_{j}>0\right)},$
$\displaystyle\varphi_{j}^{2}\sim$
$\displaystyle\;\mathcal{IG}\left(\frac{1}{2},\frac{1}{2}\right),\quad
j=1,\dots,p,$ $\displaystyle\boldsymbol{\gamma}_{d}\mid\kappa_{d}^{2}\sim$
$\displaystyle\;\mathcal{N}_{m_{d}}\left(\mathbf{0},\kappa_{d}^{2}\mathbf{I}_{m_{d}}\right),$
$\displaystyle\kappa_{d}^{2}\sim$
$\displaystyle\;\mathcal{IG}\left(a_{d},b_{d}\right),\quad d=1,\dots,D.$

Let
$\mu_{j}=\sum_{d=1}^{D}\left(\mathbf{z}_{j}^{d}\right)^{\intercal}\boldsymbol{\gamma}_{d}$,
$\boldsymbol{\lambda}=\left[\lambda_{1}\dots\lambda_{p}\right]^{\intercal}$
and $\boldsymbol{\Phi}^{2}=\text{diag}\left(\boldsymbol{\varphi}\right)$, with
$\boldsymbol{\varphi}=\left[\varphi_{1}^{2},\dots,\varphi_{p}^{2}\right]^{\intercal}$.
The full-conditional distributions of $\lambda_{j}$, $\boldsymbol{\gamma}$,
$\varphi_{j}^{2}$ and $\kappa_{d}^{2}$ in the augmented model are

$\displaystyle\pi\left(\lambda_{j}\mid\mathbf{Z},\beta_{j},\sigma^{2},\tau^{2},\boldsymbol{\gamma},\varphi_{j}^{2}\right)\propto$
$\displaystyle\;\lambda_{j}^{-1}e^{-\frac{\beta_{j}^{2}}{2\sigma^{2}\tau^{2}\lambda_{j}^{2}}-\frac{\lambda_{j}^{2}}{2s_{0}^{2}\varphi_{j}^{2}}+\frac{\mu_{j}\lambda_{j}}{s_{0}^{2}\varphi_{j}^{2}}}\cdot\mathbb{I}_{\left(\lambda_{j}>0\right)},$
(4)
$\displaystyle\varphi_{j}^{2}\mid\mathbf{Z},\lambda_{j},\boldsymbol{\gamma}\sim$
$\displaystyle\;\mathcal{IG}\left(1,\frac{1}{2}+\frac{(\lambda_{j}-\mu_{j})^{2}}{2s_{0}^{2}}\right),$
$\displaystyle\boldsymbol{\gamma}\mid\mathbf{Z},\boldsymbol{\lambda},\boldsymbol{\varphi},\boldsymbol{\kappa}^{2}\sim$
$\displaystyle\;\mathcal{N}_{M}\left(\boldsymbol{\Sigma}_{\gamma}^{\star}\left(\mathbf{Z}^{\intercal}\boldsymbol{\Phi}^{-2}\boldsymbol{\lambda}\right),s_{0}^{2}\boldsymbol{\Sigma}_{\gamma}^{\star}\right),$
$\displaystyle\kappa_{d}^{2}\mid\boldsymbol{\gamma}_{d}\sim$
$\displaystyle\;\mathcal{IG}\left(a_{d}+\frac{m_{d}}{2},b_{d}+\frac{\boldsymbol{\gamma}_{d}^{\intercal}\boldsymbol{\gamma}_{d}}{2}\right),$

where
$\boldsymbol{\Sigma}_{\gamma}^{\star}=\left(\mathbf{Z}^{\intercal}\boldsymbol{\Phi}^{-2}\mathbf{Z}+s_{0}^{2}\mathbf{D}_{\kappa}^{-1}\right)^{-1}$
and
$\mathbf{D}_{\kappa}=\text{diag}\left(\kappa_{1}^{2}\mathbf{1}_{m_{1}},\dots,\kappa_{D}^{2}\mathbf{1}_{m_{D}}\right)$.
Details for sampling parameters $\lambda_{j}$ without computing the unknown
normalizing constant are given in Section B.1 of Supporting Information.

The introduced framework presents a computational advantage: the local
variances $\lambda_{j}^{2}$ can be computed in parallel, potentially improving
the efficiency of the model.

3\. Sampling $\tau^{2}$. The half-Cauchy prior for the global scale parameter
$\tau$ is not conjugated to the prior variance of $\beta_{j}$, $j=0,\dots,p$,
in a linear regression model with gaussian errors. Therefore, we rely on the
strategy proposed by Makalic and Schmidt, (2016) in order to update $\tau$.
The authors point out that the half-Cauchy distribution can be written as a
scale mixture of inverse-Gamma distributions, which allows conjugate updates
of $\tau^{2}$.

######  Proposition 3.

If $\tau^{2}\mid\zeta\sim\mathcal{IG}\left(\frac{1}{2},\frac{1}{\zeta}\right)$
and $\zeta\sim\mathcal{IG}\left(\frac{1}{2},1\right)$, then
$\tau\sim\mathcal{C}^{+}(0,1)$.

Following proposition 3 the full-conditional distributions for $\tau^{2}$ and
$\zeta$ are available in closed-form. In particular the parameters are updated
by sampling from the following densities:

$\displaystyle\tau^{2}\mid\boldsymbol{\beta},\sigma^{2},\lambda_{0}^{2},\boldsymbol{\lambda},\zeta\sim$
$\displaystyle\;\mathcal{IG}\left(\frac{p}{2}+1,\frac{1}{\zeta}+\frac{\beta_{0}^{2}}{2\sigma^{2}\lambda_{0}^{2}}+\frac{1}{2\sigma^{2}}\sum_{j=1}^{p}\frac{\beta_{j}^{2}}{\lambda_{j}^{2}}\right),$
$\displaystyle\zeta\mid\tau^{2}\sim$
$\displaystyle\;\mathcal{IG}\left(1,1+\frac{1}{\tau^{2}}\right).$

Section A.1 of Supporting Information discusses the update of parameter
$\lambda_{0}^{2}$, the local scale parameter of the intercept $\beta_{0}$,
using a similar proposition as 3.

##  4 Variational Bayes approximation

When the number of covariates is huge the method introduced in Section 2
becomes computationally infeasible. In this section an efficient approximation
of the joint posterior distribution is discussed. We refer to Section C.1 of
Supporting Information for the derivation of the Variational approximation.

Under the assumptions of the model introduced in Section 2, the mean field
approximation yields

$\displaystyle q\left(\boldsymbol{\theta}\right)=$
$\displaystyle\;q\left(\boldsymbol{\beta}\right)\cdot
q\left(\lambda_{0}^{2}\right)\cdot q\left(\lambda_{1}\right)\cdot\ldots\cdot
q\left(\lambda_{p}\right)\cdot q\left(\psi_{0}\right)\cdot
q\left(\varphi^{2}_{1}\right)\cdot\ldots\cdot
q\left(\varphi^{2}_{p}\right)\cdot$
$\displaystyle\;\;q\left(\boldsymbol{\gamma}\right)\cdot
q\left(\kappa^{2}_{1}\right)\cdot\ldots\cdot q\left(\kappa^{2}_{D}\right)\cdot
q\left(\tau^{2}\right)\cdot q\left(\zeta\right)\cdot
q\left(\sigma^{2}\right),$

where $q\left(\boldsymbol{\beta}\right)$ and
$q\left(\boldsymbol{\gamma}\right)$ denote the joint variational distribution
of $\beta_{0},\dots,\beta_{p}$ and $\gamma_{1},\dots,\gamma_{M}$,
respectively. At each iteration of the algorithm, the variational factors are
updated as

$\displaystyle q^{\star}\left(\boldsymbol{\beta}\right)=$
$\displaystyle\;\mathcal{N}_{p+1}\left(\boldsymbol{\mu}_{\beta}^{\star},\mathbb{E}_{\sigma^{2}}\left[\sigma^{2}\right]\boldsymbol{\Sigma}_{\beta}^{\star}\right),$
$\displaystyle\;\boldsymbol{\mu}_{\beta}^{\star}=\boldsymbol{\Sigma}_{\beta}^{\star}\mathbf{X}^{\intercal}\mathbf{y},\quad\boldsymbol{\Sigma}_{\beta}^{\star}=\left(\mathbf{X}^{\intercal}\mathbf{X}+\mathbb{E}_{\lambda_{0}^{2}\cdot\lambda\cdot\tau^{2}}\left[\tau^{-2}\boldsymbol{\Lambda}^{-2}\right]\right)^{-1},$
$\displaystyle q^{\star}\left(\lambda_{0}^{2}\right)=$
$\displaystyle\;\mathcal{IG}\left(1,a_{0}^{\star}\right),$
$\displaystyle\;a_{0}^{\star}=\mathbb{E}_{\zeta}\left[\psi_{0}^{-1}\right]+\frac{1}{2}\mathbb{E}_{\beta_{0}\cdot\sigma^{2}\cdot\tau^{2}}\left[\frac{\beta_{0}^{2}}{\sigma^{2}\tau^{2}}\right],$
$\displaystyle q^{\star}\left(\psi_{0}\right)=$
$\displaystyle\;\mathcal{IG}\left(1,k_{0}^{\star}\right),$
$\displaystyle\;k_{0}^{\star}=1+\mathbb{E}_{\lambda_{0}^{2}}\left[\lambda_{0}^{-2}\right],$
$\displaystyle q^{\star}\left(\lambda_{j}\right)\propto$
$\displaystyle\;\lambda_{j}^{-1}\exp\left\\{-\frac{a_{j}^{\star}}{\lambda_{j}^{2}}-b_{j}^{\star}\lambda_{j}^{2}+c_{j}^{\star}\lambda_{j}\right\\}\cdot\mathbb{I}_{\left(\lambda_{j}>0\right)},\quad
j=1,\dots,p,$
$\displaystyle\;a_{j}^{\star}=\frac{1}{2}\mathbb{E}_{\beta\cdot\sigma^{2}\cdot\tau^{2}}\left[\frac{\beta_{j}^{2}}{\sigma^{2}\tau^{2}}\right],\quad
b_{j}^{\star}=\frac{1}{2s_{0}^{2}}\mathbb{E}_{\varphi^{2}}\left[\varphi_{j}^{-2}\right],\quad
c_{j}^{\star}=\frac{1}{s_{0}^{2}}\mathbf{z}_{j}^{\intercal}\mathbb{E}_{\gamma\cdot\varphi^{2}}\left[\frac{\boldsymbol{\gamma}}{\varphi_{j}^{2}}\right],$
$\displaystyle q^{\star}\left(\varphi_{j}^{2}\right)=$
$\displaystyle\;\mathcal{IG}\left(1,d_{j}^{\star}\right),$
$\displaystyle\;d_{j}^{\star}=\frac{1}{2}+\frac{1}{2s_{0}^{2}}\mathbb{E}_{\lambda\cdot\gamma}\left[\left(\lambda_{j}-\mathbf{z}_{j}^{\intercal}\boldsymbol{\gamma}\right)^{2}\right],$
$\displaystyle q^{\star}\left(\boldsymbol{\gamma}\right)=$
$\displaystyle\;\mathcal{N}_{M}\left(\boldsymbol{\mu}_{\gamma}^{\star},s_{0}^{2}\boldsymbol{\Sigma}_{\gamma}^{\star}\right),$
$\displaystyle\;\boldsymbol{\mu}_{\gamma}^{\star}=\boldsymbol{\Sigma}_{\gamma}^{\star}\mathbf{Z}^{\intercal}\mathbb{E}_{\varphi^{2}\cdot\gamma}\left[\boldsymbol{\Phi}^{-2}\boldsymbol{\lambda}\right],\quad\boldsymbol{\Sigma}_{\gamma}^{\star}=\left(\mathbf{Z}^{\intercal}\mathbb{E}_{\varphi^{2}}\left[\boldsymbol{\Phi}^{-2}\right]\mathbf{Z}+s_{0}^{2}\mathbb{E}_{\kappa^{2}}\left[\mathbf{D}_{\kappa}^{-1}\right]\right)^{-1},$
$\displaystyle q^{\star}\left(\kappa_{d}^{2}\right)=$
$\displaystyle\;\mathcal{IG}\left(e_{d}^{\star},f_{d}^{\star}\right),\quad
d=1,\dots,D,$ $\displaystyle\;e_{d}^{\star}=a_{d}+\frac{m_{d}}{2},\quad
f_{d}^{\star}=b_{d}+\frac{1}{2}\mathbb{E}_{\gamma}\left[\boldsymbol{\gamma}_{d}^{\intercal}\boldsymbol{\gamma}_{d}\right],$
$\displaystyle q^{\star}\left(\tau^{2}\right)=$
$\displaystyle\;\mathcal{IG}\left(\frac{p}{2}+1,g^{\star}\right),$
$\displaystyle\;g^{\star}=\mathbb{E}_{\zeta}\left[\zeta^{-1}\right]+\frac{1}{2}\mathbb{E}_{\beta_{0}\cdot\lambda_{0}^{2}\cdot\tau^{2}}\left[\frac{\beta_{0}^{2}}{\tau^{2}\lambda_{0}^{2}}\right]+\frac{1}{2}\sum_{j=1}^{p}\mathbb{E}_{\beta\cdot\sigma^{2}\cdot\lambda}\left[\frac{\beta_{j}^{2}}{\sigma^{2}\lambda_{j}^{2}}\right],$

$\displaystyle q^{\star}\left(\zeta\right)=$
$\displaystyle\;\mathcal{IG}\left(1,h^{\star}\right),$
$\displaystyle\;h^{\star}=1+\mathbb{E}_{\tau^{2}}\left[\tau^{-2}\right],$
$\displaystyle q^{\star}\left(\sigma^{2}\right)=$
$\displaystyle\;\mathcal{IG}\left(v+\frac{n+p+1}{2},l^{\star}\right),$
$\displaystyle\;l^{\star}=q+\frac{1}{2}\left(\mathbb{E}_{\beta}\left[\|\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\|_{2}^{2}\right]+\mathbb{E}_{\beta_{0}\cdot\lambda_{0}^{2}\cdot\tau^{2}}\left[\frac{\beta_{0}^{2}}{\tau^{2}\lambda_{0}^{2}}\right]+\sum_{j=1}^{p}\mathbb{E}_{\beta\cdot\lambda\cdot\tau^{2}}\left[\frac{\beta_{j}^{2}}{\tau^{2}\lambda_{j}^{2}}\right]\right),$

where the expectations are taken with respect to the updated variational
family $q^{\star}\left(\boldsymbol{\theta}\right)$. The variational lower
bound $\mathcal{L}$ of the marginal distribution $p\left(\mathbf{y}\right)$ is
computed as

$\displaystyle\mathcal{L}\propto$
$\displaystyle\;\frac{1}{2}\left(\log|\boldsymbol{\Sigma}_{\beta}^{\star}|+\log|\boldsymbol{\Sigma}_{\gamma}^{\star}|\right)+\frac{p+1}{2}\mathbb{E}_{\sigma^{2}}\left[\log\sigma^{2}\right]-\left(v+\frac{n+p+1}{2}\right)\log
l^{\star}+$ $\displaystyle\;\sum_{j=1}^{p}\left(\log s_{j}-\log
k_{j}+a_{j}^{\star}\mathbb{E}_{\lambda}\left[\lambda_{j}^{-2}\right]+b_{j}^{\star}\mathbb{E}_{\lambda}\left[\lambda_{j}^{2}\right]-c_{j}^{\star}\mathbb{E}_{\lambda}\left[\lambda_{j}\right]-\log
d_{j}^{\star}\right)-$ $\displaystyle\;\sum_{d=1}^{D}e_{d}^{\star}\log
f_{d}^{\star}-\left(\frac{p}{2}+1\right)\log g^{\star}-\log h^{\star}-\log
a_{0}^{\star}-\log k_{0}^{\star},$ (5)

where
$k_{j}=1-\mathbb{P}\left(\mathcal{N}\left(\mathbf{z}_{j}^{\intercal}\boldsymbol{\gamma},s_{0}^{2}\varphi_{j}^{2}\right)<0\right)$
and $s_{j}$ is the unknown normalizing constant of parameters $\lambda_{j}$ in
(4). Each component is derived in Section D.1 of Supporting Information.

As opposed to the Gibbs sampler, the Variational algorithm does not require
the explicit inversion of the $p\times p$ matrix
$\boldsymbol{\Sigma}_{\beta}^{\star}$, as only quantities
$\boldsymbol{\Sigma}_{\beta}^{\star}\mathbf{X}^{\intercal}\mathbf{y}$,
$\text{diag}\left(\boldsymbol{\Sigma}_{\beta}^{\star}\right)$ and
$\mathbf{X}\boldsymbol{\Sigma}_{\beta}^{\star}\mathbf{X}^{\intercal}$ are
needed. These can be efficiently evaluated with complexity
$\mathcal{O}\left(n^{2}p\right)$ following the work of Münch et al., (2019).
The algebraic details are shown in Section E.1 of Supporting Information. The
bottleneck of the algorithm is the computation of the terms $\log s_{j}$,
$\mathbb{E}_{\lambda}\left[\lambda_{j}\right]$,
$\mathbb{E}_{\lambda}\left[\lambda_{j}^{2}\right]$ and
$\mathbb{E}_{\lambda}\left[\lambda_{j}^{-2}\right]$, since no closed-form is
available. These can be evaluated in parallel following the procedure
described in Section C.1 of Supporting Informatiion.

##  5 Simulation study

In this section we empirically show the quality of the Variational
approximation to the joint posterior distribution and the benefit of using co-
data for variable selection. We use a simulation scheme that resembles the one
by Johnson, (2013). In particular, we assess the effectiveness of the
approximation on low to moderate $p$ problems, whereas variable selection is
evaluated on higher dimensional settings. We refer to Section F.1 of
Supporting Information for further details on the data generation scheme.

Variable selection is evaluated with the average area under the ROC curve
($AUC$) of the selected variable, where the posterior inclusion probability of
the generic covariate $j$ is evaluated with the quantity
$\mathbb{E}_{\lambda_{j}}\left[\lambda_{j}^{2}/\left(1+\lambda_{j}^{2}\right)\right]$
(Carvalho et al.,, 2010).

The hyperparameters for $\sigma^{2}$ and $\kappa_{1}^{2}$ are set to
$(v,q)=(1,10)$ and $(a_{1},b_{1})=(1,10)$, respectively. The Gibbs sampler is
run for $B=5000$ iterations with a burn-in period of $bn=2500$ iterations,
while the Variational algorithm is stopped either if the increase of the lower
bound is less than $\epsilon=0.001$ or if the maximum number of iteration
$B=1000$ is reached.

1\. Performance of variable selection.

Here we analyse how the model behave as the degree of prior information
varies. We study the cases $n=(50,100,150)$ and $p=(500,1000,1500)$ and we set
the number of true non-zero coefficients to $p_{0}=30$. To study how the model
borrows information from the co-data, we simulate co-data sources representing
five degrees of prior information:

  * $G0\big{)}$ no co-data set-up: we include in the co-data matrix only the intercept, therefore $\mathbf{Z}=\mathbf{1}_{p}$;

  * $G1\big{)}$ non-informative set-up: a binary co-data source is included in the model by randomly selecting $100$ regressors, therefore the co-data matrix $\mathbf{Z}$ is created from the binary vector $\mathbf{z}\in\left\\{0,1\right\\}^{p}$, where $z_{j}=1$ if the $j$-th variable is selected, $z_{j}=0$ otherwise;

  * $G2\big{)}$ weakly informative set-up: a binary co-data source is included in the model by randomly selecting $20$ of the true non-zero regressors and $80$ of the true zero regressors, therefore the co-data matrix $\mathbf{Z}$ is created from the binary vector $\mathbf{z}\in\left\\{0,1\right\\}^{p}$, where $z_{j}=1$ if the $j$-th variable is selected, $z_{j}=0$ otherwise;

  * $G3\big{)}$ informative set-up: a binary co-data source is included in the model by randomly selecting $20$ of the true non-zero regressors and $10$ of the true zero regressors, therefore the co-data matrix $\mathbf{Z}$ is created from the binary vector $\mathbf{z}\in\left\\{0,1\right\\}^{p}$, where $z_{j}=1$ if the $j$-th variable is selected, $z_{j}=0$ otherwise;

  * $G4\big{)}$ perfect co-data information set-up: the co-data matrix $\mathbf{Z}$ is created from the binary vector $\mathbf{z}\in\left\\{0,1\right\\}^{p}$, where $z_{j}=1$ if $\beta_{j}^{0}\neq 0$, $z_{j}=0$ otherwise.

Figure 1: Variable selection with Variational algorithm for co-data scenarios
$G0$ to $G4$; $AUC$ is averaged over $25$ replicates of the experiments.

The results are shown in Figure 1. The model is able to learn from the
auxiliary information and the variable selection performance improves when the
co-data is actually informative. The $AUC$, indeed, increases alongside the
magnitude of prior information and the variables are perfectly selected when
the co-data provides perfect information (case $G4$). On the other hand,
random co-data is associated with a small loss of performance and leads to
lower scores of $AUC$. In particular, if we compare the case with random co-
data ($G1$) and no co-data ($G0$), the latter performs slightly better. This
difference, however, vanishes as the number of covariates increases. Finally,
the smallest increase in $AUC$, as expected, is between $G2$ and $G3$, since
these cases represent the most similar degree of co-data information.

2\. Gibbs sampler vs Variational inference.

The accuracy of the Variational approximation is evaluated by comparing it to
the Gibbs sampler in terms of variable selection and mean squared error
($MSE$) between $\boldsymbol{\beta}$ estimates. We analysed the cases
$n=(50,100)$ and $p=(75,125,200)$ and we set $p_{0}=30$. We consider four
degrees of prior information and we refer to Section F.2 of Supporting
Information for further details on co-data simulation.

The results of variable selection are shown in Figure 5 of Supporting
Information. The two methods behave similarly in all the considered cases. In
particular, the results are equal for the lowest dimensional case, whereas the
Variational approximation approaches the performance of the Gibbs sampler when
$p$ increases. Table 1 reports the details of the similarity between the two
methods in terms of $\boldsymbol{\beta}$ estimate. Gibbs sampling and
Variational algorithm provide similar estimates for the regression parameters,
as the mean of the $MSE$ is close to zero for all the considered cases. In
particular, when $n$ is large and the co-data are strongly informative the two
methods provide almost the same estimate.

In Figure 6 of Supporting Information we show the average estimated standard
deviations of regression parameters $\boldsymbol{\beta}$ across the different
cases of co-data information. The results follows the considerations made in
MacKay et al., (2003); Wang and Titterington, (2005); Giordano et al., (2017),
as the Variational approximation underestimates the posterior variances when
compared to Gibbs sampling. However, as highlighted in Blei and Jordan, (2006)
and from the simulation results, this lack of accuracy does not necessarily
affect the global performance of the model.

The Gibbs sampler presents one main drawback: for some settings of $n$ and
$p$, the sampling method presented in Section B.1 of Supporting Information
suffers from low acceptance probability. Therefore, the efficiency of the
Gibbs sampler is heavily affected in a negative way. Given the good results of
the Variational approximation and the poor performance in terms of
computational efficiency of the Gibbs sampler, we rely on the former for the
following large $p$ applications.

Table 1: Mean of $MSE_{\boldsymbol{\beta}}$ between Gibbs sampling and the
Variational algorithm; the mean is evaluated over $10$ replicates of each
experiment.
$MSE_{\boldsymbol{\beta}}=\|\boldsymbol{\beta}_{\text{FB}}-\boldsymbol{\beta}_{\text{GS}}\|_{2}^{2}$
---
| $p=75$ | $p=125$ | $p=200$ |
$G0$ | 0.0079 | 0.2592 | 0.0365 | $n=50$
$G1$ | 0.0074 | 0.1034 | 0.0070
$G2$ | 0.0091 | 0.0647 | 0.0621
$G3$ | 0.0032 | 0.1153 | 0.0787
$G0$ | 0.0006 | 0.0268 | 0.0217 | $n=100$
$G1$ | 0.0007 | 0.0023 | 0.0212
$G2$ | 0.0005 | 0.0046 | 0.0023
$G3$ | $<0.0001$ | 0.0037 | 0.0001

##  6 Application to real data

In this section we present two genomics applications. The first allows to
evaluate our method for binary co-data in a multiple linear regression
context. The second involves multiple co-data sources in a classification
setting with a large number of variables.

In Section 5 we show that the thresholding variable selection works well when
the covariates are sampled independently. However, the posterior probabilities
are treated separately and the optimal threshold is (almost) never $0.5$. For
this reason in this section we also consider the more recent selection
procedure called decoupling shrinkage and selection (DSS) proposed by Hahn and
Carvalho, (2015). This method was introduced to deal with potentially very
strong correlations, as they are present in many genomics datasets. The
authors propose a posterior variable selection summary based on the posterior
mean of the predictors which results in a sequence of sparse models. Shrinkage
can be achieved with any prior distribution and is ‘decoupled’ from the
selection approach based on the posterior distribution. DSS method relies on
the optimization of a loss function which balances the prediction error and
the sparseness of the solution. Given an estimate $\hat{\boldsymbol{\beta}}$,
the solution is obtained via adaptive LASSO by solving the following
optimization problem,

$\boldsymbol{\theta}^{DSS}=\text{argmin}_{\boldsymbol{\theta}}\;\frac{1}{n}\|\mathbf{X}\hat{\boldsymbol{\beta}}-\mathbf{X}\boldsymbol{\theta}\|_{2}^{2}+\lambda\sum_{j=0}^{p}\frac{|\theta_{j}|}{|\hat{\beta}_{j}|},$

where the smoothing parameter $\lambda$ operates as a thresholding parameter
and can be estimated with cross-validation over a set of values (grid search).
The authors advocate this method over thresholding mainly because it naturally
handles multi-collinearity.

###  6.1 Case study 1: p38MAPK pathway

Following Kpogbezan et al., (2019), the model is tested on the p38MAPK pathway
dataset. We investigate the effect of single nucleotide polymorphisms (SNPs)
on the genes in the pathway. A subset of the data collected in the GEUVADIS
project (Lappalainen et al.,, 2013) is considered. For each of the $99$ genes
we collect a different number $p_{t}$ of SNPs, $t=1,\dots,99$, with minimum
$p_{t}=56$ and maximum $p_{t}=1169$. After a first sub-selection, $n=373$ RNA-
Seq samples are obtained from the 1000 Genomes Project. Since it is believed
that SNPs within the gene’s range have stronger influence on the gene’s
expression, a binary co-data source
$\mathbf{z}_{t}\in\left\\{0,1\right\\}^{p_{t}}$ is included in the analysis,
where $z_{t,j}=1$ if the $j$-th SNP related to gene $t$ is located inside the
gene’s range and $z_{t,j}=0$ otherwise.

The hyperparameters for $\sigma^{2}$ and $\kappa_{1}^{2}$ are set to
$(v,q)=(1,10)$ and $(a_{1},b_{1})=(1,10)$. The algorithm is stopped either if
the increase of the variational lower bound is less than $\epsilon=0.001$ or
if the maximum number of iteration $B=1000$ is reached.

Figure 2: Relative reduction of MSE (rrMSE) of LASSO (darkgrey dots), infHS
regression with DSS selection procedure (black dots) and infHS regression with
thresholding selection procedure (lightgrey dots) for all the $99$ genes. For
each panel, the maximum number of selected SNPs is fixed to $1$, $3$, $5$, and
$10$.

To evaluate the prediction performance of the model, the data are divided in
training-set ($n_{1}=249$) and test-set ($n_{2}=124$). We rely on the relative
reduction of the mean squared error (rrMSE) for all the $99$ genes. The rrMSE
for the $t$-th regression model is defined as

$\text{rrMSE}_{t}=1-\frac{\text{MSE}_{t}}{\text{MSE}_{0}},$

where $\text{MSE}_{0}$ and $\text{MSE}_{t}$ are the mean squared error of the
null model and the mean squared error of the $t$-th linear model,
respectively. The mean of $\mathbf{y}$ is centered around zero, therefore we
do not include the intercept in the null model. Note that larger values of
$\text{rrMSE}_{t}$ are associated with better predictive performance. On the
contrary, negative values denote worse predictive performance than the null
model.

We test our method against LASSO regression (Tibshirani,, 1996). The results
are shown is Figure 2 for different degrees of sparseness. For the infHS
regression, the SNPs are selected with both the thresholding procedure used in
Section 5 (lightgrey dots) and the DSS method (black dots). The latter
provides a better prediction for all genes when compared with the thresholding
variable selection. This is justified by the strong correlations among the
SNPs. When comparing infHS and LASSO regressions, instead, the considerations
made in Kpogbezan et al., (2019) hold. The results of the LASSO are more
noisy, likely due to less shrinkage of the near-zero coefficients, whereas
most of the $\text{rrMSE}_{t}$ are concentrated around zero for the infHS
model. The effect of the SNPs are relevant only for a small number of genes:
the largest values of $\text{rrMSE}_{t}$ are associated to genes $50$, $61$
and $93$. Both infHS and LASSO are able to capture these effects, however the
latter is more prone to negative values of $\text{rrMSE}_{t}$. The LASSO gives
better results for gene $98$. Note that the method by Kpogbezan et al., (2019)
gives fairly similar results to ours, as it is also based on the Horseshoe.
That method, however, is much more limited in use, as it only handles one
discrete co-data source and only continuous outcomes.

Our method take $68$ minutes to estimate all the $99$ regressions and to
evaluate the different variable selection approaches. The algorithm is run on
a x64 Windows 11 operating system.

###  6.2 Case study 2: methylation data

The model is tested on the methylation dataset of Verlaat et al., (2018),
which contains methylation profiles of self-collected cervicovaginal lavages
corresponding to $28$ women with normal cervix and $36$ women with high-grade
precursor lesions (CIN3), for a total of $n=64$ samples. In order to improve
the diagnostic classification we include in the analysis $D=5$ co-data
sources: the standard deviations of the probes, p-values from the previous
study, a binary variable differentiating the hypo-methylated and hyper-
methylated probes, a categorical variable with $6$ categories denoting the
genomic region of each probe and the probes’ means in cancer cells. See
Section C.2 of Supporting Information for theorical details of the Variational
approximation to the joint posterior distribution in a probit regression
model. Based on the previous results, we consider the probes with an external
p-value $p_{\alpha}<0.005$, resulting in a total of $p=11251$ probes, where
each probe measures the methylation of a unique location on the genome.

(a)

(b)

Figure 3: Results of the LOOCV. Left panel: ROC curves for Ridge regression
(lightgrey), LASSO (darkgrey), the informative Horseshoe regression (black)
and the informative Horseshoe with DSS variable selection procedure (dotted-
black); right panel: a and b. predicted probabilities for cases ($y_{i}=1$)
and controls ($y_{i}=0$) with infHS-DSS (black) and LASSO (darkgrey) in
decreasing order of LASSO prediction; c and d. predicted probabilities for
cases and controls with infHS-DSS (black) and ordinary infHS (darkgrey) in
decreasing order of infHS prediction.

We test our model against the ordinary Ridge regression and the LASSO. The
hyperparameters are set to $\left(a_{d},b_{d}\right)=(1,10)$ for $d=1,\dots,5$
and the predictive performances are evaluated by leave-one-out cross-
validation (LOOCV). As for case study 1, the algorithm is stopped either if
the increase of the variational lower bound is less than $\epsilon=0.001$ or
if the maximum number of iteration $B=1000$ is reached. We include in the
analysis also the thresholding and the DSS versions of infHS model. The former
does not exclude any of the covariates from the model, since all the posterior
inclusion probabilities are greater than 0.5. Therefore the results of this
approach are not shown. We also fit our model without co-data; due to many
cases ($y_{i}=1$) corresponding to very low predicted probabilities, the
results are not competitive to those of infHS, and hence not shown here.

The left panel of Figure 3 shows the ROC curves for the considered
competitors, where the predicted probabilities for infHS are computed as
$p_{i}=\Phi\left(\mathbf{x}_{i}^{\intercal}\boldsymbol{\beta}\right)$. InfHS-
DSS gives the best forecasting results with $\text{AUC}=0.903$, greatly
improving compared to its ordinary version with $\text{AUC}=0.752$; the LASSO
performs better than infHS with $\text{AUC}=0.773$ and provides the sparsest
model, with an average of $8$ variables included as opposed to the $49$
estimated by infHS-DSS. Among these $49$ probes, $37$ of them are selected in
more than $70\%$ of the folds, whereas only $10$ appear in all. The ordinary
Ridge does not compete in terms of predictive performance
($\text{AUC}=0.671$). Note that GRridge co-data method (Van de Wiel et al.,,
2016), which is one of the few methods that also tackles multiple co-data
sources, was previously reported to achieve an $\text{AUC}=0.77$ (Verlaat et
al.,, 2018).

The right panel of Figure 3 shows the predicted probabilities of infHS-DSS
model against the LASSO and the ordinary infHS. When compared to the LASSO
(figures a and b, values in decreasing order of LASSO prediction), infHS-DSS
gives similar predictions for the cases ($y_{i}=1$), with some of the
observations achieving higher predictive scores and other assuming lower ones;
on the other hand, it provides significantly lower scores for the controls
($y_{i}=0$), reducing the false positive rate and improving the overall
performance. The impact of the DSS procedure on the prediction can be assessed
by analyzing figures c and d, which compares infHS-DSS to infHS model (values
in decreasing order of infHS prediction). The former completely overwhelms the
latter, as it gives higher probabilities for almost all cases and reduces the
predicted probabilities for the highest scores of the controls. To sum up
these results we compute the mean absolute difference between the true labels
and the predicted probabilities. The results are shown in Table 2, where
infHS-DSS provides the best performance.

Table 2: Mean absolute difference between true labels ($y_{i}=0$ or $y_{i}=1$) and the predicted probabilities ($p_{i}=\Phi\left(\mathbf{x}_{i}^{\intercal}\boldsymbol{\beta}\right)$). | Ridge | LASSO | infHS | infHS-DSS
---|---|---|---|---
$\sum_{i=1}^{64}\frac{|y_{i}-p_{i}|}{64}$ | 0.466 | 0.410 | 0.444 | 0.299

We summarize the co-data information in Table 3, which shows the distribution
of the overall co-data sources and the distribution of the $p_{sel}=37$ most
selected probes. The co-data related to the the probes’ means in cancer cells
show the strongest difference in distribution between the original and the
selected features, suggesting that this source is the most informative in
terms of prior information. No evident differences are found in the other co-
data sources.

The computational time is around 5 minutes for the estimation of infHS. The
method is tested on a x64 Windows 11 operating system and the local variances
are evaluated in parallel with 4 cores.

Table 3: Summary of the co-data estimates and the co-data distribution in the $p_{sel}=37$ most selected probes and in the original population. Co-data | Selected probes distribution | Original distribution
---|---|---
($p_{sel}=37$) | ($p=11251$)
external $p$-value | mean: $0.0016$ | mean: $0.0019$
sd: $0.0014$ | sd: $0.0015$
range: ($2.4\cdot 10^{-6},0.005$) | range: ($9.9\cdot 10^{-11},0.005$)
probes sd | mean: $0.035$ | mean: $0.038$
sd: $0.022$ | sd: $0.019$
range: ($0.014,0.129$) | range: ($0.009,0.320$)
genomic region | Distant: $35.1\%$ | Distant: $31.7\%$
Island: $43.2\%$ | Island: $37.3\%$
N shelf: $5.4\%$ | N shelf: $4.6\%$
N shore $10.8\%$ | N shore: $12.4\%$
S shelf: $0\%$ | S shelf: $4\%$
S shore: $5.4\%$ | S shore: $10\%$
degree of methylation | Hypo-methylated: $57\%$ | Hypo-methylated: $56\%$
Hyper-methylated: $43\%$ | Hyper-methylated: $44\%$
cancer cells mean | mean: $0.150$ | mean: $-0.099$
sd: $1.116$ | sd: $1.137$
range: ($-1.28,3.34$) | range: ($-4.89,8.10$)

##  7 Discussion

We introduced a novel Bayesian regression approach for high-dimensional data
able to learn from auxiliary prior information, i.e. co-data. We showed that
both prediction and variable selection benefit from the inclusion of co-data,
when these are actually informative. The model allows for both continuous and
binary outcome, as well as both continuous and discrete co-data sources. In
particular, we provided a flexible method that estimates multiple co-data
coefficients jointly, contrary to the previous method of Van Nee et al.,
(2021) that models each source separately.

We discussed a full Bayesian approach, for which we developed a Gibbs sampler
to update each parameter iteratively by sampling from their full-conditional
distributions. To do so, we introduced a novel rejection sampling method to
sample the local variances, which showed a non-closed form target density.
Eventually, we proposed a Variational approximation to the joint posterior
distribution and applied the CAVI algorithm to optimize the target density.
This latter method is suited for high-dimensional problems, as it does not
require the explicit inversion of the $p\times p$ posterior covariance matrix
$\boldsymbol{\Sigma}_{\beta}^{\star}$. In particular, only its diagonal is
required and we implemented the methods proposed in Münch et al., (2019) to
efficiently achieve this. Beside this, another computational advantage is the
parallel evaluation of the local variances. The evaluation of these parameters
represents the computational bottleneck of the method: the acceptance
probability of the rejection sampling is small for some settings of the
parameters, whereas the numerical integration required in the Variational
approximation is the most expensive step of the algorithm.

The Variational inference is less useful than the Gibbs sampler for posterior
inference, as it provides posterior means point estimates and underestimates
the posterior variances (MacKay et al.,, 2003; Wang and Titterington,, 2005;
Turner and Sahani,, 2011; Giordano et al.,, 2017). This lack of accuracy, on
the other hand, does not necessarily affect the predictive performance of the
model (Blei and Jordan,, 2006). Therefore, the Variational algorithm should be
used for large $p$ datasets, whereas the Gibbs sampling could be useful when
the interest is in the posterior inference (credible intervals) and the number
of covariates is moderate.

A possible limitation of infHS for some applications occurs when a categorical
co-data source contains one very strong, relatively small co-data group. Under
this circumstance, the sparsity assumption may not be realistic for this
particular group of variables. For such applications, an interesting extension
of infHS would be to allow a dense prior for a small group of variables, for
which one expects a particularly relevant prior evidence.

Another extension to allow more flexibility would be the specification of the
prior $\lambda_{j}\sim\text{Half-}t(v)$, with $v>1$. Here $v=1$ since the
purpose is to apply the Horseshoe prior. In Biswas et al., (2021) the authors
develop coupling techniques for high-dimensional regression with this
particular prior and argue that larger values of $v$ can affect the
statistical and computational performance of the model.

To conclude, we provided one of the fastest implementation of Horseshoe
regression able to learn from auxiliary information. The R code for infHS and
the datasets used are available at https://github.com/cbusatto/infHS.

## Acknowledgements

The first author was supported also by the project THE - Tuscany Health
Ecosystem. The authors thank professor Francesco C. Stingo for his help with
the start of this project.

## References

  * Albert and Chib, (1993) Albert, J. H. and Chib, S. (1993).  Bayesian analysis of binary and polychotomous response data.  Journal of the American Statistical Association, 88(422):669–679.
  * Beal, (2003) Beal, M. J. (2003).  Variational algorithms for approximate Bayesian inference.  University of London, University College London (United Kingdom).
  * Bhattacharya et al., (2016) Bhattacharya, A., Chakraborty, A., and Mallick, B. K. (2016).  Fast sampling with Gaussian scale-mixture priors in high-dimensional regression.  Biometrika, 103 4:985–991.
  * Bishop and Nasrabadi, (2006) Bishop, C. M. and Nasrabadi, N. M. (2006).  Pattern recognition and Machine Learning, volume 4.  Springer.
  * Biswas et al., (2021) Biswas, N., Bhattacharya, A., Jacob, P. E., and Johndrow, J. E. (2021).  Coupling-based convergence assessment of some Gibbs samplers for high-dimensional Bayesian regression with shrinkage priors.  Journal of the Royal Statistical Society. Series B (Statistical Methodology).
  * Blei and Jordan, (2006) Blei, D. M. and Jordan, M. I. (2006).  Variational inference for Dirichlet process mixtures.  Bayesian Analysis, 1(1):121–143.
  * Blei et al., (2017) Blei, D. M., Kucukelbir, A., and McAuliffe, J. D. (2017).  Variational inference: A review for statisticians.  Journal of the American Statistical Association, 112(518):859–877.
  * Boonstra et al., (2013) Boonstra, P. S., Taylor, J. M. G., and Mukherjee, B. (2013).  Incorporating auxiliary information for improved prediction in high-dimensional datasets: an ensemble of shrinkage approaches.  Biostatistics (Oxford, England), 14(2):259–272.
  * Busatto and Stingo, (2023) Busatto, C. and Stingo, F. C. (2023).  Inference of multiple high-dimensional networks with the graphical horseshoe prior.
  * Carvalho et al., (2010) Carvalho, C. M., Polson, N. G., and Scott, J. G. (2010).  The horseshoe estimator for sparse signals.  Biometrika, 97:465–480.
  * Geweke, (1993) Geweke, J. (1993).  Bayesian treatment of the independent Student-t linear model.  Journal of Applied Econometrics, 8:S19–S40.
  * Giordano et al., (2017) Giordano, R., Broderick, T., and Jordan, M. (2017).  Covariances, robustness, and variational bayes.  Journal of Machine Learning Research, 19.
  * Hahn and Carvalho, (2015) Hahn, P. R. and Carvalho, C. M. (2015).  Decoupling shrinkage and selection in Bayesian linear models: A posterior summary perspective.  Journal of the American Statistical Association, 110(509):435–448.
  * Johndrow et al., (2020) Johndrow, J. E., Orenstein, P., and Bhattacharya, A. (2020).  Scalable approximate MCMC algorithms for the horseshoe prior.  Journal of Machine Learning Research, 21:73:1–73:61.
  * Johnson, (2013) Johnson, V. E. (2013).  On numerical aspects of Bayesian model selection in high and ultrahigh-dimensional settings.  Bayesian analysis (Online), 8(4):741.
  * Jordan et al., (1999) Jordan, M., Ghahramani, Z., Jaakkola, T., and Saul, L. (1999).  An introduction to variational methods for graphical models.  Machine Learning, 37:183–233.
  * Kpogbezan et al., (2019) Kpogbezan, G., van de Wiel, M., van Wieringen, W., and van der Vaart, A. (2019).  Incorporating prior information and borrowing information in high-dimensional sparse regression using the horseshoe and Variational Bayes.
  * Lappalainen et al., (2013) Lappalainen, T., Sammeth, M., Friedländer, M., Hoen, P., Monlong, J., Rivas, M., Gonzàlez Porta, M., Kurbatova, N., Griebel, T., Ferreira, P., Barann, M., Wieland, T., Greger, L., Iterson, M., Almlöf, J., Ribeca, P., Pulyakhina, I., Esser, D., Giger, T., and Dermitzakis, E. (2013).  Transcriptome and genome sequencing uncovers functional variation in humans.  Nature, 501.
  * Lee, (2022) Lee, S. Y. (2022).  Gibbs sampler and coordinate ascent variational inference: A set-theoretical review.  Communications in Statistics - Theory and Methods, 51(6):1549–1568.
  * MacKay et al., (2003) MacKay, D. J., Mac Kay, D. J., et al. (2003).  Information theory, inference and learning algorithms.  Cambridge university press.
  * Makalic and Schmidt, (2016) Makalic, E. and Schmidt, D. F. (2016).  A simple sampler for the horseshoe estimator.  IEEE Signal Processing Letters, 23(1):179–182.
  * Münch et al., (2019) Münch, M., Peeters, C., van der Vaart, A., and van de Wiel, M. (2019).  Adaptive group-regularized logistic elastic net regression.  Biostatistics (Oxford, England), 22.
  * Neuenschwander et al., (2016) Neuenschwander, B., Roychoudhury, S., and Schmidli, H. (2016).  On the use of co-data in clinical trials.  Statist. Biopharm. Res., 8(3):345–354.
  * Polson et al., (2013) Polson, N. G., Scott, J. G., and Windle, J. (2013).  Bayesian inference for logistic models using Pólya–Gamma latent variables.  Journal of the American Statistical Association, 108(504):1339–1349.
  * Salimans et al., (2015) Salimans, T., Kingma, D. P., and Welling, M. (2015).  Markov Chain Monte Carlo and variational inference: bridging the gap.  In Proceedings of the 32nd International Conference on International Conference on Machine Learning - Volume 37, ICML’15, page 1218–1226.
  * Simon et al., (2013) Simon, N., Friedman, J., Hastie, T., and Tibshirani, R. (2013).  A sparse-group Lasso.  J Comput Graph Stat, 22(2):231–245.
  * Tai and Pan, (2007) Tai, F. and Pan, W. (2007).  Incorporating prior knowledge of predictors into penalized classifiers with multiple penalty terms.  Bioinformatics, 23(14):1775–1782.
  * Tibshirani, (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the Lasso.  Journal of the Royal Statistical Society. Series B (Methodological), 58(1):267–288.
  * Turner and Sahani, (2011) Turner, R. E. and Sahani, M. (2011).  Two problems with variational expectation maximisation for time-series models.  In Bayesian Time series models, pages 109–130. Cambridge University Press.
  * Van de Wiel et al., (2016) Van de Wiel, M. A., Lien, T. G., Verlaat, W., van Wieringen, W. N., and Wilting, S. M. (2016).  Better prediction by use of co-data: adaptive group-regularized ridge regression.  Statistics in medicine, 35(3):368–381.
  * Van de Wiel et al., (2019) Van de Wiel, M. A., Te Beest, D. E., and Münch, M. M. (2019).  Learning from a lot: Empirical Bayes for high-dimensional model-based prediction.  Scandinavian Journal of Statistics, 46(1):2–25.
  * Van Nee et al., (2021) Van Nee, M. M., Wessels, L. F., and van de Wiel, M. A. (2021).  Flexible co-data learning for high-dimensional prediction.  Statistics in Medicine, 40(26):5910–5925.
  * Velten and Huber, (2019) Velten, B. and Huber, W. (2019).  Adaptive penalization in high-dimensional regression and classification with external covariates using Variational Bayes.  Biostatistics, 22(2):348–364.
  * Verlaat et al., (2018) Verlaat, W., Snoek, B. C., Heideman, D. A., Wilting, S. M., Snijders, P. J., Novianti, P. W., van Splunter, A. P., Peeters, C. F., van Trommel, N. E., Massuger, L. F., et al. (2018).  Identification and validation of a 3-gene methylation classifier for hpv-based cervical screening on self-samples.  Clinical Cancer Research, 24(14):3456–3464.
  * Wang and Titterington, (2005) Wang, B. and Titterington, D. M. (2005).  Inadequacy of interval estimates corresponding to variational Bayesian approximations.  In Proceedings of the Tenth International Workshop on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research, pages 373–380. PMLR.
  * Yuan and Lin, (2006) Yuan, M. and Lin, Y. (2006).  Model selection and estimation in regression with grouped variables.  Journal of the Royal Statistical Society Series B, 68:49–67.

##  Appendix A Web Appendix A

###  A.1 Gibbs step for local variance $\lambda_{0}^{2}$

The half-Cauchy prior for the local shrinkage parameter $\lambda_{0}$ of the
intercept $\beta_{0}$ is not conjugated to the variance in a linear regression
model with normal errors. We rely on the data-augmentation step proposed in
Makalic and Schmidt, (2016) in order to easily and efficiently update
parameter $\lambda_{0}$.

######  Proposition 4.

If
$\lambda_{0}^{2}\mid\psi_{0}\sim\mathcal{IG}\left(\frac{1}{2},\frac{1}{\psi_{0}}\right)$
and $\psi_{0}\sim\mathcal{IG}\left(\frac{1}{2},1\right)$, then
$\lambda_{0}\sim\mathcal{C}^{+}(0,1)$ (Makalic and Schmidt,, 2016).

The inverse-Gamma distribution is conjugated to itself and to the local scale
parameter, therefore a closed-form full-conditional is available and a Gibbs
step can be implemented. The full-conditional distributions of
$\lambda_{0}^{2}$ and $\psi_{0}$ are

$\displaystyle\lambda_{0}^{2}\mid\beta_{0},\tau^{2},\sigma^{2}\sim$
$\displaystyle\;\mathcal{IG}\left(1,\frac{1}{\psi_{0}}+\frac{\beta_{0}^{2}}{2\sigma^{2}\tau^{2}}\right),$
$\displaystyle\psi_{0}\mid\lambda_{0}^{2}\sim$
$\displaystyle\;\mathcal{IG}\left(1,1+\frac{1}{\lambda_{0}^{2}}\right).$

##  Appendix B Web Appendix B

###  B.1 Rejection sampling for parameters $\lambda_{j}$

In this section we propose a novel rejection sampling algorithm to sample
shrinkage parameters $\lambda_{j}$, $j=1,\dots,p$, from the full-conditional
distribution without knowing the normalizing constant. The rejection sampling
allows to draw a new value $x_{\star}$ from a density $f$ by sampling it from
a proposal distribution $g$ and accepting it with a probability proportional
to the ratio $r(x_{\star})=f(x_{\star})/g(x_{\star})$. In particular, $g$ must
be chosen such that the support of $f$ is a subset of the support of $g$. Let

$k=\text{sup}_{x}\frac{f(x)}{g(x)}<\infty$

and accept $x_{\star}$ with probability $a=f(x_{\star})/(kg(x_{\star}))$. It
can be shown that the acceptance probability of the algorithm is $1/k$.
Therefore, the goal is to find a proposal distribution $g$ such that $k$ is
small. The rejection sampling works also if the normalizing constant of $f$ is
unknown, as long as a sampling method for $g$ is available.

Consider the following density

$f(x)=c_{f}x^{-1}e^{-\psi/x^{2}-\alpha^{2}x^{2}+\beta
x}\cdot\mathbb{I}_{(x>0)},$

where $\psi,\alpha^{2}>0$, $\beta\in\mathbb{R}$ and $c_{f}$ is the unknown
normalizing constant. New values are sampled from the proposal distribution
$g(x)\sim\mathcal{G}_{3p}\left(\gamma,\alpha,\beta\right)$, with density
$g(x)=c_{g}x^{\gamma}e^{-\alpha^{2}x^{2}+\beta x}$, where
$\gamma\in\mathbb{N}^{+}$. We set parameters $\alpha$ and $\beta$ equal in
both $f$ and $g$ in order for the ratio $r(x)$ to be analytically tractable
and easily maximized. Eventually, the acceptance probability is optimized
through the choice of parameter $\gamma$. This choice of $\alpha$ and $\beta$
yields

$r(x)=\frac{c_{f}}{c_{g}}x^{-\gamma-1}e^{-\psi/x^{2}}\cdot\mathbb{I}_{(x>0)},$

which has one positive maximum at $\dot{x}=\sqrt{2\psi/(\gamma+1)}$.
Theoretical properties and details for sampling from a $\mathcal{G}_{3p}$
distribution are shown in Busatto and Stingo, (2023).

The best theoretical way to choose the parameter $\gamma$ is to differentiate
$k=r(\dot{x})c_{f}/c_{g}$ with respect to $\gamma$ and find the optimal value.
However, this can only be done with iterative methods, which negatively
affects the computational efficiency of the method. An alternative solution is
to set parameter $\gamma$ such that distributions $f$ and $g$ have the maximum
at the same value $x_{max}$. That is, the maximum of $f$ is computed by
solving a quartic equation and the parameter $\gamma$ is estimated as

$\gamma=x_{max}\left(2\alpha^{2}x_{max}-\beta\right)>-1.$

Since $\gamma\in\mathbb{N}_{0}^{+}$, the closest non-negative integer is
chosen. An example of the method is shown in Figure 4 of Supporting
Information, where the normalizing constant $c_{f}$ is computed by numerical
integration.

For some settings of the parameters $\psi$, $\alpha$ and $\beta$ the
acceptance probability of the proposed algorithm decreases toward $0$,
affecting the efficiency of the algorithm. For this reason, and for avoiding
the explicit inversion of covariance matrix
$\boldsymbol{\Sigma}_{\beta}^{\star}$, in the next section we present a
Variational Bayes algorithm which overcomes these problems.

Figure 4: Comparison of $f(x)$ and $g(x)$ with $\psi=2$, $\alpha^{2}=2.25$ and
$\beta=-2$. Estimated parameter $\gamma=5$; the acceptance probability is
$a=0.65$.

##  Appendix C Web Appendix C

###  C.1 Variational Inference

Variational inference (VI) is a deterministic optimization method to
approximate the target density
$\pi\left(\boldsymbol{\theta}\mid\mathbf{y},\mathbf{X},\mathbf{Z}\right)$,
with
$\boldsymbol{\theta}=\left(\boldsymbol{\beta},\lambda_{0}^{2},\boldsymbol{\lambda},\psi_{0},\boldsymbol{\varphi},\boldsymbol{\gamma},\boldsymbol{\kappa},\tau^{2},\zeta,\sigma^{2}\right)$,
with another variational distribution $q\left(\boldsymbol{\theta}\right)$ and
reduces Bayesian inference to an optimization problem (Salimans et al.,, 2015;
Lee,, 2022). The goal is to find $q(\boldsymbol{\theta})$ that minimizes the
Kullback-Leibler divergence (KL) between the target density and the
variational distribution. The minimization problem in is eventually reduced to
the maximization of the variational lower bound, defined as
$\mathcal{L}=\mathbb{E}_{q}\left[\log\pi\left(\boldsymbol{\theta},\mathbf{y},\mathbf{X},\mathbf{Z}\right)\right]-\mathbb{E}_{q}\left[\log
q\left(\boldsymbol{\theta}\right)\right]$.

A common factorization for $q(\boldsymbol{\theta})$ is the so-called mean-
field Variational approximation (Jordan et al.,, 1999; Beal,, 2003), which is
a compromise between computational tractability and accuracy of the
performances. The variational family $q(\boldsymbol{\theta})$ is assumed to be
the product of independent marginal variational factors
$q_{k}\left(\theta_{k}\right)$, $k=1,\dots,K$. We rely on the Coordinate
Ascent Variational Inference algorithm (CAVI) (Bishop and Nasrabadi,, 2006;
Blei et al.,, 2017) to solve the optimization problem above. Until convergence
of the lower bound $\mathcal{L}$, the CAVI algorithm iteratively updates the
parameters of the variational factors $q_{k}\left(\theta_{k}\right)$,
$k=1,\dots,K$, based on prior distributions’ hyperparameters and the current
expectation of factor $q_{-k}\left(\theta_{-k}\right)$, considered fixed
(Lee,, 2022). This way the model is able to account for non-linear
dependencies among the parameters. Under the mean field approximation, where
the components are assumed to be independent, the optimal solution is given by

$q^{\star}\left(\theta_{k}\right)\propto\exp\left\\{\mathbb{E}_{q_{-k}}\left[\log\pi\left(\theta_{k}\mid\theta_{-k},\mathbf{y},\mathbf{X},\mathbf{Z}\right)\right]\right\\}.$

While the assumption of independence between factors is particularly strict,
the CAVI algorithm provides a flexible approach and is ensured to converge to
a local optimum (Blei et al.,, 2017). Note that, when working with exponential
families in a conjugated framework, variational factor $q(\theta_{k})$ has the
same kernel of the tractable distribution
$\pi\left(\theta_{k}\mid\theta_{-k},\mathbf{y},\mathbf{X},\mathbf{Z}\right)$.

The variational lower bound is computed as

$\displaystyle\mathcal{L}=$
$\displaystyle\;\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\theta},\mathbf{y},\mathbf{X},\mathbf{Z}\right)\right]-\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\theta}\right)\right]$ $\displaystyle=$
$\displaystyle\;\mathbb{E}_{q^{\star}}\left[\log
p\left(\mathbf{y}\mid\mathbf{X},\boldsymbol{\beta},\sigma^{2}\right)\right]+\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\beta}\mid\sigma^{2},\tau^{2},\lambda_{0}^{2},\boldsymbol{\lambda}\right)\right]+\mathbb{E}_{q^{\star}}\left[\log\pi\left(\lambda_{0}^{2}\mid\psi_{0}\right)\right]+$
$\displaystyle\;\mathbb{E}_{q^{\star}}\left[\log\pi\left(\psi_{0}\right)\right]+\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\lambda}\mid\mathbf{Z},\boldsymbol{\gamma},\boldsymbol{\varphi}^{2}\right)\right]+\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\varphi}^{2}\right)\right]+\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\gamma}\mid\boldsymbol{\kappa}^{2}\right)\right]+$
$\displaystyle\;\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\kappa}^{2}\right)\right]+\mathbb{E}_{q^{\star}}\left[\log\pi\left(\tau^{2}\mid\zeta\right)\right]+\mathbb{E}_{q^{\star}}\left[\log\pi\left(\zeta\right)\right]+\mathbb{E}_{q^{\star}}\left[\log\pi\left(\sigma^{2}\right)\right]-$
$\displaystyle\;\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\beta}\right)\right]-\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\lambda_{0}^{2}\right)\right]-\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\psi_{0}\right)\right]-\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\lambda}\right)\right]-\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\varphi}^{2}\right)\right]-$
$\displaystyle\;\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\gamma}\right)\right]-\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\kappa}^{2}\right)\right]-\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\tau^{2}\right)\right]-\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\zeta\right)\right]-\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\sigma^{2}\right)\right]$ $\displaystyle\propto$
$\displaystyle\;\frac{1}{2}\left(\log|\boldsymbol{\Sigma}_{\beta}^{\star}|+\log|\boldsymbol{\Sigma}_{\gamma}^{\star}|\right)+\frac{p+1}{2}\mathbb{E}_{\sigma^{2}}\left[\log\sigma^{2}\right]-\left(v+\frac{n+p+1}{2}\right)\log
l^{\star}+$ $\displaystyle\;\sum_{j=1}^{p}\left(\log s_{j}-\log
k_{j}+a_{j}^{\star}\mathbb{E}_{\lambda}\left[\lambda_{j}^{-2}\right]+b_{j}^{\star}\mathbb{E}_{\lambda}\left[\lambda_{j}^{2}\right]-c_{j}^{\star}\mathbb{E}_{\lambda}\left[\lambda_{j}\right]-\log
d_{j}^{\star}\right)-$ $\displaystyle\;\sum_{d=1}^{D}e_{d}^{\star}\log
f_{d}^{\star}-\left(\frac{p}{2}+1\right)\log g^{\star}-\log h^{\star}-\log
a_{0}^{\star}-\log k_{0}^{\star},$ (6)

where
$k_{j}=1-\mathbb{P}\left(\mathcal{N}\left(\mathbf{z}_{j}^{\intercal}\boldsymbol{\gamma},s_{0}^{2}\varphi_{j}^{2}\right)<0\right)$
and $s_{j}$ is the unknown normalizing constant of parameters $\lambda_{j}$ in
(4) of the main article. Each component is derived in Section D.1 of
Supporting Information.

As opposed to the Gibbs sampler, the Variational algorithm does not require
the explicit inversion of the $p\times p$ matrix
$\boldsymbol{\Sigma}_{\beta}^{\star}$, as only the quantities
$\boldsymbol{\Sigma}_{\beta}^{\star}\mathbf{X}^{\intercal}\mathbf{y}$,
$\text{diag}\left(\boldsymbol{\Sigma}_{\beta}^{\star}\right)$ and
$\mathbf{X}\boldsymbol{\Sigma}_{\beta}^{\star}\mathbf{X}^{\intercal}$ are
needed. These can be efficiently evaluated with complexity
$\mathcal{O}\left(n^{2}p\right)$ following the work of Münch et al., (2019).
The details are shown in Section E.1 of Supporting Information. The bottleneck
of the algorithm is the computation of the terms $\log s_{j}$,
$\mathbb{E}_{\lambda}\left[\lambda_{j}\right]$,
$\mathbb{E}_{\lambda}\left[\lambda_{j}^{2}\right]$ and
$\mathbb{E}_{\lambda}\left[\lambda_{j}^{-2}\right]$, since no closed-form is
available. These can be evaluated by numerical integration, for example with
the Gaussian quadrature rule or its adaptive variation called Gauss-Kronrod
quadrature formula. Depending on the posterior parameters $a_{j}^{\star}$,
$b_{j}^{\star}$ and $c_{j}^{\star}$, this strategy is prone to numerical
instability. To avoid overflow problem, we evaluate these integrals with the
following steps:

  1. Step 0: assume we have to evaluate $\int_{0}^{\infty}f(x)dx$, where $f(x)\propto x^{\nu}\exp\left\\{-dx^{-2}-bx^{2}+cx\right\\}\cdot\mathbb{I}_{\left(\lambda_{j}>0\right)}$, where $\nu=\left\\{-3,-1,0,1\right\\}$;

  2. Step 1: compute the maximum $\dot{x}$ by solving a quartic equation and evaluate $f(\dot{x})$ on the log-scale;

  3. Step 2: evaluate $i_{0}=\int_{0}^{\infty}\exp\left\\{\log f(x)-\log f(\dot{x})\right\\}dx$ numerically;

  4. Step 3: return $\exp\left\\{\log i_{0}+\log f(\dot{x})\right\\}$.

###  C.2 Variational Inference in probit regression

Let $\mathbf{y}\in\left\\{0,1\right\\}^{n}$ be a binary vector and
$\mathbf{X}\in\mathbb{R}^{n\times(p+1)}$ be the design matrix. Following
Albert and Chib, (1993) we introduce the latent variable
$\mathbf{w}\in\mathbb{R}^{n}$ to reach a conjugated framework for updating
regression coefficient vector $\boldsymbol{\beta}$. The assumptions of the
probit model are

$\displaystyle y_{i}\mid\mathbf{w}=$
$\displaystyle\;\begin{cases}1\quad\text{if}\;w_{i}>0\\\
0\quad\text{if}\;w_{i}\leq 0,\end{cases}\,$ $\displaystyle
w_{i}\mid\mathbf{X},\boldsymbol{\beta}\sim$
$\displaystyle\;\mathcal{N}\left(\mathbf{x}_{i}^{\intercal}\boldsymbol{\beta},1\right),\quad
i=1,\dots,n,$ $\displaystyle\beta_{0}\mid\tau^{2},\lambda_{0}^{2}\sim$
$\displaystyle\;\mathcal{N}\left(0,\tau^{2}\lambda_{0}^{2}\right)$
$\displaystyle\beta_{j}\mid\tau,\lambda_{j}\sim$
$\displaystyle\;\mathcal{N}\left(0,\tau^{2}\lambda_{j}^{2}\right),\quad
j=1,\dots,p.$

The prior distributions for parameters
$\boldsymbol{\lambda},\boldsymbol{\gamma}$ and $\tau$ are the same of Section
2 and the algorithm follows the same steps of Section 3 of the main article.
The introduction of the latent variable $\mathbf{w}$ allows to formulate the
problem as a normal regression model on the latent variables $w_{i}$. The
normal prior distribution for the regression parameters vector
$\boldsymbol{\beta}$ is conjugated to the the normal distribution of
$\mathbf{w}$. The joint posterior distribution of the model is

$\pi\left(\boldsymbol{\theta},\mathbf{w},\mathbf{y},\mathbf{X},\mathbf{Z}\right)\propto
p\left(\mathbf{y}\mid\mathbf{w}\right)\cdot\pi\left(\mathbf{w}\mid\mathbf{X},\boldsymbol{\theta}\right)\cdot\pi\left(\boldsymbol{\theta}\mid\mathbf{Z}\right),$

where
$\boldsymbol{\theta}=\left(\boldsymbol{\beta},\lambda_{0}^{2},\boldsymbol{\lambda},\psi_{0},\boldsymbol{\varphi},\boldsymbol{\gamma},\boldsymbol{\kappa},\tau^{2},\zeta\right)$.
Under the mean field approximation, the variational factors for parameters
$\mathbf{w}$ and $\boldsymbol{\beta}$ are updated as

$\displaystyle q^{\star}\left(w_{i}\right)=$
$\displaystyle\;\begin{cases}\mathcal{N}_{+}\left(\mu_{i}^{\star},1\right)\quad\text{if}\;y_{i}=0\\\
\mathcal{N}_{-}\left(\mu_{i}^{\star},1\right)\quad\text{if}\;y_{i}=0,\end{cases}\,$
$\displaystyle q^{\star}\left(\boldsymbol{\beta}\right)=$
$\displaystyle\;\mathcal{N}_{p}\left(\boldsymbol{\mu}_{\beta}^{\star},\boldsymbol{\Sigma}_{\beta}^{\star}\right),$
$\displaystyle\;\boldsymbol{\mu}_{\beta}^{\star}=\boldsymbol{\Sigma}_{\beta}^{\star}\mathbf{X}^{\intercal}\mathbb{E}_{w}\left[\mathbf{w}\right],\quad\boldsymbol{\Sigma}_{\beta}^{\star}=\left(\mathbf{X}^{\intercal}\mathbf{X}+\mathbb{E}_{\lambda_{0}^{2}\cdot\lambda\cdot\tau^{2}}\left[\tau^{-2}\boldsymbol{\Lambda}^{-2}\right]\right)^{-1},$

where
$\mu_{i}^{\star}=\mathbf{x}_{i}^{\intercal}\mathbb{E}_{\beta}\left[\boldsymbol{\beta}\right]$
and $\mathcal{N}_{+}$ and $\mathcal{N}_{-}$ denote a normal distribution left
and right truncated at 0, respectively. The expectation of latent variables
$w_{i}$ is

$\mathbb{E}_{w}\left[w_{i}\right]=\begin{cases}\mu_{i}^{\star}+\frac{\phi_{i}}{1-\Phi_{i}}\quad\text{if}\;y_{i}=0\\\
\mu_{i}^{\star}-\frac{\phi_{i}}{\Phi_{i}}\quad\quad\text{if}\;y_{i}=0,\end{cases}\,$

where $\phi_{i}=\phi\left(-\mu_{i}^{\star}\right)$ and
$\Phi_{i}=\Phi\left(-\mu_{i}^{\star}\right)$ are the normal density and the
normal cumulative density functions, respectively. The variational lower bound
is computed as

$\displaystyle\mathcal{L}=$
$\displaystyle\;\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\theta},\mathbf{w},\mathbf{y},\mathbf{X},\mathbf{Z}\right)\right]-\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\theta},\mathbf{w}\right)\right]$ $\displaystyle=$
$\displaystyle\;\mathbb{E}_{q^{\star}}\left[\log
p\left(\mathbf{y}\mid\mathbf{w}\right)\right]+\mathbb{E}_{q^{\star}}\left[\log\pi\left(\mathbf{w}\mid\mathbf{X},\boldsymbol{\beta}\right)\right]+\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\beta}\mid\tau^{2},\lambda_{0}^{2},\boldsymbol{\lambda}\right)\right]+\mathbb{E}_{q^{\star}}\left[\log\pi\left(\lambda_{0}^{2}\mid\psi_{0}\right)\right]+$
$\displaystyle\;\mathbb{E}_{q^{\star}}\left[\log\pi\left(\psi_{0}\right)\right]+\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\lambda}\mid\mathbf{Z},\boldsymbol{\gamma},\boldsymbol{\varphi}^{2}\right)\right]+\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\varphi}^{2}\right)\right]+\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\gamma}\mid\boldsymbol{\kappa}^{2}\right)\right]+$
$\displaystyle\;\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\kappa}^{2}\right)\right]+\mathbb{E}_{q^{\star}}\left[\log\pi\left(\tau^{2}\mid\zeta\right)\right]+\mathbb{E}_{q^{\star}}\left[\log\pi\left(\zeta\right)\right]-\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\mathbf{w}\right)\right]-\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\beta}\right)\right]-$
$\displaystyle\;\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\lambda_{0}^{2}\right)\right]-\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\psi_{0}\right)\right]-\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\lambda}\right)\right]-\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\varphi}^{2}\right)\right]-\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\gamma}\right)\right]-$
$\displaystyle\;\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\kappa}^{2}\right)\right]-\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\tau^{2}\right)\right]-\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\zeta\right)\right]$ $\displaystyle\propto$
$\displaystyle\;\sum_{i=1}^{n}\left(y_{i}\log\left(1-\Phi_{i}\right)+(1-y_{i})\log\left(\Phi_{i}\right)\right)+\frac{1}{2}\left(\log|\boldsymbol{\Sigma}_{\beta}^{\star}|+\log|\boldsymbol{\Sigma}_{\gamma}^{\star}|\right)+$
$\displaystyle\;\sum_{j=1}^{p}\left(\log s_{j}-\log
k_{j}\right)+\sum_{j=1}^{p}\left(a_{j}^{\star}\mathbb{E}_{\lambda}\left[\lambda_{j}^{-2}\right]+b_{j}^{\star}\mathbb{E}_{\lambda}\left[\lambda_{j}^{2}\right]-c_{j}^{\star}\mathbb{E}_{\lambda}\left[\lambda_{j}\right]-\log
d_{j}^{\star}\right)-$ $\displaystyle\;\sum_{d=1}^{D}e_{d}^{\star}\log
f_{d}^{\star}-\left(\frac{p}{2}+1\right)\log g^{\star}-\log h^{\star}-\log
a_{0}^{\star}-\log k_{0}^{\star},$ (7)

where all the quantities are defined in Section D.2 of Supporting Information.
The method requires the computation of
$\boldsymbol{\Sigma}_{\beta}^{\star}\mathbf{X}^{\intercal}\mathbb{E}_{w}\left[w\right]$
and $\text{diag}\left(\boldsymbol{\Sigma}_{\beta}^{\star}\right)$, which can
efficiently be evaluated by applying the strategy in Section E.1 of Supporting
Information.

##  Appendix D Web Appendix D

Here we report the components for the evaluation of the lower bound
$\mathcal{L}$ in (6) and (7) of Supporting Information. The results below rely
on the property
$\mathbb{E}_{x}\left[\mathbf{x}^{\intercal}\mathbf{A}\mathbf{x}\right]=\mathbb{E}_{x}(\mathbf{x})^{\intercal}\mathbf{A}\mathbb{E}_{x}(\mathbf{x})+tr\left(\mathbf{A}\mathbf{V}\right)$,
where $\mathbf{V}$ is the correlation matrix of generical random variable
$\mathbf{x}$.

###  D.1 Linear regression

The variational lower bound $\mathcal{L}$ in (6) is computed with the
following components:

$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
p\left(\mathbf{y}\mid\boldsymbol{\beta},\sigma^{2},\boldsymbol{\lambda}\right)\right]\propto-\frac{n}{2}\mathbb{E}_{\sigma^{2}}\left[\log\sigma^{2}\right]-\frac{1}{2}\mathbb{E}_{\beta\cdot\sigma^{2}}\left[\frac{\|\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\|_{2}^{2}}{\sigma^{2}}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\beta}\mid\sigma^{2},\tau^{2},\lambda_{0}^{2},\boldsymbol{\lambda}\right)\right]\propto-\frac{p+1}{2}\mathbb{E}_{\sigma^{2}}\left[\log\sigma^{2}\right]-\frac{p+1}{2}\mathbb{E}_{\tau^{2}}\left[\log\tau^{2}\right]-\frac{1}{2}\mathbb{E}_{\lambda_{0}^{2}}\left[\log\lambda_{0}^{2}\right]-$
$\displaystyle\hskip
150.79968pt\sum_{j=1}^{p}\mathbb{E}_{\lambda}\left[\log\lambda_{j}\right]-\frac{1}{2}\sum_{j=0}^{p}\mathbb{E}_{\beta_{0}\cdot\beta\cdot\lambda_{0}^{2}\cdot\lambda\cdot\tau^{2}\cdot\sigma^{2}}\left[\frac{\beta_{j}^{2}}{\sigma^{2}\tau^{2}\lambda_{j}^{2}}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log\pi\left(\lambda_{0}^{2}\mid\psi_{0}\right)\right]\propto-\frac{1}{2}\mathbb{E}_{\psi_{0}}\left[\log\psi_{0}\right]-\frac{3}{2}\mathbb{E}_{\lambda_{0}^{2}}\left[\log\lambda_{0}^{2}\right]-\mathbb{E}_{\lambda_{0}^{2}\cdot\psi_{0}}\left[\frac{1}{\psi_{0}\lambda_{0}^{2}}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log\pi\left(\psi_{0}\right)\right]\propto-\frac{3}{2}\mathbb{E}_{\psi_{0}}\left[\log\psi_{0}\right]-\mathbb{E}_{\psi_{0}}\left[\psi_{0}^{-1}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\lambda}\mid\boldsymbol{\gamma},\boldsymbol{\varphi}^{2}\right)\right]\propto\sum_{j=1}^{p}\Bigg{(}-\log
k_{j}-\frac{1}{2}\mathbb{E}_{\varphi^{2}}\left[\log\varphi_{j}^{2}\right]-\frac{1}{2s_{0}^{2}}\mathbb{E}_{\lambda\cdot\gamma\cdot\varphi^{2}}\left[\frac{\left(\lambda_{j}-\mathbf{z}_{j}^{\intercal}\boldsymbol{\gamma}\right)^{2}}{\varphi_{j}^{2}}\right]\Bigg{)}$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\varphi}^{2}\right)\right]\propto-\frac{1}{2}\sum_{j=1}^{p}\Bigg{(}3\mathbb{E}_{\varphi^{2}}\left[\log\varphi_{j}^{2}\right]+\mathbb{E}_{\varphi^{2}}\left[\varphi_{j}^{-2}\right]\Bigg{)}$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\gamma}\mid\boldsymbol{\kappa}^{2}\right)\right]\propto-\frac{1}{2}\sum_{d=1}^{D}\Bigg{(}m_{d}\mathbb{E}_{\kappa^{2}}\left[\log\kappa_{d}^{2}\right]+\mathbb{E}_{\gamma\cdot\kappa^{2}}\left[\frac{\boldsymbol{\gamma}_{d}^{\intercal}\boldsymbol{\gamma}_{d}}{\kappa_{d}^{2}}\right]\Bigg{)}$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\kappa}^{2}\right)\right]\propto-\sum_{d=1}^{D}\Bigg{(}(a_{d}+1)\mathbb{E}_{\kappa^{2}}\left[\log\kappa_{d}^{2}\right]+b_{d}\mathbb{E}_{\kappa^{2}}\left[\kappa_{d}^{-2}\right]\Bigg{)}$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log\pi\left(\tau^{2}\mid\zeta\right)\right]\propto-\frac{1}{2}\mathbb{E}_{\zeta}\left[\log\zeta\right]-\frac{3}{2}\mathbb{E}_{\tau^{2}}\left[\log\tau^{2}\right]-\mathbb{E}_{\tau^{2}\cdot\zeta}\left[\frac{1}{\zeta\tau^{2}}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log\pi\left(\zeta\right)\right]\propto-\frac{3}{2}\mathbb{E}_{\zeta}\left[\log\zeta\right]-\mathbb{E}_{\zeta}\left[\zeta^{-1}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log\pi\left(\sigma^{2}\right)\right]\propto-\left(v+1\right)\mathbb{E}_{\sigma^{2}}\left[\log\sigma^{2}\right]-q\mathbb{E}_{\sigma^{2}}\left[\sigma^{-2}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\beta}\right)\right]\propto-\frac{p+1}{2}\mathbb{E}_{\sigma^{2}}\left[\log\sigma^{2}\right]-\frac{1}{2}\log|\boldsymbol{\Sigma}_{\beta}^{\star}|-\frac{p+1}{2}$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\lambda_{0}^{2}\right)\right]\propto\log\left(a_{0}^{\star}\right)-2\mathbb{E}_{\lambda_{0}^{2}}\left[\log\lambda_{0}^{2}\right]-\mathbb{E}_{\lambda_{0}^{2}\cdot\psi_{0}}\left[\frac{1}{\psi_{0}\lambda_{0}^{2}}\right]-\frac{1}{2}\mathbb{E}_{\beta_{0}\cdot\lambda_{0}^{2}\cdot\tau^{2}\cdot\sigma^{2}}\left[\frac{\beta_{0}^{2}}{\sigma^{2}\tau^{2}\lambda_{0}^{2}}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\psi_{0}\right)\right]\propto\log
k_{0}^{\star}-2\mathbb{E}_{\psi_{0}}\left[\log\psi_{0}\right]-\mathbb{E}_{\psi_{0}}\left[\psi_{0}^{-1}\right]-\mathbb{E}_{\lambda_{0}^{2}\cdot\psi_{0}}\left[\frac{1}{\psi_{0}\lambda_{0}^{2}}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\lambda}\right)\right]\propto\sum_{j=1}^{p}\Bigg{(}-\log
s_{j}-\mathbb{E}_{\lambda}\left[\log\lambda_{j}\right]-\frac{1}{2}\mathbb{E}_{\beta\cdot\lambda\cdot\sigma^{2}\cdot\tau^{2}}\left[\frac{\beta_{j}^{2}}{\sigma^{2}\tau^{2}\lambda_{j}^{2}}\right]-$
$\displaystyle\hskip
85.35826pt\frac{1}{2s_{0}^{2}}\mathbb{E}_{\lambda\cdot\varphi^{2}}\left[\frac{\lambda_{j}^{2}}{\varphi_{j}^{2}}\right]+\frac{1}{s_{0}^{2}}\mathbf{z}_{j}^{\intercal}\mathbb{E}_{\lambda\cdot\gamma\cdot\varphi^{2}}\left[\frac{\boldsymbol{\gamma}\lambda_{j}}{\varphi_{j}^{2}}\right]\Bigg{)}$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\varphi}^{2}\right)\right]\propto\sum_{j=1}^{p}\Bigg{(}\log
d_{j}^{\star}-2\mathbb{E}_{\varphi^{2}}\left[\log\varphi_{j}^{2}\right]-\frac{1}{2}\mathbb{E}_{\varphi^{2}}\left[\varphi_{j}^{-2}\right]-\frac{1}{2s_{0}^{2}}\mathbb{E}_{\lambda\cdot\gamma\cdot\varphi^{2}}\left[\frac{\left(\lambda_{j}-\mathbf{z}_{j}^{\intercal}\boldsymbol{\gamma}\right)^{2}}{\varphi_{j}^{2}}\right]\Bigg{)}$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\gamma}\right)\right]\propto-\frac{1}{2}\log|\boldsymbol{\Sigma}_{\gamma}^{\star}|-\frac{M}{2s_{0}^{2}}$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\kappa^{2}\right)\right]\propto\sum_{d=1}^{D}\Bigg{(}\left(a_{d}+\frac{m_{d}}{2}\right)\log
f_{d}^{\star}-b_{d}\mathbb{E}_{\kappa^{2}}\left[\kappa_{d}^{-2}\right]-\left(a_{d}+\frac{m_{d}}{2}+1\right)\mathbb{E}_{\kappa^{2}}\left[\log\kappa_{d}^{2}\right]-$
$\displaystyle\hskip
105.2751pt\frac{1}{2}\mathbb{E}_{\gamma\cdot\kappa^{2}}\left[\frac{\boldsymbol{\gamma}_{d}^{\intercal}\boldsymbol{\gamma}_{d}}{\kappa_{d}^{2}}\right]\Bigg{)}$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\tau^{2}\right)\right]\propto\left(\frac{p}{2}+1\right)\log
g^{\star}-\left(\frac{p}{2}+2\right)\mathbb{E}_{\tau^{2}}\left[\log\tau^{2}\right]-\mathbb{E}_{\tau^{2}\cdot\zeta}\left[\frac{1}{\zeta\tau^{2}}\right]-$
$\displaystyle\hskip
91.04872pt\frac{1}{2}\sum_{j=0}^{p}\mathbb{E}_{\beta_{0}\cdot\beta\cdot\lambda_{0}^{2}\cdot\lambda\cdot\tau^{2}\cdot\sigma^{2}}\left[\frac{\beta_{j}^{2}}{\sigma^{2}\tau^{2}\lambda_{j}^{2}}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\zeta\right)\right]\propto\log
h^{\star}-2\mathbb{E}_{\zeta}\left[\log\zeta\right]-\mathbb{E}_{\zeta}\left[\zeta^{-1}\right]-\mathbb{E}_{\tau^{2}\cdot\zeta}\left[\frac{1}{\zeta\tau^{2}}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\sigma^{2}\right)\right]\propto\left(v+\frac{n+p+1}{2}\right)\log
l^{\star}-q\mathbb{E}_{\sigma^{2}}\left[\sigma^{-2}\right]-\left(v+\frac{n+p+3}{2}\right)\mathbb{E}_{\sigma^{2}}\left[\log\sigma^{2}\right]-$
$\displaystyle\hskip
99.58464pt\frac{1}{2}\mathbb{E}_{\beta\cdot\sigma^{2}}\left[\frac{\|\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\|_{2}^{2}}{\sigma^{2}}\right]-\frac{1}{2}\sum_{j=0}^{p}\mathbb{E}_{\beta_{0}\cdot\beta\cdot\lambda_{0}^{2}\cdot\lambda\cdot\sigma^{2}\cdot\tau^{2}}\left[\frac{\beta_{j}^{2}}{\sigma^{2}\tau^{2}\lambda_{j}^{2}}\right]$

###  D.2 Probit regression

The variational lower bound $\mathcal{L}$ in (7) of Supporting Information is
computed with the following components:

$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
p\left(\mathbf{y}\mid\mathbf{w}\right)\right]\propto\sum_{i=1}^{n}\Bigg{(}\mathbb{E}_{w}\left[y_{i}\log(\mathbb{I}_{\left(w_{i}>0\right)}\right]+\mathbb{E}_{w}\left[(1-y_{i})\log(\mathbb{I}_{\left(w_{i}<0\right)}\right]\Bigg{)}=0$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
p\left(\mathbf{w}\mid\boldsymbol{\beta}\right)\right]\propto-\frac{1}{2}\sum_{i=1}^{n}\mathbb{E}_{w}\left[\|w_{i}-\mathbf{x}_{i}^{\intercal}\boldsymbol{\beta}\|_{2}^{2}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\beta}\mid\sigma^{2},\tau^{2},\lambda_{0}^{2},\boldsymbol{\lambda}\right)\right]\propto-\frac{p+1}{2}\mathbb{E}_{\tau^{2}}\left[\log\tau^{2}\right]-\frac{1}{2}\mathbb{E}_{\lambda_{0}^{2}}\left[\log\lambda_{0}^{2}\right]-\sum_{j=1}^{p}\mathbb{E}_{\lambda}\left[\log\lambda_{j}\right]-$
$\displaystyle\hskip
159.3356pt\frac{1}{2}\sum_{j=0}^{p}\mathbb{E}_{\beta_{0}\cdot\beta\cdot\lambda_{0}^{2}\cdot\lambda\cdot\tau^{2}}\left[\frac{\beta_{j}^{2}}{\tau^{2}\lambda_{j}^{2}}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log\pi\left(\lambda_{0}^{2}\mid\psi_{0}\right)\right]\propto-\frac{1}{2}\mathbb{E}_{\psi_{0}}\left[\log\psi_{0}\right]-\frac{3}{2}\mathbb{E}_{\lambda_{0}^{2}}\left[\log\lambda_{0}^{2}\right]-\mathbb{E}_{\lambda_{0}^{2}\cdot\psi_{0}}\left[\frac{1}{\psi_{0}\lambda_{0}^{2}}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log\pi\left(\psi_{0}\right)\right]\propto-\frac{3}{2}\mathbb{E}_{\psi_{0}}\left[\log\psi_{0}\right]-\mathbb{E}_{\psi_{0}}\left[\psi_{0}^{-1}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\lambda}\mid\boldsymbol{\gamma},\boldsymbol{\varphi}^{2}\right)\right]\propto\sum_{j=1}^{p}\Bigg{(}-\log
k_{j}-\frac{1}{2}\mathbb{E}_{\varphi^{2}}\left[\log\varphi_{j}^{2}\right]-\frac{1}{2s_{0}^{2}}\mathbb{E}_{\lambda\cdot\gamma\cdot\varphi^{2}}\left[\frac{\left(\lambda_{j}-\mathbf{z}_{j}^{\intercal}\boldsymbol{\gamma}\right)^{2}}{\varphi_{j}^{2}}\right]\Bigg{)}$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\varphi}^{2}\right)\right]\propto-\frac{1}{2}\sum_{j=1}^{p}\Bigg{(}3\mathbb{E}_{\varphi^{2}}\left[\log\varphi_{j}^{2}\right]+\mathbb{E}_{\varphi^{2}}\left[\varphi_{j}^{-2}\right]\Bigg{)}$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\gamma}\mid\boldsymbol{\kappa}^{2}\right)\right]\propto-\frac{1}{2}\sum_{d=1}^{D}\Bigg{(}m_{d}\mathbb{E}_{\kappa}\left[\log\kappa_{d}^{2}\right]+\mathbb{E}_{\gamma\cdot\kappa}\left[\frac{\boldsymbol{\gamma}_{d}^{\intercal}\boldsymbol{\gamma}_{d}}{\kappa_{d}^{2}}\right]\Bigg{)}$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log\pi\left(\boldsymbol{\kappa}^{2}\right)\right]\propto\sum_{d=1}^{D}\Bigg{(}-(a_{d}+1)\mathbb{E}_{\kappa^{2}}\left[\log\kappa_{d}^{2}\right]-b_{d}\mathbb{E}_{\kappa^{2}}\left[\kappa_{d}^{-2}\right]\Bigg{)}$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log\pi\left(\tau^{2}\mid\zeta\right)\right]\propto-\frac{1}{2}\mathbb{E}_{\zeta}\left[\log\zeta\right]-\frac{3}{2}\mathbb{E}_{\tau^{2}}\left[\log\tau^{2}\right]-\mathbb{E}_{\tau^{2}\cdot\zeta}\left[\frac{1}{\zeta\tau^{2}}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log\pi\left(\zeta\right)\right]\propto-\frac{3}{2}\mathbb{E}_{\zeta}\left[\log\zeta\right]-\mathbb{E}_{\zeta}\left[\zeta^{-1}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\mathbf{w}\right)\right]\propto\sum_{i=1}^{n}\Bigg{(}-\frac{1}{2}\mathbb{E}_{w}\left[\|w_{i}-\mathbf{x}_{i}^{\intercal}\boldsymbol{\beta}\|_{2}^{2}\right]-y_{i}\log\left(1-\Phi_{i}\right)-\left(1-y_{i}\right)\log\left(\Phi_{i}\right)\Bigg{)}$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\beta}\right)\right]\propto-\frac{1}{2}\log|\boldsymbol{\Sigma}_{\beta}^{\star}|-\frac{p+1}{2}$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\lambda_{0}^{2}\right)\right]\propto\log\left(a_{0}^{\star}\right)-2\mathbb{E}_{\lambda_{0}^{2}}\left[\log\lambda_{0}^{2}\right]-\mathbb{E}_{\lambda_{0}^{2}\cdot\psi_{0}}\left[\frac{1}{\psi_{0}\lambda_{0}^{2}}\right]-\frac{1}{2}\mathbb{E}_{\beta_{0}\cdot\tau^{2}\cdot\lambda_{0}^{2}}\left[\frac{\beta_{0}^{2}}{\tau^{2}\lambda_{0}^{2}}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\psi_{0}\right)\right]\propto\log
k_{0}^{\star}-2\mathbb{E}_{\psi_{0}}\left[\log\psi_{0}\right]-\mathbb{E}_{\psi_{0}}\left[\psi_{0}^{-1}\right]-\mathbb{E}_{\lambda_{0}^{2}\cdot\psi_{0}}\left[\frac{1}{\psi_{0}\lambda_{0}^{2}}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\lambda}\right)\right]\propto\sum_{j=1}^{p}\Bigg{(}-\log
s_{j}-\mathbb{E}_{\lambda}\left[\log\lambda_{j}\right]-\frac{1}{2}\mathbb{E}_{\beta\cdot\tau^{2}\cdot\lambda}\left[\frac{\beta_{j}^{2}}{\tau^{2}\lambda_{j}^{2}}\right]-\frac{1}{2s_{0}^{2}}\mathbb{E}_{\lambda\cdot\varphi^{2}}\left[\frac{\lambda_{j}^{2}}{\varphi_{j}^{2}}\right]+$
$\displaystyle\hskip
88.2037pt\frac{1}{s_{0}^{2}}\mathbf{z}_{j}^{\intercal}\mathbb{E}_{\lambda\cdot\gamma\cdot\varphi^{2}}\left[\frac{\boldsymbol{\gamma}\lambda_{j}}{\varphi_{j}^{2}}\right]\Bigg{)}$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\varphi}^{2}\right)\right]\propto\sum_{j=1}^{p}\Bigg{(}\log
d_{j}^{\star}-2\mathbb{E}_{\varphi^{2}}\left[\log\varphi_{j}^{2}\right]-\frac{1}{2}\mathbb{E}_{\varphi^{2}}\left[\varphi_{j}^{-2}\right]-\frac{1}{2s_{0}^{2}}\mathbb{E}_{\lambda\cdot\gamma\cdot\varphi^{2}}\left[\frac{\left(\lambda_{j}-\mathbf{z}_{j}^{\intercal}\boldsymbol{\gamma}\right)^{2}}{\varphi_{j}^{2}}\right]\Bigg{)}$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\boldsymbol{\gamma}\right)\right]\propto-\frac{1}{2}\log|\boldsymbol{\Sigma}_{\gamma}^{\star}|-\frac{M}{2s_{0}^{2}}$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\kappa^{2}\right)\right]\propto\sum_{d=1}^{D}\Bigg{(}\left(a_{d}+\frac{m_{d}}{2}\right)\log
f_{d}^{\star}-(a_{d}+\frac{m_{d}}{2}+1)\mathbb{E}_{\kappa^{2}}\left[\log\kappa_{d}^{2}\right]-b_{d}\mathbb{E}_{\kappa^{2}}\left[\kappa_{d}^{-2}\right]-$
$\displaystyle\hskip
96.73918pt\frac{1}{2}\mathbb{E}_{\gamma\cdot\kappa^{2}}\left[\frac{\boldsymbol{\gamma}_{d}^{\intercal}\boldsymbol{\gamma}_{d}}{\kappa_{d}^{2}}\right]\Bigg{)}$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\tau^{2}\right)\right]\propto\left(\frac{p}{2}+1\right)\log\left(g^{\star}\right)-\left(\frac{p}{2}+2\right)\mathbb{E}_{\tau^{2}}\left[\log\tau^{2}\right]-\mathbb{E}_{\tau^{2}}\left[\frac{1}{\zeta\tau^{2}}\right]-\frac{1}{2}\sum_{j=0}^{p}\mathbb{E}_{\beta_{0}\cdot\beta\cdot\lambda_{0}^{2}\cdot\lambda\cdot\tau^{2}}\left[\frac{\beta_{j}^{2}}{\tau^{2}\lambda_{j}^{2}}\right]$
$\displaystyle\mathbb{E}_{q^{\star}}\left[\log
q^{\star}\left(\zeta\right)\right]\propto\log
h^{\star}-2\mathbb{E}_{\zeta}\left[\log\zeta\right]-\mathbb{E}_{\zeta}\left[\zeta^{-1}\right]-\mathbb{E}_{\zeta}\left[\frac{1}{\zeta\tau^{2}}\right],$

where
$\Phi_{i}=\Phi\left(-\mathbf{x}_{i}^{\intercal}\mathbb{E}_{\beta}\left[\boldsymbol{\beta}\right]\right)$.

##  Appendix E Web Appendix E

###  E.1 Computational methods for the variational algorithm

The Variational algorithm has the advantage of not requiring the explicit
inversion of covariance matrix
$\boldsymbol{\Sigma}_{\beta}^{\star}=\left(\mathbf{X}^{\intercal}\mathbf{X}+\boldsymbol{\Delta}\right)^{-1}$,
where
$\boldsymbol{\Delta}=\mathbb{E}_{\lambda\cdot\tau^{2}}\left[\tau^{-2}\Lambda^{-2}\right]$
is diagonal, which becomes computationally infeasible when the number of
covariates $p$ increases. The quantities needed are the diagonal entries of
$\boldsymbol{\Sigma}_{\beta}^{\star}$, the posterior mean
$\boldsymbol{\mu}_{\beta}^{\star}=\boldsymbol{\Sigma}_{\beta}^{\star}\left(\mathbf{X}^{\intercal}\mathbf{y}\right)$
and
$tr\left(\mathbf{X}^{\intercal}\mathbf{X}\boldsymbol{\Sigma}_{\beta}^{\star}\right)$.
By means of the Woodbury identity, the needed quantities can be efficiently
computed as

$\displaystyle\text{diag}\left(\boldsymbol{\Sigma}_{\beta}^{\star}\right)=$
$\displaystyle\;\text{diag}\left(\boldsymbol{\Delta}^{-1}\right)-\text{diag}\left(\boldsymbol{\Delta}^{-1}\mathbf{X}^{\intercal}\left(\mathbf{I}_{n}+\mathbf{X}\boldsymbol{\Delta}^{-1}\mathbf{X}^{\intercal}\right)^{-1}\mathbf{X}\boldsymbol{\Delta}^{-1}\right)$
$\displaystyle=$
$\displaystyle\;\text{diag}\left(\boldsymbol{\Delta}^{-1}\right)-\left[\left(\boldsymbol{\Delta}^{-1}\mathbf{X}^{\intercal}\right)\left(\mathbf{I}_{n}+\mathbf{X}\boldsymbol{\Delta}^{-1}\mathbf{X}^{\intercal}\right)^{-1}\circ\left(\boldsymbol{\Delta}^{-1}\mathbf{X}^{\intercal}\right)\right]\cdot\mathbf{1}_{n\times
1},$ (8) $\displaystyle\boldsymbol{\mu}_{\beta}^{\star}=$
$\displaystyle\;\left(\boldsymbol{\Sigma}_{\beta}^{\star}\mathbf{X}^{\intercal}\right)\mathbf{y}$
$\displaystyle=$
$\displaystyle\;\left[\left(\boldsymbol{\Delta}^{-1}\mathbf{X}^{\intercal}\right)-\left(\boldsymbol{\Delta}^{-1}\mathbf{X}^{\intercal}\right)\left(\mathbf{I}_{n}+\mathbf{X}\boldsymbol{\Delta}^{-1}\mathbf{X}^{\intercal}\right)^{-1}\left(\mathbf{X}\boldsymbol{\Delta}^{-1}\mathbf{X}^{\intercal}\right)\right]\mathbf{y},$
(9) $\displaystyle
tr\left(\mathbf{X}\boldsymbol{\Sigma}_{\beta}^{\star}\mathbf{X}^{\intercal}\right)=$
$\displaystyle\;\sum_{i=1}^{n}\mathbf{x}_{i}^{\intercal}\boldsymbol{\Sigma}_{\beta}^{\star}\mathbf{x}_{i}=\sum_{i=1}^{n}\mathbf{x}_{i}^{\intercal}\left(\boldsymbol{\Sigma}_{\beta}^{\star}\mathbf{x}_{i}\right),$
(10)

where $\circ$ denotes the Hadamard product, $\mathbf{I}_{n}$ is the $n\times
n$ identity matrix and $\mathbf{1}_{n\times 1}$ denotes the $n$-dimensional
vector of ones. In (10) the quantity
$\boldsymbol{\Sigma}_{\beta}^{\star}\mathbf{x}_{i}$ has already been computed
when evaluating vector $\boldsymbol{\mu}_{\beta}^{\star}$ and all the matrix
products in (8)-(10) can be evaluated with $\mathcal{O}\left(n^{2}p\right)$
operations, which is linear in $p$. Note that in the probit regression the
last quantity (10) is not needed, thus the posterior mean can be efficiently
computed as

$\boldsymbol{\mu}_{\beta}^{\star}=\boldsymbol{\Delta}^{-1}\left(\mathbf{X}^{\intercal}\mathbb{E}_{w}\left[\mathbf{w}\right]\right)-\left(\boldsymbol{\Delta}^{-1}\mathbf{X}^{\intercal}\right)\left(\mathbf{I}_{n}+\mathbf{X}\boldsymbol{\Delta}^{-1}\mathbf{X}^{\intercal}\right)^{-1}\left(\mathbf{X}\boldsymbol{\Delta}^{-1}\right)\left(\mathbf{X}^{\intercal}\mathbb{E}_{w}\left[\mathbf{w}\right]\right),$

which only involves matrix-vector products, further improving the
computational efficiency of the algorithm. Finally, the determinant can be
efficiently evaluated as

$\left|\boldsymbol{\Sigma}_{\beta}^{\star}\right|=\left|\boldsymbol{\Delta}\right|/\left|\mathbf{I}_{n}+\mathbf{X}\boldsymbol{\Delta}^{-1}\mathbf{X}^{\intercal}\right|.$

##  Appendix F Web Appendix F

###  F.1 Details on the simulation studies

For all the considered cases we rely on the following simulation scheme. Let
$\boldsymbol{\beta}^{0}$ be the true $(p+1)$-dimensional regression parameter
vector. We set the number of true non-zero coefficients to $p_{0}$ (intercept
excluded). The entries of design matrices are sampled independently as
$x_{ij}\sim\mathcal{N}(0,1)$, $i=1,\dots,n$, $j=1,\dots,p$, whereas the
response variable $\mathbf{y}$ and the $(p+1)$-dimensional true regression
coefficient vector are sampled following a modified version of the scheme in
Johnson, (2013). Specifically

$\displaystyle y_{i}=$
$\displaystyle\;\mathbf{x}_{i}^{\intercal}\boldsymbol{\beta}^{0}+\varepsilon_{i},\quad\varepsilon_{i}\sim\mathcal{N}(0,1),\quad
i=1,\dots,n,$ $\displaystyle\beta_{j}^{0}=$
$\displaystyle\;\begin{cases}v_{0}|t|\quad\text{if}\mbox{ }j=0,\\\
\left(-1\right)^{u}\left(v^{2}\log(n)/\sqrt{n}+v|t|\right)\quad\text{if}\mbox{
}j=1,\dots,p_{0},\\\ 0\quad\text{otherwise},\end{cases}$

where $v_{0}^{2}=0.5$, $v^{2}=0.75$, $u\sim\mathcal{B}(0.4)$ and
$t\sim\mathcal{N}\left(0,1\right)$.

###  F.2 Co-data simulation

Here we give insights on how co-data information is simulated for the
assessment of the Variational approximation accuracy with respect to the Gibbs
sampler in Section 1 of the main article. Moreover we provide the graphical
representation of the AUC for variable selection for both methods.

The four degrees of co-data information are simulated as:

  * $G0\big{)}$ no co-data set-up: we include in the co-data matrix only the intercept, therefore $\mathbf{Z}=\mathbf{1}_{p}$;

  * $G1\big{)}$ non-informative set-up: a binary co-data source is included in the model by randomly selecting $30$ regressors, therefore the co-data matrix $\mathbf{Z}$ is created from the binary vector $\mathbf{z}\in\left\\{0,1\right\\}^{p}$, where $z_{j}=1$ if the $j$-th variable is selected, $z_{j}=0$ otherwise;

  * $G2\big{)}$ informative set-up: a binary co-data source is included in the model by randomly selecting $20$ of the true non-zero regressors and $10$ of the true zero regressors, therefore the co-data matrix $\mathbf{Z}$ is created from the binary vector $\mathbf{z}\in\left\\{0,1\right\\}^{p}$, where $z_{j}=1$ if the $j$-th variable is selected, $z_{j}=0$ otherwise;

  * $G3\big{)}$ perfect co-data information set-up: the co-data matrix $\mathbf{Z}$ is created from the binary vector $\mathbf{z}\in\left\\{0,1\right\\}^{p}$, where $z_{j}=1$ if $\beta_{j}^{0}\neq 0$, $z_{j}=0$ otherwise.

##  Appendix G Web Appendix G

###  G.1 Web Algorithm 1

1 Input: $B,bn\in\mathbb{N}$, $\mathbf{y}\in\mathbb{R}^{n}$,
$\mathbf{X}\in\mathbb{R}^{n\times(p+1)}$,
$\mathbf{Z}_{1},\dots,\mathbf{Z}_{D}\in\mathbb{R}^{p\times m_{d}}$,
$\mathbf{a},\mathbf{b}\in\mathbb{R}_{+}^{D}$, $s_{0}^{2}\in\mathbb{R}_{+}$;

2

// Set $\lambda_{0}=1$, $\boldsymbol{\lambda}=\mathbf{1}_{p}$,
$\boldsymbol{\gamma}=\mathbf{0}_{M}$ and $\sigma^{2},\tau^{2}=1$ and sample
parameters $\boldsymbol{\varphi}$, $\boldsymbol{\kappa}^{2}$ and $\zeta$ from
their prior distributions

3

4for  _$b=1\text{ to }B$_ do

5 sample
$\boldsymbol{\beta}\mid\mathbf{y},\mathbf{X},\sigma^{2},\tau^{2},\lambda_{0}^{2},\boldsymbol{\lambda}\sim\mathcal{N}_{p+1}\left(\boldsymbol{\Sigma}_{\beta}^{\star}\mathbf{X}^{\intercal}\mathbf{y},\sigma^{2}\boldsymbol{\Sigma}_{\beta}^{\star}\right)$,
where
$\boldsymbol{\Sigma}_{\beta}^{\star}=\left(\mathbf{X}^{\intercal}\mathbf{X}+\tau^{-2}\boldsymbol{\Lambda}^{-2}\right)^{-1}$;

6

7 sample
$\lambda_{0}^{2}\mid\beta_{0},\tau^{2},\sigma^{2}\sim\mathcal{IG}\left(1,\frac{1}{\psi_{0}}+\frac{\beta_{0}^{2}}{2\sigma^{2}\tau^{2}}\right)$;

8

9 sample
$\psi_{0}\mid\lambda_{0}^{2}\sim\mathcal{IG}\left(1,1+\frac{1}{\lambda_{0}^{2}}\right)$;

10 for  _$j=1\text{ to }p$_ do

11 sample
$\lambda_{j}\mid\mathbf{Z},\beta_{j},\sigma^{2},\tau^{2},\boldsymbol{\gamma},\varphi_{j}^{2}\propto\lambda_{j}^{-1}e^{-\frac{\beta_{j}^{2}}{2\sigma^{2}\tau^{2}\lambda_{j}^{2}}-\frac{\lambda_{j}^{2}}{2s_{0}^{2}\varphi_{j}^{2}}+\frac{\mu_{j}\lambda_{j}}{s_{0}^{2}\varphi_{j}^{2}}}\cdot\mathbb{I}_{\left(\lambda_{j}>0\right)}$
following the procedure in Section B.1 of Supporting Information;

12 sample
$\varphi_{j}^{2}\mid\mathbf{Z},\lambda_{j},\boldsymbol{\gamma}\sim\mathcal{IG}\left(1,\frac{1}{2}+\frac{(\lambda_{j}-\mu_{j})^{2}}{2s_{0}^{2}}\right)$;

13

14 end for

15 sample
$\boldsymbol{\gamma}\mid\mathbf{Z},\boldsymbol{\lambda},\boldsymbol{\varphi},\boldsymbol{\kappa}^{2}\sim\mathcal{N}_{M}\left(\boldsymbol{\Sigma}_{\gamma}^{\star}\left(\mathbf{Z}^{\intercal}\boldsymbol{\Phi}^{-2}\boldsymbol{\lambda}\right),s_{0}^{2}\boldsymbol{\Sigma}_{\gamma}^{\star}\right)$,
where
$\boldsymbol{\Sigma}_{\gamma}^{\star}=\left(\mathbf{Z}^{\intercal}\boldsymbol{\Phi}^{-2}\mathbf{Z}+s_{0}^{2}\mathbf{D}_{\kappa}^{-1}\right)^{-1}$
and
$\mathbf{D}_{\kappa}=\text{diag}\left(\kappa_{1}^{2}\mathbf{1}_{m_{1}},\dots,\kappa_{D}^{2}\mathbf{1}_{m_{D}}\right)$;

16 for  _$d=1\text{ to }D$_ do

17 sample
$\kappa_{d}^{2}\mid\boldsymbol{\gamma}_{d}\sim\mathcal{IG}\left(a_{d}+\frac{m_{d}}{2},b_{d}+\frac{\boldsymbol{\gamma}_{d}^{\intercal}\boldsymbol{\gamma}_{d}}{2}\right)$;

18

19 end for

20 sample
$\tau^{2}\mid\boldsymbol{\beta},\sigma^{2},\lambda_{0}^{2},\boldsymbol{\lambda},\zeta\sim\mathcal{IG}\left(\frac{p}{2}+1,\frac{1}{\zeta}+\frac{\beta_{0}^{2}}{2\sigma^{2}\lambda_{0}^{2}}+\frac{1}{2\sigma^{2}}\sum_{j=1}^{p}\frac{\beta_{j}^{2}}{\lambda_{j}^{2}}\right)$;

21 sample
$\zeta\mid\tau^{2}\sim\mathcal{IG}\left(1,1+\frac{1}{\tau^{2}}\right)$;

22 sample
$\sigma^{2}\mid\mathbf{y},\mathbf{X},\boldsymbol{\beta},\tau^{2},\lambda_{0}^{2},\boldsymbol{\lambda}\sim\mathcal{IG}\left(v+\frac{n+p+1}{2},q+\frac{1}{2}\|\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\|_{2}^{2}+\frac{\beta_{0}^{2}}{2\tau^{2}\lambda_{0}^{2}}+\frac{1}{2\tau^{2}}\sum_{j=1}^{p}\frac{\beta_{j}^{2}}{\lambda_{j}^{2}}\right)$;

23 if  _$b >bn$_ then

24 save $\boldsymbol{\beta}$, $\sigma^{2}$, $\tau^{2}$, $\lambda_{0}^{2}$,
$\boldsymbol{\lambda}$ , $\boldsymbol{\gamma}$ and $\boldsymbol{\kappa}^{2}$

25 end if

26

27 end for

28return saved values;

Algorithm 1 Gibbs sampler for Informative Horseshoe regression

###  G.2 Web Algorithm 2

1 Input: $\mathbf{y}\in\mathbb{R}^{n}$,
$\mathbf{X}\in\mathbb{R}^{n\times(p+1)}$,
$\mathbf{Z}_{1},\dots,\mathbf{Z}_{D}\in\mathbb{R}^{p\times m_{d}}$,
$v,q\in\mathbb{R}^{+}$, $\mathbf{a},\mathbf{b}\in\mathbb{R}_{+}^{D}$,
$s_{0}^{2}\in\mathbb{R}_{+}$;

2

// Set $b=1$, $\epsilon=10^{-3}$, $\mathcal{L}^{(0)}=-\infty$ and initialize
all the needed moments

3

4while  _$\mathcal{L}^{(b)}-\mathcal{L}^{(b-1)} >\epsilon$_ do

5 Update parameter $\boldsymbol{\mu}_{\beta}^{\star}$ and compute the
quantities $\text{diag}\left(\boldsymbol{\Sigma}_{\beta}^{\star}\right)$,
$|\boldsymbol{\Sigma}_{\beta}^{\star}|$ and
$\mathbf{X}\boldsymbol{\Sigma}_{\beta}^{\star}\mathbf{X}^{\intercal}$ as in
Section E.1 of Supporting Information;

6 Update parameters $a_{0}^{\star}$ and $k_{0}^{\star}$;

7 Update parameters $a_{j}^{\star}$, $b_{j}^{\star}$, $c_{j}^{\star}$ and
$d_{j}^{\star}$ and evaluate the normalizing constant $s_{j}$,
$\mathbb{E}_{\lambda}\left[\lambda_{j}\right]$,
$\mathbb{E}_{\lambda}\left[\lambda_{j}^{2}\right]$ and
$\mathbb{E}_{\lambda}\left[\lambda_{j}^{-2}\right]$ with numerical
integration, for $j=1,\dots,p$;

8 Update parameters $\boldsymbol{\mu}_{\gamma}^{\star}$ and
$\boldsymbol{\Sigma}_{\gamma}^{\star}$;

9 Update parameters $e_{d}^{\star}$ and $f_{d}^{\star}$, for $d=1,\dots,D$;

10 Update parameters $g^{\star}$ and $h^{\star}$;

11 Update parameter $l^{\star}$;

12 Compute $\mathcal{L}^{(b)}$ and set $b=b+1$;

13

14 end while

15return $\boldsymbol{\mu}_{\beta}^{\star}$,
$\boldsymbol{\Sigma}_{\beta}^{\star}$, $a_{0}^{\star}$, $\mathbf{a}^{\star}$,
$\mathbf{b}^{\star}$, $\mathbf{c}^{\star}$,
$\boldsymbol{\mu}_{\gamma}^{\star}$, $\boldsymbol{\Sigma}_{\gamma}^{\star}$,
$\mathbf{e}^{\star}$, $\mathbf{f}^{\star}$, $g^{\star}$ and $l^{\star}$;

Algorithm 2 Variational Bayes approximation for informative Horseshoe
regression

###  G.3 Web Algorithm 3

1 Input: $\mathbf{y}\in\left\\{0,1\right\\}^{n}$,
$\mathbf{X}\in\mathbb{R}^{n\times(p+1)}$,
$\mathbf{Z}_{1},\dots,\mathbf{Z}_{D}\in\mathbb{R}^{p\times m_{d}}$,
$\mathbf{a},\mathbf{b}\in\mathbb{R}_{+}^{D}$, $s_{0}^{2}\in\mathbb{R}_{+}$;

2

// Set $b=1$, $\epsilon=10^{-3}$, $\mathcal{L}^{(0)}=-\infty$ and initialize
all the needed moments

3

4while  _$\mathcal{L}^{(b)}-\mathcal{L}^{(b-1)} >\epsilon$_ do

5 Update parameter $\mu_{i}^{\star}$, for $i=1,\dots,n$;

6 Update parameter $\boldsymbol{\mu}_{\beta}^{\star}$ and compute the
quantities $\text{diag}\left(\boldsymbol{\Sigma}_{\beta}^{\star}\right)$ and
$|\boldsymbol{\Sigma}_{\beta}^{\star}|$ as in Section E.1 of Supporting
Information;

7 Update parameters $a_{0}^{\star}$ and $k_{0}^{\star}$;

8 Update parameters $a_{j}^{\star}$, $b_{j}^{\star}$, $c_{j}^{\star}$ and
$d_{j}^{\star}$ and evaluate the normalizing constant $s_{j}$,
$\mathbb{E}_{\lambda}\left[\lambda_{j}\right]$,
$\mathbb{E}_{\lambda}\left[\lambda_{j}^{2}\right]$ and
$\mathbb{E}_{\lambda}\left[\lambda_{j}^{-2}\right]$ with numerical
integration, for $j=1,\dots,p$;

9 Update parameters $\boldsymbol{\mu}_{\gamma}^{\star}$ and
$\boldsymbol{\Sigma}_{\gamma}^{\star}$;

10 Update parameters $e_{d}^{\star}$ and $f_{d}^{\star}$, for $d=1,\dots,D$;

11 Update parameters $g^{\star}$ and $h^{\star}$;

12 Compute $\mathcal{L}^{(b)}$ and set $b=b+1$;

13

14 end while

15return $\boldsymbol{\mu}_{\beta}^{\star}$, $a_{0}^{\star}$,
$\mathbf{a}^{\star}$, $\mathbf{b}^{\star}$, $\mathbf{c}^{\star}$,
$\boldsymbol{\mu}_{\gamma}^{\star}$, $\boldsymbol{\Sigma}_{\gamma}^{\star}$,
$\mathbf{e}^{\star}$, $\mathbf{f}^{\star}$, $g^{\star}$;

Algorithm 3 Variational Bayes approximation for probit informative Horseshoe
regression

##  Appendix H Web Appendix H

###  H.1 Web Figure 1

Figure 5: Variable selection with Gibbs sampler and Variational algorithm; the
$AUC$ is averaged over $10$ replicates of each experiments.

###  H.2 Web Figure 2

Figure 6: Estimated standard deviations of regression coefficients
$\boldsymbol{\beta}$ with Gibbs sampler and Variational algorithm; results are
evaluated over $10$ replicates of each experiments. Data have been aggregated
for the different type of co-data.