elements of $\boldsymbol{\beta}^{\text{S}}$ have been specified, with
$s_{\tau^{2}}=10^{5}$ to limit prior information about the degree of sparsity.
Along all the simulation study, the MFVB approximations were obtained by
running 200 iterations of the streamlined or naïve algorithms, i.e. four times
the number of iterations used in the numerical experiments of Nolan et al.
(2020). This number of MFVB iterations was chosen to guarantee appropriate
convergence of the variational algorithms, as well as a fair computational
time comparison between naïve and streamlined MFVB. Maestrini and Wand (2021)
suggest running variational algorithms until the relative change in the
variational parameters reaches a certain threshold. In our numerical
experiments, the 200 MFVB iterations guaranteed that the relative change in
the variational parameter estimates was smaller that $10^{-3}$ and for this or
smaller values the accuracy and fixed effect selection performance of MFVB was
not impacted. Variational algorithm convergence could also be assessed by
monitoring the lower bound at every iteration or after a certain number of
iterations; however, deriving the variational lower bound is subject to
calculation errors and its computation can be time consuming.

###  6.1 Accuracy Assessment

We measured the quality of the variational approximations using the accuracy
index proposed in Section 3.4 of Faes _et al._ (2011). For a generic
univariate parameter $\theta\in\Theta$, this is defined as

$\mbox{Accuracy}(\theta)\equiv\left(1-\dfrac{1}{2}\int_{\Theta}\left|\mathfrak{q}^{*}(\theta)-\mathfrak{p}(\theta|\boldsymbol{y})\right|\mathrm{d}\theta\right)\%,$
(27)

where $\mathfrak{p}(\theta|\boldsymbol{y})$ is the marginal posterior density
of $\theta$ and $\mathfrak{q}^{*}(\theta)$ is the associated optimal
approximating density function obtained via MFVB. This index takes values
between 0% and 100%, with a score of 100% indicating perfect matching between
$\mathfrak{q}^{*}(\theta)$ and $\mathfrak{p}(\theta|\boldsymbol{y})$, and 0%
if the densities have no overlapping mass. In practice, computation of the
true $\mathfrak{p}(\theta|\boldsymbol{y})$ is numerically challenging, so we
employed binned kernel density estimation with direct plug-in bandwidth
selection (e.g. Section 3.6.1 of Wand & Jones, 1995) and applied it to the
MCMC samples using the R package KernSmooth (Wand, 2020). The MCMC samplings
were performed using a warmup of length 5,000 followed by 100,000 iterations,
to which we applied a thinning value of 20. The integration in (27) was then
accurately performed via trapezoidal numerical quadrature.

Figure 2: Side-by-side boxplots of the accuracy scores from the simulation
study for a selection of model parameters and random effects. Outliers are
displayed as solid points.

Figure 2 shows the accuracy results through boxplots for the 50 data
replicates. The boxplots refer to all the entries of
$\boldsymbol{\beta}^{\text{R}}$ and $\boldsymbol{\beta}^{\text{A}}$, to 10
elements of $\boldsymbol{\beta}^{\text{S}}$ chosen such that 5 of them have
non-zero values ($\beta^{\text{S}}_{1}$, $\beta^{\text{S}}_{3}$,
$\beta^{\text{S}}_{4}$, $\beta^{\text{S}}_{6}$ and $\beta^{\text{S}}_{8}$) and
5 are null ($\beta^{\text{S}}_{11}$, $\beta^{\text{S}}_{23}$,
$\beta^{\text{S}}_{24}$, $\beta^{\text{S}}_{36}$ and $\beta^{\text{S}}_{48}$),
the intercept- and slope-associated elements of
$\boldsymbol{u}^{\mbox{\tiny\rm L1}}_{i}$ and $\boldsymbol{u}^{\mbox{\tiny\rm
L2}}_{ij}$, for $1\leq i\leq 3$ and $j=1$, the entries of
$\boldsymbol{\Sigma}^{\mbox{\tiny\rm L1}}$ and
$\boldsymbol{\Sigma}^{\mbox{\tiny\rm L2}}$, $\sigma^{2}$, and $\tau^{2}$. The
intercept- and slope-associated parameters are identified by the subscripts
“int” and “slp”, respectively.

Variational approximations showed high accuracy scores for all the model
parameters considered and across the different data replications. All the
fixed effects subject to selection and having non-zero true values exhibited
accuracy scores greater than 90%. In contrast, those having true values equal
to zero had lower accuracy scores between 75% and 85% due to the spiky
marginal posterior densities being approximated by Gaussian variational
densities. All the other fixed effects parameters, random effects and variance
parameters showed accuracy scores greater than 90%. The accuracy scores of the
global variance parameter $\tau^{2}$ were under 50%. Despite not directly
shown, all the $\mathfrak{q}^{*}(\zeta_{h})$’s had accuracy scores between 75%
and 80%. Similar results were obtained for the alternative global-local prior
specifications.

Figure 3: Approximate posterior density functions of some of the three-lever
random effects model parameters obtained from the first replication of the
simulation study. Each plot shows the optimal approximate posterior density
function $\mathfrak{q}^{*}(\theta)$ obtained via MFVB (black dashed curves)
and the MCMC-based $\mathfrak{p}(\theta|\boldsymbol{y})$ densities (grey
curves). A vertical line indicates the parameter true value. The percentages
of accuracy are also provided.

Figure 3 displays the MFVB and MCMC approximate posterior densities, together
with the associated accuracy scores, obtained for the first data replication
of the simulation study. These plots allow to visually assess the quality of
the approximation and show that the approximate posterior density functions
are generally concentrated around the true parameter values. Notice that the
global-local prior shrinks the MCMC marginal posterior densities of the
$\beta^{\text{S}}$-parameters having null true value towards zero. The
corresponding $\mathfrak{q}^{*}$-densities are Gaussian and, although they are
not able to capture peak and tail behaviors, they provide satisfactory
approximations. Notice also that the $\tau^{2}$ approximate posterior
densities are concentrated around a very small value, as expected for the
sparse design setting we considered in the simulation study.

###  6.2 Fixed Effects Selection Assessment

We assessed fixed effects selection performances by running Algorithm 3 on the
same 50 data replications, also experimenting the other two global-local
priors considered in this work: the Laplace and Normal-Exponential-Gamma with
$\lambda=0.25$. We also considered the Gaussian prior specification described
in Nolan _et al._ (2020) with hyperparameters
$\boldsymbol{\mu}_{\boldsymbol{\beta}^{\text{S}}}=\boldsymbol{0}$ and
$\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{\text{S}}}=10^{10}\boldsymbol{I}$.
Then, for each data replication and each of the different four prior
specifications considered we used the approximating densities corresponding to
fixed effects subject to selection to perform the SAVS procedure of Algorithm
1. We employed this procedure also for the approximate posteriors obtained
from the Gaussian prior to provide a comparison with a prior that does not
belong to the global-local family.

Let TP (true positives) denote the number of selected fixed effects having
true value different from zero and TN (true negatives) denote the number of
unselected irrelevant fixed effects having true value equal to zero. We
measured the fixed effects selection performance for each prior choice using
the $F_{1}$-score (Van Rijsbergen, 1979), which is defined as:

$F_{1}\equiv\left(\frac{2\times\text{precision}\times\text{recall}}{\text{precision}+\text{recall}}\right)\times
100\%,$ (28)

with $\text{precision}\equiv\text{TP}/(\text{TP}+\text{FP})$ and
$\text{recall}\equiv\text{TP}/(\text{TP}+\text{FN})$, where FP and FN denote
the number of false positives (incorrectly selected fixed effects) and false
negatives (relevant fixed effects that have not been selected), respectively.
This index takes values between 0% and 100%, with higher values to be
preferred.

The median $F_{1}$-score was 63.5% (1st quartile: 43.5%; 3rd quartile: 94.2%)
for the Gaussian prior and 95.24% (1st quartile: 80%; 3rd quartile: 100%) for
the Laplace prior. Both the Horseshoe and Negative-Exponential-Gamma priors
exhibited a $F_{1}$-score equal to 100% for all the 50 data replications of
the simulation study, meaning that perfect selection ($\text{TP}=10$ and
$\text{TN}=40$) was always achieved. From this simplified simulated sparse
data scenario it is apparent that global-local priors effectively provide
better variable selection performances than the Gaussian prior. The lower
performances of the Gaussian and Laplace priors are due to the SAVS procedure
being applied to optimal approximate $\mathfrak{q}^{*}$-densities that have
not been properly shrunk towards zero and so irrelevant fixed effects tend to
be selected. Nonetheless, all the four different priors gave $\text{FP}=0$,
meaning that they did not incorrectly select irrelevant fixed effects.

###  6.3 Speed and Memory Saving Assessment

Streamlined variational inference has been conceived to obtain efficient
implementations of variational algorithms. Hence, the assessment of speed and
memory savings is another important aspect. We assessed speed and memory
performances of the streamlined variational Algorithm 3 and its naïve
counterpart, whose updates are described in Section 5.1. Four different group
numbers and three different lengths for the vector of fixed effects to be
selected were considered, namely $m\in\\{10,50,100,200\\}$ and
$p_{\text{S}}\in\\{25,100,200\\}$, aiming to explore scalability of the
streamlined methodology to high dimensions. For each combination of $m$ and
$p_{\text{S}}$, we simulated 10 data replications from model (22) with a
random intercept and one slope for a single continuous predictor, choosing the
sub-group dimensions $n_{i}$ uniformly on the discrete set
$\\{10,\ldots,20\\}$ and the sub-group specific unit dimensions $o_{ij}$
uniformly on the discrete set $\\{20,\ldots,30\\}$. This setting allows to
test models with heterogeneous dimensions but the same number of groups $m$.
We considered a Horseshoe prior for the fixed effects subject to selection,
and all the other model dimensions and hyperparameters specifications were the
same as those described before.

For each simulated dataset, we collected the computational timings and the
total size of the input data required for performing both the streamlined and
naïve algorithms. We also ran MCMC (warmup of length 5,000 followed by 25,000
iterations, to which we applied a thinning value of 5) and recorded its
computational timings, although MFVB and MCMC do not admit a genuine
comparison due to the fact that both depend upon different convergence
requirements, as explained in Section 5.1 of Ormerod _et al._ (2017).
Moreover, our efficient implementation of MCMC is based on independent
resampling from the full conditional densities of $\boldsymbol{\beta}$,
$\boldsymbol{u}^{\mbox{\tiny\rm L1}}_{i}$ and $\boldsymbol{u}^{\mbox{\tiny\rm
L2}}_{ij}$, whereas MFVB uses a joint approximation
$\mathfrak{q}(\boldsymbol{\beta},\boldsymbol{u})$ for all these vectors.
Nonetheless, MCMC timings are also reported to provide an intuition on the
computational effort required for sampling from the posterior distribution
when $m$ and $p_{\text{S}}$ increase.

|  | Total runtime of the algorithm (seconds) | Total size of the required data inputs (megabytes)
---|---|---|---
$m$ | $p_{\text{S}}$ | _Streamlined MFVB_ | _Naïve MFVB_ | $\frac{\text{\it Na{\"{i}}ve MFVB}}{\text{\it Streamlined MFVB}}$ | _MCMC_ | _Streamlined MFVB_ | _Naïve MFVB_ | $\frac{\text{\it Na{\"{i}}ve MFVB}}{\text{\it Streamlined MFVB}}$
10 | 25 | $0.70_{(0.06)}$ | $3.86_{(0.50)}$ | $5.54$ | $13.03_{(0.96)}$ | $1.39_{(0.07)}$ | $11.06_{(1.00)}$ | $7.97$
| 100 | $3.90_{(0.23)}$ | $6.51_{(0.82)}$ | $1.67$ | $46.99_{(1.47)}$ | $3.68_{(0.24)}$ | $12.73_{(1.34)}$ | $3.46$
| 200 | $9.71_{(1.22)}$ | $12.09_{(2.04)}$ | $1.24$ | $176.22_{(9.37)}$ | $6.60_{(0.52)}$ | $15.09_{(1.89)}$ | $2.29$
50 | 25 | $3.29_{(0.11)}$ | $365.07_{(37.10)}$ | $111.01$ | $63.28_{(3.76)}$ | $6.69_{(0.20)}$ | $240.40_{(13.92)}$ | $35.91$
| 100 | $19.66_{(0.87)}$ | $415.29_{(51.66)}$ | $21.12$ | $138.42_{(5.72)}$ | $18.37_{(0.81)}$ | $254.04_{(20.84)}$ | $13.83$
| 200 | $49.25_{(1.47)}$ | $501.33_{(39.76)}$ | $10.18$ | $305.25_{(6.51)}$ | $34.00_{(1.04)}$ | $269.34_{(14.90)}$ | $7.92$
100 | 25 | $6.77_{(0.32)}$ | $2877.33_{(177.77)}$ | $424.72$ | $151.05_{(10.48)}$ | $13.56_{(0.34)}$ | $967.34_{(42.15)}$ | $71.33$
| 100 | $40.42_{(1.66)}$ | $3172.66_{(116.83)}$ | $78.50$ | $266.71_{(10.06)}$ | $36.77_{(1.29)}$ | $1010.25_{(26.94)}$ | $27.47$
| 200 | $97.81_{(2.07)}$ | $3403.55_{(204.38)}$ | $34.80$ | $494.15_{(9.34)}$ | $67.53_{(1.31)}$ | $1013.90_{(48.67)}$ | $15.01$
200 | 25 | $13.30_{(0.30)}$ |  $>$ 5 hours | $>1355$ | $328.87_{(20.04)}$ | $26.90_{(0.56)}$ | $3817.06_{(113.85)}$ | $141.88$
| 100 | $79.97_{(1.18)}$ |  $>$ 5 hours | $>225$ | $578.88_{(13.45)}$ | $73.74_{(1.01)}$ | $3941.15_{(97.90)}$ | $53.45$
| 200 | $197.20_{(4.13)}$ |  $>$ 5 hours | $>95$ | $904.64_{(15.77)}$ | $135.84_{(2.73)}$ | $3991.66_{(102.48)}$ | $29.38$

Table 2: Average (standard deviation of) elapsed computing times in seconds
and average (standard deviation of) total size of required data inputs in
megabytes for fitting model (22) with three-level random effects specification
and $p_{\text{S}}$ fixed effects subject to penalization via Horseshoe prior.
Results are shown for different group sizes $m$ and different values for
$p_{\text{S}}$.

The tabulated results are shown in Table 2. The “ratio” columns help
understand the gain obtained employing the streamlined MFVB methodology over
its naïve counterpart. For increasing $m$, streamlined MFVB reached
convergence faster than naïve MFVB. Notice that in the biggest scenario under
examination (last row of the table), the streamlined MFVB algorithm ran in
less then 4 minutes on average, while naïve MFVB required more than 5 hours.
Bigger scenarios are computationally demanding for the naïve implementation
and may fail to run due to excessive storage demand, whilst we did not
experience these issues with Algorithm 3. Similar comments apply to the huge
saving of memory allocation for the required input data provided by the
streamlined MFVB implementation.

We conclude the discussion by noticing that our methodology takes the
approximate joint posterior dependence between the fixed effects and random
effects parameters into account through the multivariate Gaussian
approximating density $\mathfrak{q}(\boldsymbol{\beta},\boldsymbol{u})$. This
choice allows to better capture the _a posteriori_ covariance structure
between $\boldsymbol{\beta}$ and $\boldsymbol{u}$ and ensures better results
in terms of approximation accuracy. Nevertheless, alternative and less
restrictive factorizations can be considered, especially if the quality of the
approximation can be sacrificed and faster algorithms are desired. For
instance, the additional factorization
$\mathfrak{q}(\boldsymbol{\beta},\boldsymbol{u})=\mathfrak{q}(\boldsymbol{\beta}^{\text{R}},\boldsymbol{u})\mathfrak{q}(\boldsymbol{\beta}^{\text{A}},\boldsymbol{\beta}^{\text{S}})$
could remarkably speed up the computations if $p_{\text{R}}$ is small,
regardless of the size of $p_{\text{A}}$ and $p_{\text{S}}$. For implementing
this and any other additional independence constraints, mean field
restrictions such as (24) and the Algorithms proposed in Section 5.2 need to
be modified accordingly.

##  7 Application to Data from a Perinatal Study

We present an application of the methodology and algorithms proposed in
Section 5 to the National Collaborative Perinatal Project data (Klebanoff,
2009), a multisite prospective cohort study which took place in the United
States of America between 1959 and 1974. This study was designed to identify
the effects of complications during pregnancy or the perinatal period on birth
and child outcomes. The data are publicly available from the U.S. National
Archives with identifier 606622. Many online resources already employed this
dataset or subset of it for several analysis and over the years the dataset
has become a high-quality reference for biomedical and behavioral research in
many areas such as obstetrics, perinatology, pediatrics, and developmental
psychology.

The same data were examined in Nolan _et al._ (2020) and we expressly account
for the same model specification to experiment a suitable fixed effects
selection for a moderately large set of regressors that have been excluded
from their analysis and that may have a relevant impact onto explaining the
response variable. A full-blown analysis goes beyond the scope of this paper
and we focused on predicting the height-for-age z-score for 37,257 infants
followed longitudinally over their first year of life, following indications
from Taylor (1980). The height-for-age z-score is a standardized measure of
the World Health Organization for the height of children after accounting for
age; see World Health Organization (2006) for insights on how notable
discrepancies from this index standard reference values constitutes an alarm
signal for malnutrition symptoms.

We performed a Bayesian analysis that accounts for the heterogeneity of the
evolution of such index across infants. Our analysis was performed through a
two-level random effects linear model having a random intercept, and linear
and quadratic slopes ($q=3$) to account for the quadratic evolution of that
score over time. All the fixed effects regressors, excluding the intercept,
age of the infant (in days since birth) and its square, were subject to
selection; these include characteristics of the infant at birth (e.g. weight,
length, head circumference, sex and Apgar scores), and characteristics of the
mother, father and family. We also accounted for possible interactions between
some infant characteristics and sex, for a total of 38 candidate predictors
subject to selection.

The model we fitted respects the general specification (22) and can be
expressed for the generic $i$th infant as follows:

$\begin{array}[]{c}\boldsymbol{y}_{i}|\boldsymbol{\beta},\boldsymbol{u}_{i},\sigma^{2}\stackrel{{\scriptstyle{\tiny\mbox{ind.}}}}{{\sim}}\mbox{N}\left(\boldsymbol{X}^{\text{R}}_{i}\boldsymbol{\beta}^{\text{R}}+\boldsymbol{X}^{\text{S}}_{i}\boldsymbol{\beta}^{\text{S}}+\boldsymbol{Z}_{i}\boldsymbol{u}_{i},\sigma^{2}\boldsymbol{I}\right),\quad\boldsymbol{u}_{i}|\Sigma\stackrel{{\scriptstyle{\tiny\mbox{ind.}}}}{{\sim}}\mbox{N}_{3}(\boldsymbol{0},\boldsymbol{\Sigma}),\quad
1\leq i\leq 37,257\>,\\\\[6.45831pt]
\left[\begin{array}[]{c}\boldsymbol{\beta}^{\text{R}}\\\\[4.30554pt]
\boldsymbol{\beta}^{\text{S}}\end{array}\right]\Bigg{|}\>\boldsymbol{\zeta},\tau^{2}\sim\mbox{N}\left(\left[\begin{array}[]{c}\boldsymbol{0}_{3}\\\\[4.30554pt]
\boldsymbol{0}_{38}\end{array}\right],\left[\begin{array}[]{cc}10^{10}\boldsymbol{I}_{3}&\boldsymbol{O}\\\\[4.30554pt]
\boldsymbol{O}&\tau^{2}\>\mbox{diag}(\boldsymbol{\zeta})^{-1}\end{array}\right]\right),\\\\[15.0694pt]
\sigma^{2}|a_{\sigma^{2}}\sim\text{Inverse-}\chi^{2}(1,1/a_{\sigma^{2}}),\quad
a_{\sigma^{2}}\sim\text{Inverse-}\chi^{2}(1,10^{-10}),\\\\[6.45831pt]
\boldsymbol{\Sigma}|\boldsymbol{A}_{\mbox{\tiny{$\boldsymbol{\Sigma}$}}}\sim\text{Inverse-
G-Wishart}(G_{\mbox{\tiny
full}},6,\boldsymbol{A}_{\mbox{\tiny{$\boldsymbol{\Sigma}$}}}^{-1}),\\\\[6.45831pt]
\boldsymbol{A}_{\mbox{\tiny{$\boldsymbol{\Sigma}$}}}\sim\text{Inverse-G-
Wishart}(G_{\mbox{\tiny diag}},1,2\times
10^{-10}\boldsymbol{I}_{3}),\\\\[6.45831pt]
\zeta_{h}|a_{\zeta_{h}}\stackrel{{\scriptstyle{\tiny\mbox{ind.}}}}{{\sim}}\begin{cases}\>\mbox{Inverse-}\chi^{2}(2,1)\qquad&\mbox{for
a Laplace prior}\\\ \>\mbox{Gamma}(1/2,a_{\zeta_{h}})\qquad&\mbox{for a
Horseshoe prior}\\\
\>\mbox{Inverse-}\chi^{2}(2,2a_{\zeta_{h}})\qquad&\mbox{for a NEG
prior,}\end{cases}\\\\[17.22217pt]
a_{\zeta_{h}}\stackrel{{\scriptstyle{\tiny\mbox{ind.}}}}{{\sim}}\begin{cases}\qquad\quad-\qquad&\mbox{for
a Laplace prior}\\\ \mbox{Gamma}(1/2,1)\qquad&\mbox{for a Horseshoe prior}\\\
\mbox{Gamma}(\lambda,1)\qquad&\mbox{for a NEG
prior,}\end{cases}\\\\[17.22217pt] \text{for }1\leq h\leq 38,\\\\[4.30554pt]
\tau^{2}|a_{\tau^{2}}\sim\mbox{Inverse-}\chi^{2}(1,1/a_{\tau^{2}}),\quad
a_{\tau^{2}}\sim\mbox{Inverse-}\chi^{2}(1,10^{-10}).\end{array}$ (29)

For the $i$th infant, $o_{i}$ time-point measurements were recorded, ranging
from one to four in number. The $\boldsymbol{X}_{i}^{\text{R}}$ matrix has
size $n_{i}\times 3$ with the first column being a vector of ones, the second
one consists of the time-point measurements for the $i$th infant and the third
column containing the square of the elements of the second one. We set
$\boldsymbol{X}_{i}^{\text{A}}=\boldsymbol{O}$, while the
$\boldsymbol{X}_{i}^{\text{S}}$ matrix of size $n_{i}\times 38$ consisted of
all the considered predictors subject to selection. Moreover,
$\boldsymbol{Z}_{i}=\boldsymbol{X}_{i}^{\text{R}}$ by definition.
Uninformative priors were placed over all the model parameters. We fitted the
model using the three global-local priors treated in this work, and the
Gaussian prior for $\boldsymbol{\beta}^{\text{S}}$ treated in Nolan _et al._
(2020).

Streamlined MFVB and MCMC were used for model fitting. The former was
performed running Algorithm 2 and stopping it after 200 iterations, while the
latter was performed running 25,000 iterations to which a thinning factor of 5
was applied after discarding 5,000 burnin iterations. The whole input data
required approximately 200 megabytes of memory storage, while a naïve MFVB
procedure would necessitate several gigabytes of memory to entirely store the
$\boldsymbol{Z}$ matrix, that is composed by 11,693,705,562 cells of which the
99.997% are zeros. All the covariates excluding the binary ones were
standardized and the estimates were rescaled back to the original scale before
presenting the results. The streamlined MFVB algorithms took 2 to 3 minutes to
run for each prior specification, while the associated MCMC samplers took more
than 35 minutes.

Figure 4: The $90\%$ high posterior density credible intervals for the fixed
effects subject to selection in the real data application. The four different
priors for $\boldsymbol{\beta}^{\text{S}}$ are represented with different
colors. For each fixed effect, the thicker lines correspond to the intervals
obtained from the MCMC approximate marginal posterior densities, while the
thinner lines represent those obtained from the streamlined MFVB approximating
densities.

We omit the presentation of model interpretation, goodness-of-fit analysis,
accuracy of the approximations, convergence of the MCMC chains and
visualization of the fitted height-for-age z-score trajectories over time.
Instead, in Figure 4 we present 90% high posterior density credible intervals
for all the fixed effects subject to selection. The thicker lines represent
the intervals obtained from the MCMC posterior samples through the emp.hpd
function of the TeachingDemos package (Snow, 2020), while the thinner lines
represent those calculated from the determined optimal MFVB approximate
posteriors. The two lines are superimposed to facilitate immediate comparison,
for each different prior specification.

Overall, MFVB provided very accurate high posterior density credible intervals
when compared to MCMC. The MCMC chain associated to the birth length fixed
effect showed some convergence problems, probably due to its moderate
correlation with the response variable and this reduced the overlap with MFVB.
Importantly, the Laplace global-local prior produced results very similar to
those provided by the Gaussian prior specification. On the other hand, both
the Horseshoe and Normal Exponential Gamma priors seem to have effectively
shrunk most of the fixed effects towards zero, especially those associated to
dummy variables.

The SAVS procedure applied to the approximate posterior densities classifies
the following fixed effects as relevant: indicator that the infant is male,
infant’s year of birth, birth length, birth head circumference, the Apgar
score assessed 5 minutes after the infant’s birth and the mother height. Some
effects such as the ethnicity of the mother, the place where she conceived,
her marital status and informations on previous pregnancies were not indicated
as relevant by the SAVS procedure, albeit some of their associated high
posterior density credible intervals are far from zero. Notice also that, for
some fixed effects, global-local priors drag the intervals towards the origin,
with different intensities probably depending on the variability and
correlation of the covariates. More ad-hoc analyses are firmly suggested,
although these go beyond the scope of the current work.

##  8 Conclusions

In this work, we developed streamlined mean-field variational Bayes procedures
for Gaussian response linear mixed models having nested random effects
structures and admitting fixed effects prior specifications alternative to the
classical Gaussian one. The priors we considered are amenable to automated and
hyperparameter-free fixed effects selection procedures. Simulated and real
data examples showed how streamlined variational inference can provide
impressive benefits in terms of computational time and memory saving when
compared to inefficient implementations of variational approximations. Albeit
the marginal posterior densities of fixed effects subject to selection are
approximated by bell-shaped curves, our studies showed high performances of
the automated selection procedure. It is also worth mentioning that the more
restrictive mean-field approximations discussed at the end of Section 6 can
sensibly improve the benefit of streamlined variational inference for large
$p$ scenarios and, in general, when speed is more important than approximation
accuracy.

Numerous ramifications of the proposed methodology can be envisaged. These
include the treatment of models with unit-specific errors, heteroskedastic
covariance structures for groups and sub-groups, higher levels of nesting or
crossed random effects and the extension to generalized linear mixed models.
Ormerod and Wand (2010), Wand et al. (2011), Nolan and Wand (2017), Maestrini
& Wand (2018) and McLean & Wand (2019) provide variational inference
algorithms for models with a variety of response distributions such as
Bernoulli with probit or logit link, Poisson, Negative Binomial, $t$,
Asymmetric Laplace, Skew Normal, Skew $t$ and Finite Normal Mixtures. The
variational algorithms presented in these references have been derived using
data-augmentation representations of the response distributions. The same
strategy can be adopted to derive MFVB algorithms for non-Gaussian response
mixed models and allow for the implementation of streamlined variational
inference through a relatively straightforward adaptation of the algorithms
presented in this work. Maestrini (2019) also provides explicit streamlined
variational inference algorithms for two-level mixed models with binary-
logistic and Poisson responses, and generic Gaussian priors for the fixed
effects coefficients.

Regarding the selection of fixed effects, alternative shrinkage priors
belonging to the global-local family can be accounted for with minor
modifications to the proposed algorithms. Following the prescriptions of
Neville _et al._ (2014), it is possible to study global-local prior
specifications that are not based on auxiliary variables and investigate
whether this could lead to better approximation accuracies, although at the
cost of more computationally intensive variational updates. We also mention
the possibility of admitting spike-and-slab prior formulations following, for
example, a variational inference approach similar to that of Carbonetto &
Stephens (2012). The drawback is that spike-and-slab priors require more
involved algebra for implementing streamlined variational inference. Using the
Bayesian Adaptive Lasso (BaLasso) prior proposed by Leng _et al._ (2014) it is
also possible to extend our individual fixed-effects selection approach to
account for (ordered) group selection by imposing different shrinkage levels
to different coefficients. A variational Bayes approach for generalized linear
models with priors of this type has been explored by Tung _et al._ (2019),
where the variational parameter updates of their VBGLMM algorithm could be
streamlined using the framework studied in our work.

One last possible direction to explore concerns streamlined variational
inference for models with priors for selecting random effects. Consider for
simplicity the two-level random effects case in which
$\boldsymbol{u}_{i}\stackrel{{\scriptstyle{\tiny\mbox{ind.}}}}{{\sim}}\mbox{N}(\boldsymbol{0},\boldsymbol{\Sigma})$
for $1\leq i\leq m$. If the covariance matrix $\boldsymbol{\Sigma}$ is
supposed to be diagonal, i.e. $\mbox{Cov}(u_{ih},u_{ik})=0$ for each $1\leq
h,k\leq q$ with $h\neq k$ and $1\leq i\leq m$, and a global-local prior
distribution is imposed to each random effect, then explicit streamlined MFVB
updates can be obtained with straightforward manipulation of the results
presented in the supplementary material and references therein. This applies
to three-level models in a similar way. However, assuming that
$\boldsymbol{\Sigma}$ is diagonal may sometimes be quite restrictive and the
highest computational advantage of streamlined MFVB over its naïve counterpart
is achieved when the random effects covariance matrix is non-diagonal.
Furthermore, it is typically easier to provide an interpretation to the
selection of fixed effects than explaining why certain random effects are
relevant only for a subset of clusters or sub-clusters. For these reasons the
present work provides a first extension of the streamlined MFVB approach of
Nolan _et al._ (2020) focusing on fixed-effects selection. Another relevant
extension is the derivation of (streamlined) MFVB algorithms with random
effects selection procedures admitting non-diagonal structures for
$\boldsymbol{\Sigma}$. The covariance matrix decomposition approach of Chen &
Dunson (2003) and related methods discussed in Section 1 can be exploited for
this extension.

## Acknowledgments

The research presented in this article was supported by the Australian
Research Council Discovery Project DP180100597 and the Australian Research
Council Center of Excellence for Mathematical and Statistical Frontiers (grant
CE140100049).

## References

Andrews, D. F. & Mallows, C. L. (1974). Scale mixtures of normal
distributions. _Journal of the Royal Statistical Society B_ , 36, 99–102.

Armagan, A. & Dunson, D. B. (2011). Sparse variational analysis of linear
mixed models for large data sets. _Statistics & Probability Letters_, 81,
1056–1062.

Armagan, A., Dunson, D. B. & Lee, J. (2013). Generalized double Pareto
shrinkage. _Statistica Sinica_ , 23, 119–143.

Atay-Kayis, A. & Massam, H. (2005). A Monte Carlo method for computing
marginal likelihood in nondecomposable Gaussian graphical models. _Biometrika_
, 92, 317–335.

Baltagi, B. H. (2013). _Econometric Analysis of Panel Data, Fifth Edition_.
Chichester, U.K.: John Wiley & Sons.

Barbieri, M. M. & Berger, J. O. (2004). Optimal predictive model selection.
_The Annals of Statistics_ , 32, 870–897.

Bhadra, A., Datta, J., Polson, N. G. & Willard, B. (2017). The horseshoe+
estimator of ultra-sparse signals. _Bayesian Analysis_ , 12, 1105–1131.

Bhadra, A., Datta, J., Polson, N. G. & Willard, B. (2019). Lasso meets
horseshoe: a survey. _Statistical Science_ , 34, 405–427.

Bhattacharya, A., Chakraborty, A. & Mallick, B. K. (2016). Fast sampling with
Gaussian scale mixture priors in high-dimensional regression. _Biometrika_ ,
103, 985–991.

Bhattacharya, A., Pati, D., Pillai, N. S. & Dunson, D. B. (2015). Dirichlet-
Laplace priors for optimal shrinkage. _Journal of the American Statistical
Association_ , 110, 1479–1490.

Bishop, C. M. (2006). _Pattern Recognition and Machine Learning_. New York:
Springer.

Blei, D. M., Kucukelbir, A. & McAuliffe, J. D. (2017). Variational inference:
a review for statisticians. _Journal of the American Statistical Association_
, 112, 859–877.

Bogdan, M., Chakrabarti, A., Frommlet, F. & Ghosh, J. K. (2011). Asymptotic
Bayes-optimality under sparsity of some multiple testing procedures. _The
Annals of Statistics_ , 39, 1551–1579.

Bondell, H. D. & Reich, B. J. (2012). Consistent high-dimensional Bayesian
variable selection via penalized credible regions. _Journal of the American
Statistical Association_ , 107, 1610–1624.

Boyd, S. & Vandenberghe, L. (2004). _Convex Optimization_ , Cambridge, U.K.:
Cambridge University Press.

Brown, H. & Prescott, R. (2014). _Applied Mixed Models in Medicine, Third
Edition_ , Chichester, U.K.: John Wiley & Sons.

Bürkner, P.-C. (2017). brms: an R package for Bayesian multilevel models using
Stan. _Journal of Statistical Software_ , 80, 1–28.

Carbonetto, P. & Stephens, M. (2012). Scalable variational inference for
Bayesian variable selection in regression, and its accuracy in genetic
association studies. _Bayesian Analysis_ , 7, 73–108.

Carpenter, B., Gelman, A., Hoffman, M., Lee, D., Goodrich, B., Betancourt, M.,
Brubaker, M., Guo, J., Li, P. & Riddell, A. (2017). Stan: a probabilistic
programming language. _Journal of Statistical Software_ , 76, 1–32.

Carvalho, C. M., Polson, N. G. & Scott, J. G. (2009). Handling sparsity via
the horseshoe. _Journal of Machine Learning Research_ , 5, 73–80.

Carvalho, C. M., Polson, N. G. & Scott, J. G. (2010). The horseshoe estimator
for sparse signals. _Biometrika_ , 97, 465–480.

Chen, Z. & Dunson, D. B. (2003). Random effects selection in linear mixed
models. _Biometrics_ , 59, 762–769.

Eddelbuettel, D. & Sanderson, C. (2014). RcppArmadillo: accelerating R with
high-performance C++ linear algebra. _Computational Statistics and Data
Analysis_ , 71, 1054–1063.

Efron, B. (2008). Microarrays, empirical Bayes and the two-groups model.
_Statistical Science_ , 23, 1–22.

Faes, C., Ormerod, J. & Wand, M.P. (2011). Variational Bayesian inference for
parametric and nonparametric regression with missing data. _Journal of the
American Statistical Association_ , 106, 959–971.

Fan, Y. & Li, R. (2012). Variable selection in linear mixed effects models.
_The Annals of Statistics_ , 40, 2043–2068.

Fitzmaurice, G. , Davidian, M., Verbeke, G. & Molenberghs, G. (2008).
_Longitudinal Data Analysis_. Boca Raton, Florida: Chapman & Hall/CRC.

Frank, I. E. & Friedman, J. H. (1993). A statistical view of some chemometrics
regression tools. _Technometrics_ , 35, 109–148.

Gelman, A. (2006). Prior distributions for variance parameters in hierarchical
models. _Bayesian Analysis_ , 1, 515–533.

George, E. I. & McCulloch, R. E. (1997). Approaches for Bayesian variable
selection. _Statistica Sinica_ , 7, 339–373.

Goldstein, H. (2010). _Multilevel Statistical Models, Fourth Edition_.
Chichester, U.K.: John Wiley & Sons.

Griffin, J. E. & Brown, P. J. (2010). Inference with normal-gamma prior
distributions in regression problems. _Bayesian Analysis_ , 5, 171–188.

Griffin, J. E. & Brown, P. J. (2011). Bayesian hyper-lassos with non-convex
penalization. _Australian and New Zealand Journal of Statistics_ , 53,
423–442.

Groll, A. & Tutz, G. (2012). Variable selection for generalized linear mixed
models by $l_{1}$-penalized estimation. _Statistics and Computing_ , 24,
137–154.

Hahn, P. R. & Carvalho, C. M. (2015). Decoupling shrinkage and selection in
Bayesian linear models: a posterior summary perspective. _Journal of the
American Statistical Association_ , 110, 435–448.

Hoerl, A. E. & Kennard, R. W. (1970). Ridge regression: biased estimation for
nonorthogonal problems. _Technometrics_ , 12, 55–67.

Huang, A. & Wand, M. P. (2013). Simple marginally noninformative prior
distributions for covariance matrices. _Bayesian Analysis_ , 8, 439–452.

Hughes, D. M., García-Finaña, M. & Wand, M. P. (2021). Fast approximate
inference for multivariate longitudinal data. _Biostatistics_ (volume and page
numbers pending).

Hui, F. K. C., Müller, S. & Welsh, A. H. (2017). Joint selection in mixed
models using regularized PQL. _Journal of the American Statistical
Association_ , 112, 1323–1333.

Ishwaran, H. & Rao, J. S. (2005). Spike and slab variable selection:
frequentist and Bayesian strategies. _The Annals of Statistics_ , 33, 730–773.

Johnstone, I. M. & Silverman, B. W. (2005). Bayes selection of wavelet
thresholds. _The Annals of Statistics_ , 33, 1700–1752.

Kinney, S. K. & Dunson, D. B. (2007). Fixed and random effects selection in
linear and logistic models. _Biometrics_ , 63, 690–698.

Klebanoff M. A. (2009). The collaborative perinatal project: a 50-year
retrospective. _Paediatric and Perinatal Epidemiology_ , 23, 2–8.

Korte, A., Vilhjálmsson, B., Segura, V., Platt, A., Long, Q. & Nordborg, M.
(2012). A mixed-model approach for genome-wide association studies of
correlated traits in structured populations. _Nature Genetics_ , 44,
1066–1071.

Lee, C. Y. Y. & Wand, M. P. (2016). Streamlined mean field variational Bayes
for longitudinal and multilevel data analysis. _Biometrical Journal_ 58,
868–895.

Leng, C., Tran, M.N. & Nott, D. (2014). Bayesian adaptive Lasso. _Annals of
the Institute of Statistical Mathematics_ , 66, 221–244.

Li, H. & Pati, D. (2017). Variable selection using shrinkage priors.
_Computational Statistics and Data Analysis_ , 107, 107–119.

Li, J., Wang, Z., Li, R. & Wu, R. (2015). Bayesian group Lasso for
nonparametric varying-coefficient models with application to functional
genome-wide association studies. _The Annals of Applied Statistics_ , 9,
640–664.

Lindner, C. C. & Rodger, C. A. (1997). _Design Theory_. Boca Raton, Florida:
Chapman & Hall/CRC.

Li, Y., Wang, S., Song, P. X.–K., Wang, N., Zhou, L. & Zhu, J. (2018). Doubly
regularized estimation and selection in linear mixed-effects models for high-
dimensional longitudinal data. _Statistics and Its Interface_ , 11, 721–737.

Luts, J., Broderick, T. & Wand, M. P. (2014). Real-time semiparametric
regression. _Journal of Computational and Graphical Statistics_ , 23, 589–615.

Maestrini, L. & Wand, M.P. (2018). Variational message passing for skew t
regression. _Stat_ , 7, 1–11.

Maestrini, L. (2019). _On Variational Approximations for Frequentist and
Bayesian Inference_. Doctor of Philosophy thesis, Università degli Studi di
Padova, Italy.

Maestrini, L. & Wand, M. P. (2021). The inverse G-Wishart distribution and
variational message passing. _Australian and New Zealand Journal of
Statistics_ , 63, 517–541.

McLean, M.W. & Wand, M.P. (2019). Variational message passing for elaborate
response regression models. Bayesian Analysis, 14, 371–398.

Menictas, M., Di Credico, G. & Wand, M. P. (2022). Streamlined variational
inference for linear mixed models with crossed random effects. _Journal of
Computational and Graphical Statistics_ (volume and page numbers pending).

Menictas, M., Nolan, T.H., Simpson, D. G. & Wand, M. P. (2021). Streamlined
variational inference for higher level group-specific curve models.
_Statistical Modelling_ , 21, 479–519.

Minka, T., Winn, J., Guiver, J., Zaykov, Y., Fabian, D. & Bronskill, J.
(2018). Infer.NET 0.3, Microsoft Research Cambridge, Cambridge, U.K.
_http://dotnet.github.io/infer_.

Mitchell, T. J. & Beauchamp, J. J. (1988). Bayesian variable selection in
linear regression. _Journal of the American Statistical Association_ , 83,
1023–1032.

Neville, S. E., Ormerod, J. T & Wand, M. P. (2014). Mean field variational
Bayes for continuous sparse signal shrinkage: pitfalls and remedies.
_Electronic Journal of Statistics_ , 8, 1113–1151.

Nolan, T. H., Menictas, M. & Wand, M. P. (2020). Streamlined computing for
variational inference with higher level random effects. _Journal of Machine
Learning Research_ , 21, 1–62.

Nolan, T.H. & Wand, M.P. (2017). Accurate logistic variational message
passing: algebraic and numerical details. _Stat_ , 6, 102–112.

Nolan, T. H. & Wand, M. P. (2020). Streamlined solutions to multilevel sparse
matrix problems. _ANZIAM Journal_ , 62, 18–41.

O’Hara, R. B. & Sillanpää, M. J. (2009). A review of Bayesian variable
selection methods: what, how and which. _Bayesian Analysis_ , 4, 85–118.

Ormerod, J. T. & Wand, M. P. (2010). Explaining variational approximations.
_The American Statistician_ , 64, 140–153.

Ormerod, J. T., You, C. & Müller, S. (2017). A variational Bayes approach to
variable selection. _Electronic Journal of Statistics_ , 11, 3549–3594.

Park, T. & Casella, G. (2008). The Bayesian lasso. _Journal of the American
Statistical Association_ , 103, 681–686.

Pinheiro, J. C. & Bates, D. M. (2006). _Mixed-Effects Models in_ S and S-PLUS.
New York: Springer Science & Business Media.

Polson, N. G. & Scott, J. G. (2011). Shrink globally, act locally: sparse
Bayesian regularization and prediction. _Bayesian Statistics_ , 9, 501–538.

Rao, J. N. K. & Molina, I. (2015). _Small Area Estimation, Second Edition_.
Hoboken, New Jersey: John Wiley & Sons.

Ray, P., & Bhattacharya, A. (2018). Signal Adaptive Variable Selector for the
horseshoe prior. Unpublished manuscript available at
_https://arxiv.org/abs/1810.09004_.

R Core Team. (2022). R: A language and environment for statistical computing.
R Foundation for Statistical Computing, Vienna, Austria.
_https://www.R-project.org/_.

Ruppert, D., Wand, M. P. & Carroll, R. (2003). _Semiparametric Regression_.
Cambridge, U.K.: Cambridge University Press.

Sanderson, C. & Curtin, R. (2016). Armadillo: a template-based C++ library for
linear algebra. _Journal of Open Source Software_ , 1, 26.

Schelldorger, J., Bühlmann, P. & Van De Geer, S. (2011). Estimation for high-
dimensional linear mixed-effects models using $l_{1}$-penalization.
_Scandinavian Journal of Statistics_ , 38, 197–214.

Sikorska, K., Rivadeneira, F., Groenen, P. J., Hofman, A., Uitterlinden, A.
G., Eilers, P. H. & Lesaffre, E. (2012). Fast linear mixed model computations
for genome-wide association studies with longitudinal data. _Statistics in
Medicine_ , 32, 165–180.

Snow, G. (2020). TeachingDemos: demonstrations for teaching and learning. R
package version 2.12. _https://CRAN.R-project.org/package=TeachingDemos_.

Tang, X., Xu, X., Ghosh, M. & Ghosh, P. (2018). Bayesian variable selection
and estimation based on global-local shrinkage priors. _Sankhya A_ , 80,
215–246.

Taylor, P. M. (1980). The first year of life: the collaborative perinatal
project of the National Institute of Neurological and Communicative Disorders
and Stroke. _Journal of the American Medical Association_ , 244, 1503.

Tibshirani, R. (1996). Regression shrinkage and selection via the lasso.
_Journal of the Royal Statistical Society B_ , 58, 267–288.

Tung, D. T., Tran, M.-N. & Cuong, T. M. (2019). Bayesian adaptive lasso with
variational Bayes for variable selection in high-dimensional generalized
linear mixed models. _Communications in Statistics - Simulation and
Computation_ , 48, 530–543.

Van Rijsbergen, C. J. (1979). _Information Retrieval (2nd ed.)_. London:
Butterworths.

Verbeke, G. & Molenberghs, G. (2000). _Linear Mixed Models for Longitudinal
Data_. New York: Springer.

Vonesh, E. F. & Chinchilli, V. G. (1997). _Linear and Nonlinear Models for the
Analysis of Repeated Measurements_. London: Chapman and Hall.

Wand, M. P. (2020). KernSmooth: functions for kernel smoothing supporting Wand
& Jones (1995). R package version 2.23-18.
_https://CRAN.R-project.org/package=KernSmooth_.

Wand, M. P. & Jones, M.C. (1995). _Kernel Smoothing_. London: Chapman and
Hall.

Wand, M. P., Ormerod, J. T., Padoan, S. A., & Frühwirth, R. (2011). Mean field
variational Bayes for elaborate distributions. _Bayesian Analysis_ , 6,
847–900.

Wand, M. P. & Ormerod, J. T. (2011). Penalized wavelets: embedding wavelets
into semiparametric regression. _Electronic Journal of Statistics_ , 5,
1654–1717.

Wang, S. S. J. & Wand, M. P. (2011). Using Infer.NET for statistical analyses.
_The American Statistician_ , 65, 115–126.

West, M. (1987). On scale mixtures of normal distributions. _Biometrika_ , 74,
646–648.

Whittaker, E. T. & Watson, G. N. (1990). _A Course in Modern Analysis_.
Cambridge, U.K.: Cambridge University Press.

World Health Organization (2006). WHO child growth standards: length/height-
for-age, weight-for-age, weight-for-length, weight-for-height and body mass
index-for-age: methods and development. Available at
_https://apps.who.int/iris/handle/10665/43413_.

Yang, M. (2013). Bayesian nonparametric centered random effects models with
variable selection. _Biometrical Journal_ , 55, 217–230.

Yang, M., Wang, M. & Dong, G. (2020). Bayesian variable selection for mixed
effects model with shrinkage prior. _Computational Statistics_ , 35, 227–243.

Zhang, Y. & Bondell, H. D. (2018). Variable selection via penalized credible
regions with Dirichlet-Laplace global-local shrinkage priors. _Bayesian
Analysis_ , 13, 823–844.

Zou, H. (2006). The adaptive lasso and its oracle properties. _Journal of the
American Statistical Association_ , 101, 1418–1429.

Zou, H. & Hastie, T. (2005). Regularization and variable selection via the
elastic net. _Journal of the Royal Statistical Society Series B_ , 67,
301–320.

Supplementary Material for:

Sparse Linear Mixed Model Selection via

Streamlined Variational Bayes

By Emanuele Degani†,♮, Luca Maestrini‡,

Dorota Toczydłowska♯ & Matt P. Wand♯

Università degli Studi di Padova†, Banca d’Italia – Eurosystem♮,

The Australian National University‡ and University of Technology Sydney♯

##  Appendix S.1 Distributions and Associated Useful Results

Many probability distributions are used throughout the paper. We provide
details on their probability density functions together with additional useful
results.

###  S.1.1 Inverse-Gaussian Distribution

A continuous random variable $x$ has an Inverse Gaussian distribution with
mean parameter $\mu>0$ and rate parameter $\lambda>0$, written
$x\sim\mbox{Inverse-Gaussian}(\mu,\lambda)$, if the density function of $x$ is

$\mathfrak{p}(x)=\lambda^{1/2}(2\pi
x^{3})^{-1/2}\exp\left\\{-\frac{\lambda(x-\mu)^{2}}{2\mu^{2}x}\right\\},\quad
x>0.$

###  S.1.2 Laplace Distribution

A continuous random variable $x$ has a Laplace distribution (also known as
_Double-Exponential_ or _two-tailed Exponential_ distribution) with mean
parameter $\mu$ and scale parameter $\sigma>0$, written
$x\sim\mbox{Laplace}(\mu,\sigma)$, if the density function of $x$ is

$\mathfrak{p}(x)=\frac{1}{2\sigma}\exp\left(-\frac{|x-\mu|}{\sigma}\right),\quad
x\in\mathbb{R}.$

A useful result from Andrews & Mallows (1974) and West (1987) shows that

$\begin{array}[]{c}\mbox{if}\quad
x|b\sim\text{N}(\mu,\sigma^{2}/b)\quad\mbox{and}\quad b\sim\mbox{Inverse-
Gamma}(1,1/2),\\\\[4.30554pt] \mbox{then}\quad
x\sim\mbox{Laplace}(\mu,\sigma).\end{array}$

###  S.1.3 Horseshoe Distribution

A continuous random variable $x$ has a Horseshoe distribution with mean
parameter $\mu$ and scale parameter $\sigma>0$, written
$x\sim\mbox{Horseshoe}(\mu,\sigma)$, if the density function of $x$ is

$\mathfrak{p}(x)=(2\pi^{3})^{-1/2}\sigma^{-1}\exp\left\\{\frac{(x-\mu)^{2}}{2\sigma^{2}}\right\\}E_{1}\left\\{\frac{(x-\mu)^{2}}{2\sigma^{2}}\right\\},\quad
x\in\mathbb{R}$

where $E_{1}(x)\equiv\int_{x}^{\infty}t^{-1}e^{-t}\mathrm{d}t$, with $x\neq
0$, is the exponential integral function of order 1. A useful result from
Carvalho _et al._ (2010) shows that

$\begin{array}[]{c}\mbox{if}\quad x|b\sim\text{N}(\mu,\sigma^{2}/b),\quad
b|c\sim\mbox{Gamma}(1/2,c)\quad\mbox{and}\quad
c\sim\mbox{Gamma}(1/2,1),\\\\[4.30554pt] \mbox{then}\quad
x\sim\mbox{Horseshoe}(\mu,\sigma).\end{array}$

###  S.1.4 Normal-Exponential-Gamma Distribution

A continuous random variable $x$ has a Normal-Exponential-Gamma distribution
with mean parameter $\mu$, scale parameter $\sigma>0$ and shape parameter
$\lambda>0$, written $x\sim\mbox{NEG}(\mu,\sigma,\lambda)$, if the density
function of $x$ is

$\mathfrak{p}(x)=\pi^{-1/2}\sigma^{-1}\lambda
2^{\lambda}\Gamma(\lambda+1/2)\exp\left\\{\frac{(x-\mu)^{2}}{4\sigma^{2}}\right\\}D_{-2\lambda-1}\left(\left|\frac{x-\mu}{\sigma}\right|\right),\quad
x\in\mathbb{R},$

where $D_{\nu}(x)\equiv
2^{\nu/2+1/4}W_{\nu/2+1/4,-1/4}(x^{2}/2)/\sqrt{x},\>x>0$ is the parabolic
cylinder function of order $\nu\in\mathbb{R}$ and $W_{k,m}$ is a confluent
hypergeometric function or order $k$ and $m$, as defined by Whittaker & Watson
(1990). A useful result from Griffin & Brown (2011) shows that

$\begin{array}[]{c}\mbox{if}\quad x|b\sim\text{N}(\mu,\sigma^{2}/b),\quad
b|c\sim\mbox{Inverse-Gamma}(1,c)\quad\mbox{and}\quad
c\sim\mbox{Gamma}(\lambda/2,1),\\\\[4.30554pt] \mbox{then}\quad
x\sim\mbox{NEG}(\mu,\sigma,\lambda).\end{array}$

###  S.1.5 Gamma Distribution

A continuous random variable $x$ has a Gamma distribution with shape parameter
$\alpha>0$ and scale parameter $\beta>0$, written
$x\sim\mbox{Gamma}(\alpha,\beta)$, if the density function of $x$ is

$\mathfrak{p}(x)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}\exp\\{-\beta
x\\},\quad x>0.$

###  S.1.6 Inverse-$\chi^{2}$ and Inverse-Gamma Distributions

A continuous random variable $x$ has an Inverse-$\chi^{2}$ distribution with
shape parameter $\xi>0$ and scale parameter $\lambda>0$, written
$x\sim\mbox{Inverse-}\chi^{2}(\xi,\lambda)$, if the density function of $x$ is

$\mathfrak{p}(x)=\frac{\left(\lambda/2\right)^{\xi/2}}{\Gamma(\xi/2)}x^{-(\xi/2)-1}\exp\\{-(\lambda/2)/x\\},\quad
x>0.$

A continuous random variable $x$ has an Inverse-Gamma distribution with shape
parameter $\alpha>0$ and scale parameter $\beta>0$, written
$x\sim\mbox{Inverse-Gamma}(\alpha,\beta)$, if the density function of $x$ is

$\mathfrak{p}(x)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{-\alpha-1}\exp\\{-\beta/x\\},\quad
x>0.$

Notice that

$\quad x\sim\mbox{Inverse-}\chi^{2}(\xi,\lambda)\quad\mbox{if and only
if}\quad x\sim\mbox{Inverse-Gamma}(\xi/2,\lambda/2)$

and, equivalently,

$\quad x\sim\mbox{Inverse-Gamma}(\alpha,\beta)\quad\mbox{if and only if}\quad
x\sim\mbox{Inverse-}\chi^{2}(2\alpha,2\beta).$

###  S.1.7 Half-Cauchy Distribution

A continuous random variable $x$ has a Half-Cauchy distribution with scale
parameter $\sigma>0$, written $x\sim\mbox{Half-Cauchy}(\sigma)$, if the
density function of $x$ is

$p(x)=2/[\sigma\pi\\{1+(x/\sigma)^{2}\\}],\quad x>0.$

###  S.1.8 Half-t Distribution

A continuous random variable $x$ has a $\mbox{Half-}t$ distribution with
$\nu>0$ degrees of freedom and scale parameter $\sigma>0$, written
$x\sim\mbox{Half-}t(\sigma,\nu)$, if the density function of $x$ is

$\mathfrak{p}(x)=\frac{2\Gamma(\frac{\nu+1}{2})}{\sqrt{\pi\nu}\Gamma(\nu/2)\sigma\\{1+(x/\sigma)^{2}/\nu\\}^{\frac{\nu+1}{2}}},\quad
x>0.$

If $\nu=1$, then $x\sim\mbox{Half-Cauchy}(\sigma)$. Wand _et al._ (2011) show
that

$\begin{array}[]{c}\mbox{if}\quad x|a\sim\mbox{Inverse-
Gamma}(\nu/2,\nu/a)\quad\mbox{and}\quad a\sim\mbox{Inverse-
Gamma}(1/2,1/A^{2}),\\\\[4.30554pt]
\mbox{then}\quad\sqrt{x}\sim\mbox{Half-}t(A,\nu).\end{array}$

The same result can be equivalently formulated with:

$x|a\sim\mbox{Inverse-Gamma}(\nu/2,\nu/(2a))\quad\mbox{and}\quad
a\sim\mbox{Inverse-Gamma}(1/2,1/(2A^{2})),$

or with:

$x|a\sim\mbox{Inverse-}\chi^{2}(\nu,1/a)\quad\mbox{and}\quad
a\sim\mbox{Inverse-}\chi^{2}(1,1/(\nu A^{2})).$

###  S.1.9 Inverse-G-Wishart Distribution

The Inverse-G-Wishart distribution arises from the inverse of matrices with a
_G-Wishart_ distribution (e.g. Atay-Kayis & Massam, 2005; Maestrini & Wand,
2021). For any positive integer $d$, let $G$ be an undirected graph with $d$
nodes labeled $1,\ldots,d$ and set $E$ consisting of sets of pairs of nodes
that are connected by an edge. We say that the symmetric $d\times d$ matrix
$\boldsymbol{M}$ _respects_ $G$ if

$\boldsymbol{M}_{ij}=0\quad\mbox{for all}\quad\\{i,j\\}\notin E.$

A $d\times d$ random matrix $\boldsymbol{X}$ has an Inverse G-Wishart
distribution with graph $G$ and parameters $\xi>0$ and symmetric $d\times d$
matrix $\boldsymbol{\Lambda}$, written $\boldsymbol{X}\sim\mbox{Inverse-G-
Wishart}(G,\xi,\boldsymbol{\Lambda})$, if and only if the density function of
$\boldsymbol{X}$ satisfies

$\mathfrak{p}(\boldsymbol{X})\propto|\boldsymbol{X}|^{-(\xi+2)/2}\exp\left\\{-\frac{1}{2}\mbox{tr}(\boldsymbol{\Lambda}\,\boldsymbol{X}^{-1})\right\\}$

over arguments $\boldsymbol{X}$ such that $\boldsymbol{X}$ is symmetric and
positive definite and $\boldsymbol{X}^{-1}$ respects $G$. Two important
special cases are

$G=G_{\mbox{\tiny full}}\equiv\mbox{totally connected $d$-node graph},$

for which the Inverse G-Wishart distribution coincides with the ordinary
Inverse Wishart distribution $\boldsymbol{X}\sim\mbox{Inverse-
Wishart}(\xi-d+1,\boldsymbol{\Lambda})$, and

$G=G_{\mbox{\tiny diag}}\equiv\mbox{totally disconnected $d$-node graph},$

for which the Inverse G-Wishart distribution coincides with a product of
independent Inverse Chi-Squared distributions. The subscripts of
$G_{\mbox{\tiny full}}$ and $G_{\mbox{\tiny diag}}$ reflect the fact that
$\boldsymbol{X}^{-1}$ is a full matrix and $\boldsymbol{X}^{-1}$ is a diagonal
matrix in each special case. In the $d=1$ special case, the Inverse G-Wishart
distribution coincides with the Inverse Chi-Squared distribution.

###  S.1.10 Huang-Wand Distribution

A $d\times d$ symmetric positive definite matrix $\boldsymbol{\Sigma}$ has a
Huang-Wand distribution with $\nu>0$ degrees of freedom and scale parameters
$s_{1},\ldots,s_{d}$, written $\boldsymbol{\Sigma}\sim\mbox{Huang-
Wand}(\nu;s_{1},\ldots,s_{d})$, if and only if the following augmented
representation holds:

$\begin{array}[]{c}\boldsymbol{\Sigma}|\boldsymbol{A}\sim\mbox{Inverse-G-
Wishart}(G_{\mbox{\tiny full}},\nu+2d-2,\boldsymbol{A}^{-1})\\\\[4.30554pt]
\boldsymbol{A}\sim\mbox{Inverse-G-Wishart}(G_{\mbox{\tiny
diag}},1,\\{\nu\,\mbox{diag}(s_{1}^{2},\ldots,s_{d}^{2})\\}^{-1}).\end{array}$

This distribution is defined in Huang & Wand, 2013). The formulation displayed
here is provided in Maestrini & Wand (2021).

If $\nu=2$, a marginally noninformative Huang-Wand prior specification for
$\boldsymbol{\Sigma}$ can be obtained for arbitrarily large scale parameters.
This corresponds to the standard deviation parameters
$\sigma_{j}\equiv(\boldsymbol{\Sigma})_{jj}^{1/2},1\leq j\leq d$, having
$\mbox{Half-}t$ distributions with $\nu$ degrees of freedom and scale
parameter given by $s_{j}$, and the correlation parameters
$\rho_{jj^{\prime}}\equiv(\boldsymbol{\Sigma})_{jj^{\prime}}^{1/2}/(\sigma_{j}\sigma_{j^{\prime}})$,
$1\leq j,j^{\prime}\leq d$, having a Uniform distribution on the interval
$(-1,1)$.

##  Appendix S.2 Multilevel Sparse Matrix Problem Algorithms

Algorithms 2 and 3 described in Section 5 rely upon two matrix algebraic
routines for efficiently solving the two-level and three-level versions of the
_multilevel sparse matrix problems_ defined in Nolan & Wand (2020). They
correspond to the SolveTwoLevelSparseMatrix and SolveThreeLevelSparseMatrix
algorithms that we list hereafter as Algorithms S.1 and S.2, respectively.

We briefly describe two-level and three-level sparse matrix structures, and
give explicit definition of the aforementioned routines for efficiently
solving the associated linear system problems, following Appendix A of Nolan
_et al._ (2020).

###  S.2.1 Two-Level Sparse Matrix Problems

Two-level sparse matrix problems are summarized in Section 3.3. Such problems
are efficiently solved by the SolveTwoLevelSparseMatrix routine, which is here
listed as Algorithm S.1 and is justified by Theorem 2.2 of Nolan & Wand
(2020).

  * Inputs: $\boldsymbol{a}_{1}(p\times 1),\boldsymbol{A}_{11}(p\times p),\,\big{\\{}\big{(}\boldsymbol{a}_{2,i}(q\times 1),\boldsymbol{A}_{22,i}(q\times q),\boldsymbol{A}_{12,i}(p\times q)\big{)}:\ 1\leq i\leq m\big{\\}}$

  * $\boldsymbol{\omega}\longleftarrow\boldsymbol{a}_{1}$ ; $\boldsymbol{\Omega}\longleftarrow\boldsymbol{A}_{11}$

  * For $i=1,\ldots,m$:

    * $\boldsymbol{\omega}\longleftarrow\boldsymbol{\omega}-\boldsymbol{A}_{12,i}\boldsymbol{A}_{22,i}^{-1}\boldsymbol{a}_{2,i}$ ; $\boldsymbol{\Omega}\longleftarrow\boldsymbol{\Omega}-\boldsymbol{A}_{12,i}\boldsymbol{A}_{22,i}^{-1}\boldsymbol{A}_{12,i}^{T}$

  * $\boldsymbol{A}^{11}\longleftarrow\boldsymbol{\Omega}^{-1}$ ; $\boldsymbol{x}_{1}\longleftarrow\boldsymbol{A}^{11}\boldsymbol{\omega}$

  * For $i=1,\ldots,m$:

    * $\boldsymbol{x}_{2,i}\longleftarrow\boldsymbol{A}_{22,i}^{-1}(\boldsymbol{a}_{2,i}-\boldsymbol{A}_{12,i}^{T}\boldsymbol{x}_{1})$ ; $\boldsymbol{A}^{12,i}\longleftarrow\,-(\boldsymbol{A}_{22,i}^{-1}\boldsymbol{A}_{12,i}^{T}\boldsymbol{A}^{11})^{T}$

    * $\boldsymbol{A}^{22,i}\longleftarrow\,\boldsymbol{A}_{22,i}^{-1}\big{(}\boldsymbol{I}-\boldsymbol{A}_{12,i}^{T}\boldsymbol{A}^{12,i}\big{)}$

  * Outputs: $\boldsymbol{x}_{1},\boldsymbol{A}^{11},\big{\\{}\big{(}\boldsymbol{x}_{2,i},\boldsymbol{A}^{22,i},\boldsymbol{A}^{12,i}):\ 1\leq i\leq m\big{\\}}$.

Algorithm S.1 The SolveTwoLevelSparseMatrix algorithm for solving the two-
level sparse matrix problem $\boldsymbol{x}=\boldsymbol{A}^{-1}\boldsymbol{a}$
and sub-blocks of $\boldsymbol{A}^{-1}$ corresponding to the non-zero sub-
blocks of $\boldsymbol{A}$.

###  S.2.2 Three-Level Sparse Matrix Problems

Three-level sparse matrix problems are described in Section 3 of Nolan & Wand
(2020). An illustrative three-level sparse matrix example is

$\boldsymbol{A}=\left[\begin{array}[]{cccccccc}\boldsymbol{A}_{11}&\boldsymbol{A}_{12,1}&\boldsymbol{A}_{12,11}&\boldsymbol{A}_{12,12}&\boldsymbol{A}_{12,2}&\boldsymbol{A}_{12,21}&\boldsymbol{A}_{12,22}&\boldsymbol{A}_{12,23}\\\
\boldsymbol{A}_{12,1}^{T}&\boldsymbol{A}_{22,1}&\boldsymbol{A}_{12,1,1}&\boldsymbol{A}_{12,1,2}&\boldsymbol{O}&\boldsymbol{O}&\boldsymbol{O}&\boldsymbol{O}\\\
\boldsymbol{A}_{12,11}^{T}&\boldsymbol{A}_{12,1,1}^{T}&\boldsymbol{A}_{22,11}&\boldsymbol{O}&\boldsymbol{O}&\boldsymbol{O}&\boldsymbol{O}&\boldsymbol{O}\\\
\boldsymbol{A}_{12,12}^{T}&\boldsymbol{A}_{12,1,2}^{T}&\boldsymbol{O}&\boldsymbol{A}_{22,12}&\boldsymbol{O}&\boldsymbol{O}&\boldsymbol{O}&\boldsymbol{O}\\\
\boldsymbol{A}_{12,2}^{T}&\boldsymbol{O}&\boldsymbol{O}&\boldsymbol{O}&\boldsymbol{A}_{22,2}&\boldsymbol{A}_{12,2,1}&\boldsymbol{A}_{12,2,2}&\boldsymbol{A}_{12,2,3}\\\
\boldsymbol{A}_{12,21}^{T}&\boldsymbol{O}&\boldsymbol{O}&\boldsymbol{O}&\boldsymbol{A}_{12,2,1}^{T}&\boldsymbol{A}_{22,21}&\boldsymbol{O}&\boldsymbol{O}\\\
\boldsymbol{A}_{12,22}^{T}&\boldsymbol{O}&\boldsymbol{O}&\boldsymbol{O}&\boldsymbol{A}_{12,2,2}^{T}&\boldsymbol{O}&\boldsymbol{A}_{22,22}&\boldsymbol{O}\\\
\boldsymbol{A}_{12,23}^{T}&\boldsymbol{O}&\boldsymbol{O}&\boldsymbol{O}&\boldsymbol{A}_{12,2,3}^{T}&\boldsymbol{O}&\boldsymbol{O}&\boldsymbol{A}_{22,23}\end{array}\right],$

which corresponds to level 2 group sizes $n_{1}=2$ and $n_{2}=3$, and a level
3 group size $m=2$. A generic three-level sparse matrix $\boldsymbol{A}$
consists of the following components:

  * •

A $p\times p$ matrix $\boldsymbol{A}_{11}$, which is assigned the
$(1,1)$-block position.

  * •

A set of partitioned matrices
$\left\\{\left[\begin{array}[]{c|c|c|c}\boldsymbol{A}_{12,i}&\boldsymbol{A}_{12,ij}&\dots&\boldsymbol{A}_{12,in_{i}}\end{array}\right]:1\leq
i\leq m\right\\}$, which is assigned the $(1,2)$-block position. For each
$1\leq i\leq m$, $\boldsymbol{A}_{12,i}$ is $p\times q_{1}$, and for each
$1\leq j\leq n_{i}$, $\boldsymbol{A}_{12,ij}$ is $p\times q_{2}$.

  * •

A $(2,1)$-block, which is simply the transpose of the $(1,2)$-block.

  * •

A block diagonal structure along the $(2,2)$-block position, where each sub-
block is a two-level sparse matrix, as defined in Section 3.3. For each $1\leq
i\leq m$, $\boldsymbol{A}_{22,i}$ is $q_{1}\times q_{1}$, and for each $1\leq
j\leq n_{i}$, $\boldsymbol{A}_{12,i,j}$ is $q_{1}\times q_{2}$ and
$\boldsymbol{A}_{22,ij}$ is $q_{2}\times q_{2}$.

The three-level sparse matrix problem arising from the illustrative example
above is defined as finding the vector $\boldsymbol{x}$ such that
$\boldsymbol{A}\,\boldsymbol{x}=\boldsymbol{a}$, where:

$\boldsymbol{a}\equiv\left[\begin{array}[]{c}\boldsymbol{a}_{1}\\\
\boldsymbol{a}_{2,1}\\\ \boldsymbol{a}_{2,11}\\\ \boldsymbol{a}_{2,12}\\\
\boldsymbol{a}_{2,2}\\\ \boldsymbol{a}_{2,21}\\\ \boldsymbol{a}_{2,22}\\\
\boldsymbol{a}_{2,23}\end{array}\right]\qquad\text{and}\qquad\boldsymbol{x}\equiv\left[\begin{array}[]{c}\boldsymbol{x}_{1}\\\
\boldsymbol{x}_{2,1}\\\ \boldsymbol{x}_{2,11}\\\ \boldsymbol{x}_{2,12}\\\
\boldsymbol{x}_{2,2}\\\ \boldsymbol{x}_{2,21}\\\ \boldsymbol{x}_{2,22}\\\
\boldsymbol{x}_{2,23}\end{array}\right],$

and determining the sub-blocks of $\boldsymbol{A}^{-1}$ corresponding to the
non-zero sub-blocks of $\boldsymbol{A}$. The structure of
$\boldsymbol{A}^{-1}$ is

$\boldsymbol{A}^{-1}=\left[\begin{array}[]{cccccccc}\boldsymbol{A}^{11}&\boldsymbol{A}^{12,1}&\boldsymbol{A}^{12,11}&\boldsymbol{A}^{12,12}&\boldsymbol{A}^{12,2}&\boldsymbol{A}^{12,21}&\boldsymbol{A}^{12,22}&\boldsymbol{A}^{12,23}\\\
\boldsymbol{A}^{12,1\,T}&\boldsymbol{A}^{22,1}&\boldsymbol{A}^{12,1,1}&\boldsymbol{A}^{12,1,2}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}\\\
\boldsymbol{A}^{12,11\,T}&\boldsymbol{A}^{12,1,1\,T}&\boldsymbol{A}^{22,11}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}\\\
\boldsymbol{A}^{12,12\,T}&\boldsymbol{A}^{12,1,2\,T}&{\LARGE\mbox{$\times$}}&\boldsymbol{A}^{22,12}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}\\\
\boldsymbol{A}^{12,2\,T}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&\boldsymbol{A}^{22,2}&\boldsymbol{A}^{12,2,1}&\boldsymbol{A}^{12,2,2}&\boldsymbol{A}^{12,2,3}\\\
\boldsymbol{A}^{12,21\,T}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&\boldsymbol{A}^{12,2,1\,T}&\boldsymbol{A}^{22,21}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}\\\
\boldsymbol{A}^{12,22\,T}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&\boldsymbol{A}^{12,2,2\,T}&{\LARGE\mbox{$\times$}}&\boldsymbol{A}^{22,22}&{\LARGE\mbox{$\times$}}\\\
\boldsymbol{A}^{12,23\,T}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&\boldsymbol{A}^{12,2,3\,T}&{\LARGE\mbox{$\times$}}&{\LARGE\mbox{$\times$}}&\boldsymbol{A}^{22,23}\end{array}\right].$

For a general three-level sparse linear system problem, $\boldsymbol{a}_{1}$
and $\boldsymbol{x}_{1}$ are $p\times 1$ vectors. For each $1\leq i\leq m$,
$\boldsymbol{a}_{2,i}$ and $\boldsymbol{x}_{2,i}$ are $q_{1}\times 1$ vectors.
For each $1\leq i\leq m$ and $1\leq j\leq n_{i}$ the vectors
$\boldsymbol{a}_{2,ij}$ and $\boldsymbol{x}_{2,ij}$ have dimension
$q_{2}\times 1$.

Such problems are efficiently solved by the SolveThreeLevelSparseMatrix
routine, which is here listed as Algorithm S.2 and is justified by Theorem 3.2
of Nolan & Wand (2020).

  * Input: $\boldsymbol{a}_{1}(p\times 1),\boldsymbol{A}_{11}(p\times p),\,\big{\\{}\big{(}\boldsymbol{a}_{2,i}(q_{1}\times 1),\boldsymbol{A}_{22,i}(q_{1}\times q_{1}),\boldsymbol{A}_{12,i}(p\times q_{1}):\ 1\leq i\leq m\big{\\}}$,
$\qquad\qquad\big{\\{}\boldsymbol{a}_{2,ij}(q_{2}\times
1),\boldsymbol{A}_{22,ij}(q_{2}\times q_{2}),\boldsymbol{A}_{12,ij}(p\times
q_{2}),\boldsymbol{A}_{12,i,j}(q_{1}\times q_{2})\big{)}:\ 1\leq i\leq m,\
1\leq j\leq n_{i}\big{\\}}$.

  * $\boldsymbol{\omega}\longleftarrow\boldsymbol{a}_{1}$ ; $\boldsymbol{\Omega}\longleftarrow\boldsymbol{A}_{11}$

  * For $i=1,\ldots,m$:

    * $\boldsymbol{h}_{2,i}\longleftarrow\boldsymbol{a}_{2,i}$ ; $\boldsymbol{H}_{12,i}\longleftarrow\boldsymbol{A}_{12,i}$ ; $\boldsymbol{H}_{22,i}\longleftarrow\boldsymbol{A}_{22,i}$

    * For $j=1,\ldots,n_{i}$:

      * $\boldsymbol{h}_{2,i}\longleftarrow\boldsymbol{h}_{2,i}-\boldsymbol{A}_{12,i,j}\boldsymbol{A}_{22,ij}^{-1}\boldsymbol{a}_{2,ij}$ ; $\boldsymbol{H}_{12,i}\longleftarrow\boldsymbol{H}_{12,i}-\boldsymbol{A}_{12,ij}\boldsymbol{A}_{22,ij}^{-1}\boldsymbol{A}_{12,i,j}^{T}$
$\boldsymbol{H}_{22,i}\longleftarrow\boldsymbol{H}_{22,i}-\boldsymbol{A}_{12,i,j}\boldsymbol{A}_{22,ij}^{-1}\boldsymbol{A}_{12,i,j}^{T}$

      * $\boldsymbol{\omega}\longleftarrow\boldsymbol{\omega}-\boldsymbol{A}_{12,ij}\boldsymbol{A}_{22,ij}^{-1}\boldsymbol{a}_{2,ij}$ ; $\boldsymbol{\Omega}\longleftarrow\boldsymbol{\Omega}-\boldsymbol{A}_{12,ij}\boldsymbol{A}_{22,ij}^{-1}\boldsymbol{A}_{12,ij}^{T}$

    * $\boldsymbol{\omega}\longleftarrow\boldsymbol{\omega}-\boldsymbol{H}_{12,i}\boldsymbol{H}_{22,i}^{-1}\boldsymbol{h}_{2,i}$ ; $\boldsymbol{\Omega}\longleftarrow\boldsymbol{\Omega}-\boldsymbol{H}_{12,i}\boldsymbol{H}_{22,i}^{-1}\boldsymbol{H}_{12,i}^{T}$

  * $\boldsymbol{A}^{11}\longleftarrow\boldsymbol{\Omega}^{-1}$ ; $\boldsymbol{x}_{1}\longleftarrow\boldsymbol{A}^{11}\boldsymbol{\omega}$

  * For $i=1,\ldots,m$:

    * $\boldsymbol{x}_{2,i}\longleftarrow\,\boldsymbol{H}_{22,i}^{-1}(\boldsymbol{h}_{2,i}-\boldsymbol{H}_{12,i}^{T}\boldsymbol{x}_{1})$ ; $\boldsymbol{A}^{12,i}\longleftarrow\,-(\boldsymbol{H}_{22,i}^{-1}\boldsymbol{H}_{12,i}^{T}\boldsymbol{A}^{11})^{T}$

    * $\boldsymbol{A}^{22,i}\longleftarrow\,\boldsymbol{H}_{22,i}^{-1}(\boldsymbol{I}-\boldsymbol{H}_{12,i}^{T}\boldsymbol{A}^{12,i})$

    * For $j=1,\ldots,n_{i}$:

      * $\boldsymbol{x}_{2,ij}\longleftarrow\boldsymbol{A}_{22,ij}^{-1}\big{(}\boldsymbol{a}_{2,ij}-\boldsymbol{A}_{12,ij}^{T}\boldsymbol{x}_{1}-\boldsymbol{A}_{12,i,j}^{T}\boldsymbol{x}_{2,i}\big{)}$

      * $\boldsymbol{A}^{12,ij}\longleftarrow\,-\big{\\{}\boldsymbol{A}_{22,ij}^{-1}\big{(}\boldsymbol{A}_{12,ij}^{T}\boldsymbol{A}^{11}+\boldsymbol{A}_{12,i,j}^{T}\,\boldsymbol{A}^{12,i\,T}\big{)}\big{\\}}^{T}$

      * $\boldsymbol{A}^{12,i,j}\longleftarrow\,-\big{\\{}\boldsymbol{A}_{22,ij}^{-1}\big{(}\boldsymbol{A}_{12,ij}^{T}\,\boldsymbol{A}^{12,i}+\boldsymbol{A}_{12,i,j}^{T}\,\boldsymbol{A}^{22,i}\big{)}\big{\\}}^{T}$

      * $\boldsymbol{A}^{22,ij}\longleftarrow\boldsymbol{A}_{22,ij}^{-1}\big{(}\boldsymbol{I}-\boldsymbol{A}_{12,ij}^{T}\boldsymbol{A}^{12,ij}-\boldsymbol{A}_{12,i,j}^{T}\,\boldsymbol{A}^{12,i,j}\big{)}$

  * Outputs: $\boldsymbol{x}_{1},\boldsymbol{A}^{11},\big{\\{}\big{(}\boldsymbol{x}_{2,i},\boldsymbol{A}^{22,i},\boldsymbol{A}^{12,i}):\ 1\leq i\leq m\big{\\}},$
$\hbox{}\qquad\qquad\big{\\{}\big{(}\boldsymbol{x}_{2,ij},\boldsymbol{A}^{22,ij},\boldsymbol{A}^{12,ij},\boldsymbol{A}^{12,i,j},\big{)}:\
1\leq i\leq m,\ 1\leq j\leq n_{i}\big{\\}}$.

Algorithm S.2 The SolveThreeLevelSparseMatrix algorithm for solving the three-
level sparse matrix problem $\boldsymbol{x}=\boldsymbol{A}^{-1}\boldsymbol{a}$
and sub-blocks of $\boldsymbol{A}^{-1}$ corresponding to the non-zero sub-
blocks of $\boldsymbol{A}$.

##  Appendix S.3 Derivations

We derive explicit updates for the parameters of the optimal density functions
in (21) using arguments similar to those provided in Neville _et al._ (2014)
with some adjustments. Such updates are then combined with the derivations in
Appendix B of Nolan _et al._ (2020) for deriving the streamlined MFVB
algorithms for our two- and three- level linear mixed-effects models, i.e.
Algorithms 2 and 3.

###  S.3.1 Derivation of $\mathfrak{q}^{*}(\beta_{0},\boldsymbol{\beta})$ and
Associated Parameter Updates

Given the model formulation (22), the full conditional density function for
$(\beta_{0},\boldsymbol{\beta})$ is:

$\displaystyle\mathfrak{p}(\beta_{0},\boldsymbol{\beta}|\text{rest})$
$\displaystyle\propto\mathfrak{p}(\boldsymbol{y}|\beta_{0},\boldsymbol{\beta},\sigma^{2})\>\mathfrak{p}(\beta_{0})\>\mathfrak{p}(\boldsymbol{\beta}|\boldsymbol{\zeta},\tau^{2})$
$\displaystyle\propto\exp\Bigg{\\{}-\frac{1}{2}\Bigg{[}\left[\begin{array}[]{c}\beta_{0}\\\
\boldsymbol{\beta}\end{array}\right]^{T}\left(\frac{1}{\sigma^{2}}\big{[}\>\boldsymbol{1}\>\big{|}\>\boldsymbol{X}\>\big{]}^{T}\big{[}\>\boldsymbol{1}\>\big{|}\>\boldsymbol{X}\>\big{]}+\left[\begin{array}[]{cc}\sigma^{-2}_{\beta_{0}}&\boldsymbol{0}^{T}\\\\[4.30554pt]
\boldsymbol{0}&\tau^{-2}\>\mbox{diag}(\boldsymbol{\zeta})\end{array}\right]\right)\left[\begin{array}[]{c}\beta_{0}\\\
\boldsymbol{\beta}\end{array}\right]$
$\displaystyle\qquad-2\left[\begin{array}[]{c}\beta_{0}\\\
\boldsymbol{\beta}\end{array}\right]^{T}\left(\frac{1}{\sigma^{2}}\big{[}\>\boldsymbol{1}\>\big{|}\>\boldsymbol{X}\>\big{]}^{T}\boldsymbol{y}+\left[\begin{array}[]{cc}\sigma^{-2}_{\beta_{0}}&\boldsymbol{0}^{T}\\\\[4.30554pt]
\boldsymbol{0}&\tau^{-2}\>\mbox{diag}(\boldsymbol{\zeta})\end{array}\right]\left[\begin{array}[]{c}\mu_{\beta_{0}}\\\
\boldsymbol{0}\end{array}\right]\right)\Bigg{]}\Bigg{\\}}$

and application of (8) results in:

$\displaystyle\mathfrak{q}^{*}(\beta_{0},\boldsymbol{\beta})\propto\exp\left\\{E_{\mathfrak{q}}\\{\log\mathfrak{p}(\beta_{0},\boldsymbol{\beta}|\text{rest})\\}\right\\}$
$\displaystyle\propto\exp\Bigg{\\{}-\frac{1}{2}\Bigg{[}\left[\begin{array}[]{c}\beta_{0}\\\
\boldsymbol{\beta}\end{array}\right]^{T}\left(E_{\mathfrak{q}}\\{\sigma^{-2}\\}\big{[}\>\boldsymbol{1}\>\big{|}\>\boldsymbol{X}\>\big{]}^{T}\big{[}\>\boldsymbol{1}\>\big{|}\>\boldsymbol{X}\>\big{]}+\left[\begin{array}[]{cc}\sigma^{-2}_{\beta_{0}}&\boldsymbol{0}^{T}\\\\[4.30554pt]
\boldsymbol{0}&\quad
E_{\mathfrak{q}}\\{\tau^{-2}\\}\>\mbox{diag}(E_{\mathfrak{q}}\\{\boldsymbol{\zeta}\\})\end{array}\right]\right)\left[\begin{array}[]{c}\beta_{0}\\\
\boldsymbol{\beta}\end{array}\right]$
$\displaystyle\qquad-2\left[\begin{array}[]{c}\beta_{0}\\\
\boldsymbol{\beta}\end{array}\right]^{T}\left(E_{\mathfrak{q}}\\{\sigma^{-2}\\}\big{[}\>\boldsymbol{1}\>\big{|}\>\boldsymbol{X}\>\big{]}^{T}\boldsymbol{y}+\left[\begin{array}[]{cc}\sigma^{-2}_{\beta_{0}}&\boldsymbol{0}^{T}\\\\[4.30554pt]
\boldsymbol{0}&E_{\mathfrak{q}}\\{\tau^{-2}\\}\>\mbox{diag}(E_{\mathfrak{q}}\\{\boldsymbol{\zeta}\\})\end{array}\right]\left[\begin{array}[]{c}\mu_{\beta_{0}}\\\
\boldsymbol{0}\end{array}\right]\right)\Bigg{]}\Bigg{\\}}.$

After _completion of the square_ manipulations for the multivariate Gaussian
distribution and standard algebraic manipulations, it follows immediately
that:

$\mathfrak{q}^{*}(\beta_{0},\boldsymbol{\beta})\mbox{ is a
N}\left(\boldsymbol{\mu}_{\mathfrak{q}(\beta_{0},\boldsymbol{\beta})},\boldsymbol{\Sigma}_{\mathfrak{q}(\beta_{0},\boldsymbol{\beta})}\right)\text{
density function}$

with

$\boldsymbol{\Sigma}_{\mathfrak{q}(\beta_{0},\boldsymbol{\beta})}\longleftarrow\left(\mu_{\mathfrak{q}(1/\sigma^{2})}\>\big{[}\>\boldsymbol{1}\>\big{|}\>\boldsymbol{X}\>\big{]}^{T}\big{[}\>\boldsymbol{1}\>\big{|}\>\boldsymbol{X}\>\big{]}+\left[\begin{array}[]{cc}\sigma^{-2}_{\beta_{0}}&\boldsymbol{0}^{T}\\\\[4.30554pt]
\boldsymbol{0}&\mu_{\mathfrak{q}(1/\tau^{2})}\>\mbox{diag}(\boldsymbol{\mu}_{\mathfrak{q}(\boldsymbol{\zeta})})\end{array}\right]\right)^{-1}$

and

$\boldsymbol{\mu}_{\mathfrak{q}(\beta_{0},\boldsymbol{\beta})}\longleftarrow\boldsymbol{\Sigma}_{\mathfrak{q}(\beta_{0},\boldsymbol{\beta})}\left(\mu_{\mathfrak{q}(1/\sigma^{2})}\>\big{[}\>\boldsymbol{1}\>\big{|}\>\boldsymbol{X}\>\big{]}^{T}\boldsymbol{y}+\left[\begin{array}[]{c}\mu_{\beta_{0}}/\sigma^{2}_{\beta_{0}}\\\
\boldsymbol{0}\end{array}\right]\right).$

###  S.3.2 Derivation of $\mathfrak{q}^{*}(\sigma^{2})$ and Associated
Parameter Updates

Given the model formulation (22), the full conditional density function for
$\sigma^{2}$ is:

$\displaystyle\mathfrak{p}(\sigma^{2}|\text{rest})$
$\displaystyle\propto\mathfrak{p}(\boldsymbol{y}|\beta_{0},\boldsymbol{\beta},\sigma^{2})\>\mathfrak{p}(\sigma^{2}|a_{\sigma^{2}})$
$\displaystyle\propto(\sigma^{2})^{-(\nu_{\sigma^{2}}+n)/2-1}\exp\left\\{-\frac{1}{\sigma^{2}}\frac{a_{\sigma^{2}}^{-1}+\lVert\boldsymbol{y}-\boldsymbol{1}\beta_{0}-\boldsymbol{X}\boldsymbol{\beta}\rVert^{2}}{2}\right\\}$

and application of (8) results in:

$\displaystyle\mathfrak{q}^{*}(\sigma^{2})$
$\displaystyle\propto\exp\left\\{E_{\mathfrak{q}}\\{\log\mathfrak{p}(\sigma^{2}|\text{rest})\\}\right\\}$
$\displaystyle\propto(\sigma^{2})^{-(\nu_{\sigma^{2}}+n)/2-1}\exp\left\\{-\frac{1}{\sigma^{2}}\frac{E_{\mathfrak{q}}\\{a_{\sigma^{2}}^{-1}\\}+E_{\mathfrak{q}}\\{\lVert\boldsymbol{y}-\boldsymbol{1}\beta_{0}-\boldsymbol{X}\boldsymbol{\beta}\rVert^{2}\\}}{2}\right\\}.$

After standard algebraic manipulations, it follows immediately that:

$\mathfrak{q}^{*}(\sigma^{2})\mbox{ is an
Inverse-}\chi^{2}\left(\xi_{\mathfrak{q}(\sigma^{2})},\lambda_{\mathfrak{q}(\sigma^{2})}\right)\text{
density function}$

with

$\xi_{\mathfrak{q}(\sigma^{2})}\longleftarrow n+\nu_{\sigma^{2}}$

and

$\lambda_{\mathfrak{q}(\sigma^{2})}\longleftarrow\mu_{\mathfrak{q}(1/a_{\sigma^{2}})}+\left\lVert\boldsymbol{y}-\big{[}\>\boldsymbol{1}\>\big{|}\>\boldsymbol{X}\>\big{]}\boldsymbol{\mu}_{\mathfrak{q}(\beta_{0},\boldsymbol{\beta})}\right\rVert^{2}+\mbox{tr}{\left\\{\boldsymbol{\Sigma}_{\mathfrak{q}(\beta_{0},\boldsymbol{\beta})}\big{[}\>\boldsymbol{1}\>\big{|}\>\boldsymbol{X}\>\big{]}^{T}\big{[}\>\boldsymbol{1}\>\big{|}\>\boldsymbol{X}\>\big{]}\right\\}}.$

###  S.3.3 Derivation of $\mathfrak{q}^{*}(a_{\sigma^{2}})$ and Associated
Parameter Updates

Given the model formulation (22), the full conditional density function for
$a_{\sigma^{2}}$ is:

$\displaystyle\mathfrak{p}(a_{\sigma^{2}}|\text{rest})$
$\displaystyle\propto\mathfrak{p}(\sigma^{2}|a_{\sigma^{2}})\>\mathfrak{p}(a_{\sigma^{2}})$
$\displaystyle\propto(a_{\sigma^{2}})^{-(\nu_{\sigma^{2}}+1)/2-1}\exp\left\\{-\frac{1}{a_{\sigma^{2}}}\frac{\sigma^{-2}+(\nu_{\sigma^{2}}s_{\sigma^{2}}^{2})^{-1}}{2}\right\\}$

and application of (8) results in:

$\displaystyle\mathfrak{q}^{*}(a_{\sigma^{2}})$
$\displaystyle\propto\exp\left\\{E_{\mathfrak{q}}\\{\log\mathfrak{p}(a_{\sigma^{2}}|\text{rest})\\}\right\\}$
$\displaystyle\propto(a_{\sigma^{2}})^{-(\nu_{\sigma^{2}}+1)/2-1}\exp\left\\{-\frac{1}{a_{\sigma^{2}}}\frac{E_{\mathfrak{q}}\\{\sigma^{-2}\\}+(\nu_{\sigma^{2}}s_{\sigma^{2}}^{2})^{-1}}{2}\right\\}.$

After standard algebraic manipulations, it follows immediately that:

$\mathfrak{q}^{*}(a_{\sigma^{2}})\mbox{ is an
Inverse-}\chi^{2}\left(\xi_{\mathfrak{q}(a_{\sigma^{2}})},\lambda_{\mathfrak{q}(a_{\sigma^{2}})}\right)\text{
density function}$

with

$\xi_{\mathfrak{q}(a_{\sigma^{2}})}\longleftarrow
1+\nu_{\sigma^{2}}\quad\text{and}\quad\lambda_{\mathfrak{q}(a_{\sigma^{2}})}\longleftarrow\mu_{\mathfrak{q}(1/\sigma^{2})}+(\nu_{\sigma^{2}}s_{\sigma^{2}}^{2})^{-1}.$

###  S.3.4 Derivation of $\mathfrak{q}^{*}(\tau^{2})$ and Associated Parameter
Updates

Given the model formulation (22), the full conditional density function for
$\tau^{2}$ is:

$\displaystyle\mathfrak{p}(\tau^{2}|\text{rest})$
$\displaystyle\propto\left\\{\prod_{h=1}^{H}\mathfrak{p}(\beta_{h}|\tau^{2})\right\\}\mathfrak{p}(\tau^{2}|a_{\tau^{2}})$
$\displaystyle\propto(\tau^{2})^{-(H+1)/2-1}\exp\left\\{-\frac{1}{\tau^{2}}\frac{a_{\tau^{2}}^{-1}+\sum_{h=1}^{H}\zeta_{h}(\beta_{h})^{2}}{2}\right\\}$

and application of (8) results in:

$\displaystyle\mathfrak{q}^{*}(\tau^{2})$
$\displaystyle\propto\exp\left\\{E_{\mathfrak{q}}\\{\log\mathfrak{p}(\tau^{2}|\text{rest})\\}\right\\}$
$\displaystyle\propto(\tau^{2})^{-(H+1)/2-1}\exp\left\\{-\frac{1}{\tau^{2}}\frac{E_{\mathfrak{q}}\\{a_{\tau^{2}}^{-1}\\}+\sum_{h=1}^{H}E_{\mathfrak{q}}\\{\zeta_{h}\\}E_{\mathfrak{q}}\\{\beta_{h}^{2}\\}}{2}\right\\}.$

After standard algebraic manipulations, it follows immediately that:

$\mathfrak{q}^{*}(\tau^{2})\mbox{ is an
Inverse-}\chi^{2}\left(\xi_{\mathfrak{q}(\tau^{2})},\lambda_{\mathfrak{q}(\tau^{2})}\right)\text{
density function}$

with

$\xi_{\mathfrak{q}(\tau^{2})}\longleftarrow
H+1\quad\text{and}\quad\lambda_{\mathfrak{q}(\tau^{2})}\longleftarrow\mu_{\mathfrak{q}(1/a_{\tau^{2}})}+\boldsymbol{\mu}_{\mathfrak{q}(\boldsymbol{\zeta})}^{T}\boldsymbol{\mu}_{\mathfrak{q}((\boldsymbol{\beta})^{2})}.$

###  S.3.5 Derivation of $\mathfrak{q}^{*}(a_{\tau^{2}})$ and Associated
Parameter Updates

Given the model formulation (22), the full conditional density function for
$a_{\tau^{2}}$ is:

$\displaystyle\mathfrak{p}(a_{\tau^{2}}|\text{rest})$
$\displaystyle\propto\mathfrak{p}(\tau^{2}|a_{\tau^{2}})\mathfrak{p}(a_{\tau^{2}})$
$\displaystyle\propto(a_{\tau^{2}})^{-2}\exp\left\\{-\frac{1}{a_{\tau^{2}}}\frac{\tau^{-2}+s_{\tau^{2}}^{-2}}{2}\right\\}$

and application of (8) results in:

$\displaystyle\mathfrak{q}^{*}(a_{\tau^{2}})\propto\exp\left\\{E_{\mathfrak{q}}\\{\log\mathfrak{p}(a_{\tau^{2}}|\text{rest})\\}\right\\}\propto(a_{\tau^{2}})^{-2}\exp\left\\{-\frac{1}{a_{\tau^{2}}}\frac{E_{\mathfrak{q}}\\{\tau^{-2}\\}+s_{\tau^{2}}^{-2}}{2}\right\\}.$

After standard algebraic manipulations, it follows immediately that:

$\mathfrak{q}^{*}(a_{\tau^{2}})\mbox{ is an
Inverse-}\chi^{2}\left(\xi_{\mathfrak{q}(a_{\tau^{2}})},\lambda_{\mathfrak{q}(a_{\tau^{2}})}\right)\text{
density function}$

with

$\xi_{\mathfrak{q}(a_{\tau^{2}})}\longleftarrow
2\quad\text{and}\quad\lambda_{\mathfrak{q}(a_{\tau^{2}})}\longleftarrow\mu_{\mathfrak{q}(1/\tau^{2})}+s^{-2}_{\tau^{2}}.$

###  S.3.6 Derivation of $\mathfrak{q}^{*}(\zeta_{h})$ and Associated
Parameter Updates

Given the model formulation (22), the full conditional density function for
$\zeta_{h}$, $1\leq h\leq H$, is:

$\mathfrak{p}(\zeta_{h}|\text{rest})\propto\mathfrak{p}(\beta_{h}|\zeta_{h},\tau)\>\mathfrak{p}(\zeta_{h}|a_{\zeta_{h}}),$

where $\mathfrak{p}(\beta_{h}|\zeta_{h})$ is the density function of a
$\mbox{N}(0,\tau^{2}/\zeta_{h})$ distribution, while
$\mathfrak{p}(\zeta_{h}|a_{\zeta_{h}})$ depends on the hierarchical
specification given in Table 1 for each of the considered global-local priors.

####  S.3.6.1 Laplace Prior Specification

If $\beta_{h}|\tau\sim\mbox{Laplace}(0,\tau)$ then:

$\displaystyle\mathfrak{p}(\zeta_{h}|\text{rest})\propto(\zeta_{h})^{-3/2}\exp\left\\{-\zeta_{h}\frac{\beta_{h}^{2}}{2\tau^{2}}-\frac{1}{2\zeta_{h}}\right\\}$

and application of (8) results in:

$\displaystyle\mathfrak{q}^{*}(\zeta_{h})\propto\exp\left\\{E_{\mathfrak{q}}\\{\log\mathfrak{p}(\zeta_{h}|\text{rest})\\}\right\\}\propto(\zeta_{h})^{-3/2}\exp\left\\{-\zeta_{h}\frac{E_{\mathfrak{q}}\\{\beta_{h}^{2}\\}E_{\mathfrak{q}}\\{\tau^{-2}\\}}{2}-\frac{1}{2\zeta_{h}}\right\\}.$

After standard algebraic manipulations, it follows immediately that:

$\mathfrak{q}^{*}(\zeta_{h})\mbox{ is an Inverse-
Gaussian}(\mu_{\mathfrak{q}(\zeta_{h})},1)\mbox{ density function }$

with

$\mu_{\mathfrak{q}(\zeta_{h})}\longleftarrow\sqrt{1\Big{/}\left(\mu_{\mathfrak{q}(1/\tau^{2})}\mu_{\mathfrak{q}(\beta_{h}^{2})}\right)}.$

####  S.3.6.2 Horseshoe Prior Specification

If $\beta_{h}^{|}\tau\sim\mbox{Horseshoe}(0,\tau)$ then:

$\displaystyle\mathfrak{p}(\zeta_{h}|\text{rest})\propto\exp\left\\{-\zeta_{h}\left(\frac{\beta_{h}^{2}}{2\tau^{2}}+a_{\zeta_{h}}\right)\right\\}$

and application of (8) results in:

$\displaystyle\mathfrak{q}^{*}(\zeta_{h})\propto\exp\\{E_{\mathfrak{q}}\\{\log\mathfrak{p}(\zeta_{h}|\text{rest})\\}\\}\propto\exp\left\\{-\zeta_{h}\left(\frac{E_{\mathfrak{q}}\\{\beta_{h}^{2}\\}E_{\mathfrak{q}}\\{\tau^{-2}\\}}{2}+E_{\mathfrak{q}}\\{a_{\zeta_{h}}\\}\right)\right\\}.$

After standard algebraic manipulations, it follows immediately that:

$\mathfrak{q}^{*}(\zeta_{h})\mbox{ is a
Gamma}(1,\lambda_{\mathfrak{q}(\zeta_{h})})\mbox{ density function }$

with

$\lambda_{\mathfrak{q}(\zeta_{h})}\longleftarrow\frac{\mu_{\mathfrak{q}(1/\tau^{2})}\mu_{\mathfrak{q}(\beta_{h}^{2})}}{2}+\mu_{\mathfrak{q}(a_{\zeta_{h}})}.$

####  S.3.6.3 Normal-Exponential-Gamma Prior Specification

If $\beta_{h}|\tau\sim\mbox{NEG}(0,\tau,\lambda)$ then:

$\displaystyle\mathfrak{p}(\zeta_{h}|\text{rest})\propto(\zeta_{h})^{-3/2}\exp\left\\{-\zeta_{h}\frac{\beta_{h}^{2}}{2\tau^{2}}-\frac{a_{\zeta_{h}}}{\zeta_{h}}\right\\}$

and application of (8) results in:

$\displaystyle\mathfrak{q}^{*}(\zeta_{h})\propto\exp\left\\{E_{\mathfrak{q}}\\{\log\mathfrak{p}(\zeta_{h}|\text{rest})\\}\right\\}\propto(\zeta_{h})^{-3/2}\exp\left\\{-\zeta_{h}\frac{E_{\mathfrak{q}}\\{\beta_{h}^{2}\\}E_{\mathfrak{q}}\\{\tau^{-2}\\}}{2}-\frac{E_{\mathfrak{q}}\\{a_{\zeta_{h}}\\}}{\zeta_{h}}\right\\}.$

After standard algebraic manipulations, it follows immediately that:

$\mathfrak{q}^{*}(\zeta_{h})\mbox{ is an Inverse-
Gaussian}(\mu_{\mathfrak{q}(\zeta_{h})},\lambda_{\mathfrak{q}(\zeta_{h})})\mbox{
density function }$

with

$\mu_{\mathfrak{q}(\zeta_{h})}\longleftarrow\sqrt{2\mu_{\mathfrak{q}(a_{\zeta_{h}})}\Big{/}\left(\mu_{\mathfrak{q}(1/\tau^{2})}\mu_{\mathfrak{q}(\beta_{h}^{2})}\right)}\quad\text{and}\quad\lambda_{\mathfrak{q}(\zeta_{h})}\longleftarrow
2\mu_{\mathfrak{q}(a_{\zeta_{h}})}.$

###  S.3.7 Derivation of $\mathfrak{q}^{*}(a_{\zeta_{h}})$ and Associated
Parameter Updates

Given the model formulation (22), the full conditional density function for
$\zeta_{h}$, $1\leq h\leq H$, is:

$\mathfrak{p}(a_{\zeta_{h}}|\text{rest})\propto\mathfrak{p}(\zeta_{h}|a_{\zeta_{h}})\>\mathfrak{p}(a_{\zeta_{h}}),$

where both $\mathfrak{p}(\zeta_{h}|a_{\zeta_{h}})$ and
$\mathfrak{p}(a_{\zeta_{h}})$ depend on the hierarchical specification given
in Table 1 for each of the considered global-local priors.

####  S.3.7.1 Laplace Prior Specification

According to Table 1, for the Laplace prior it is not necessary to introduce
the $a_{\zeta_{h}}$ auxiliary variables. Hence $\mathfrak{q}(a_{\zeta_{h}})$
is not present in (24) and does not need to be determined.

####  S.3.7.2 Horseshoe Prior Specification

If $\beta_{h}|\tau\sim\mbox{Horseshoe}(0,\tau)$ then:

$\displaystyle\mathfrak{p}(a_{\zeta_{h}}|\text{rest})\propto\exp\\{-(\zeta_{h}+1)a_{\zeta_{h}}\\}$

and application of (8) results in:

$\displaystyle\mathfrak{q}^{*}(a_{\zeta_{h}})\propto\exp\\{E_{\mathfrak{q}}\\{\log\mathfrak{p}(a_{\zeta_{h}}|\text{rest})\\}\\}\propto\exp\left\\{-a_{\zeta_{h}}\left(E_{\mathfrak{q}}\\{\zeta_{h}\\}+1\right)\right\\}.$

After standard algebraic manipulations, it follows immediately that:

$\mathfrak{q}^{*}(a_{\zeta_{h}})\mbox{ is a
Gamma}(1,\lambda_{\mathfrak{q}(a_{\zeta_{h}})})\mbox{ density function }$

with

$\lambda_{\mathfrak{q}(\zeta_{h})}\longleftarrow\mu_{\mathfrak{q}(\zeta_{h})}+1.$

####  S.3.7.3 Normal-Exponential-Gamma Prior Specification

If $\beta_{h}|\tau\sim\mbox{NEG}(0,\tau,\lambda)$ then:

$\displaystyle\mathfrak{p}(a_{\zeta_{h}}|\text{rest})\propto(a_{\zeta_{h}})^{(\lambda+1)-1}\exp\\{a_{\zeta_{h}}(\zeta_{h}^{-1}+1)\\}$

and application of (8) results in:

$\displaystyle\mathfrak{q}^{*}(a_{\zeta_{h}})\propto\exp\left\\{E_{\mathfrak{q}}\\{\log\mathfrak{p}(a_{\zeta_{h}}|\text{rest})\\}\right\\}\propto(a_{\zeta_{h}})^{(\lambda+1)-1}\exp\\{a_{\zeta_{h}}(E_{\mathfrak{q}}\\{\zeta_{h}^{-1}\\}+1)\\}.$

After standard algebraic manipulations, it follows immediately that:

$\mathfrak{q}^{*}(\zeta_{h})\mbox{ is a
Gamma}(\lambda+1,\lambda_{\mathfrak{q}(a_{\zeta_{h}})})\mbox{ density function
}$

with

$\lambda_{\mathfrak{q}(\zeta_{h})}\longleftarrow\mu_{\mathfrak{q}(1/\zeta_{h})}+1.$