# Tumor Radiogenomics with Bayesian Layered Variable Selection

Shariq Mohammeda,b,∗, Sebastian Kurtekc, Karthik Bharathd,
Arvind Raoa,b and Veerabhadran Baladandayuthapania

aDepartment of Biostatistics, University of Michigan
bDepartment of Computational Medicine & Bioinformatics, University of Michigan
cDepartment of Statistics, Ohio State University
dSchool of Mathematical Sciences, University of Nottingham
∗corresponding author. E-mail<EMAIL_ADDRESS>
###### Abstract

We propose a statistical framework to integrate radiological magnetic
resonance imaging (MRI) and genomic data to identify the underlying
radiogenomic associations in lower grade gliomas (LGG). We devise a novel
imaging phenotype by dividing the tumor region into concentric spherical
layers that mimics the tumor evolution process. MRI data within each layer is
represented by voxel–intensity-based probability density functions which
capture the complete information about tumor heterogeneity. Under a
Riemannian-geometric framework these densities are mapped to a vector of
principal component scores which act as imaging phenotypes. Subsequently, we
build Bayesian variable selection models for each layer with the imaging
phenotypes as the response and the genomic markers as predictors. Our novel
hierarchical prior formulation incorporates the interior-to-exterior structure
of the layers, and the correlation between the genomic markers. We employ a
computationally-efficient Expectation–Maximization-based strategy for
estimation. Simulation studies demonstrate the superior performance of our
approach compared to other approaches. With a focus on the cancer driver genes
in LGG, we discuss some biologically relevant findings. Genes implicated with
survival and oncogenesis are identified as being associated with the spherical
layers, which could potentially serve as early-stage diagnostic markers for
disease monitoring, prior to routine invasive approaches.

Keywords: cancer driver genes, lower grade gliomas, radiogenomic associations,
spike-and-slab prior.

##  1 Introduction

Low-grade gliomas (LGG) are infiltrative brain neoplasms characterized as
World Health Organization (WHO) grade II and III neoplasms. LGG is a uniformly
fatal disease of young adults (mean age: 41 years), with survival times
averaging approximately seven years (Claus et al.,, 2015). LGG patients
usually have better survival than patients with high-grade (WHO grade IV)
gliomas. However, some of the LGG tumors recur after treatment and progress
into high-grade tumors, making it essential to understand the underlying
etiology, and to improve the treatment management and monitoring for LGG
patients.

Due to recent technological advances, large multi-modal datasets are being
produced; these include imaging as well as multi-platform genomics data. In
complex disease systems such as cancer, integrative analyses of such data can
reveal important biological insights into specific disease mechanisms and its
subsequent clinical translation (Wang et al.,, 2013; Morris and
Baladandayuthapani,, 2017). Genomic profiling technologies, such as
microarrays, next-generation sequencing, methylation arrays and proteomic
analyses, have facilitated thorough investigations at the molecular level.
Such multi-platform genomic data resources have been used to develop models to
better understand the molecular characterization across cancers and especially
in gliomas (Ohgaki et al.,, 2004; Ohgaki and Kleihues,, 2007). While genomic
data provides information on the molecular characterization of the disease,
radiological imaging, such as magnetic resonance imaging (MRI), computed
tomography (CT) and positron emission tomography (PET), provide complementary
information about the structural aspects of the disease. In this context,
_radiomic_ analysis involves mining and extraction of various types of
quantitative imaging features from different modalities obtained through high-
throughput radiology images (Bakas et al., 2017b, ). These image-derived
features describe various characteristics such as morphology and texture,
among others. Integrative analyses of genomic and radiomic features, commonly
referred to as radiogenomic analysis, capture complementary characteristics of
the underlying tumor (Gevaert et al.,, 2014; Mazurowski et al.,, 2017).
However, such integrative modeling approaches present multiple analytical and
computational challenges including incorporating complex biological structure,
and high-dimensionality of quantitative imaging/genomic markers, necessitating
principled and biologically-informed dimension reduction and information
extraction techniques.

Existing radiomic and radiogenomic studies typically consider voxel-level
summary features (based on intensity histogram, geometry, and texture
analyses) as representations of a tumor image, to assess the progression (or
regression) of tumors using statistical or machine learning models (Lu et
al.,, 2018; Elshafeey et al.,, 2019). The main drawbacks of using voxel-
intensity summary statistics/features are: (i) they do not capture the entire
information in the voxel intensity histogram; (ii) there is subjectivity in
the choice of the number of features and the location where they are computed.
Consequently, the radiomic/radiogenomic analyses using summary features are
unable to detect potential small-scale and sensitive changes in the tumor
(Just,, 2014). To address these challenges, we propose a de novo strategy for
quantifying the tumor heterogeneity using voxel-intensity based imaging
phenotypes that mimic the tumor evolution process.

###  1.1 De Novo Quantification of Imaging Phenotypes

Tumors typically evolve from a single (or group of) cancerous cell(s) and grow
multiplicatively via cross-talk with other nearby cells, thereby producing a
central necrosis region (dead cells) (Swanson et al.,, 2003). Hence, as we
move from the inner core of a tumor toward the exterior, different layers of
the tumor possess distinct tissue characteristics. Therefore, it is important
to construct imaging phenotypes which (a) mimic the underlying tumor evolution
process, and (b) efficiently harness the tumor heterogeneity from the imaging
scans, for downstream modeling and interpretations. To address (a), we divide
the tumor region into one inner sphere followed by concentric spherical shells
which potentially capture the tumor growth process. In Figure 1, we show a
visualization of an example tumor region divided into three spherical layers.
Tumor spheroids (cell cultures with necrotic core and peripheral layer of
proliferative cells) have been used as reliable models of in vitro solid
tumors (Breslin and O’Driscoll,, 2013; Costa et al.,, 2016). To address (b),
we propose to work with the probability density functions (smoothed
histograms) constructed using layer-wise tumor voxel intensities from three-
dimensional (3D) MRI scans. These probability density functions (PDFs) are the
data objects of interest and capture the entire information about the
distribution of voxel intensities, no longer requiring any summaries of the
distribution. Recent studies have successfully utilized density-based
approaches in unsupervised clustering (Saha et al.,, 2016) and regression
models for radiogenomic analyses (Mohammed et al.,, 2021; Yang et al.,, 2020).

---
Figure 1: A graphic to visualize an example tumor region divided into three
spherical layers. The object on the left shows three spherical layers colored
in shades of red, blue and green. For each of these layers, the curves on the
right depict corresponding PDFs.

Briefly, we propose a novel approach for integrating multimodal data to
identify the genomic markers that have significant associations with the PDF-
based radiomic phenotypes. The imaging phenotype we consider arises from MRI,
as it provides high-resolution images with a wide range of contrasts through
multiple imaging sequences (see Section 2). The radiomic phenotypes are
constructed by considering the aforementioned topological structure of the
tumor and mimicking its evolution process. Our statistical framework
identifies significant genomic associations on the different concentric
spherical shells. These associations better inform the complex interplay
between the molecular signatures and imaging characteristics in LGGs, and can
provide effective non-invasive diagnostic options by monitoring the associated
genomic markers prior to invasive alternatives such as biopsy. A schematic
representation of our approach is presented in Figure 2 (see Section 2 for
details). This figure depicts an end-to-end workflow, which includes (i)
constructing spherical layers from the tumor voxels in MRI scans and computing
the corresponding imaging phenotypes, and (ii) building a statistical
framework which includes the structural information from the data to identify
radiogenomic associations, elements of which are summarized below.

Figure 2: Schematic representation of the proposed modeling approach. The
tumor voxels are divided into spherical layers across all of the MRI
sequences. Under a Riemannian-geometric framework, we construct principal
component (PC) scores for each spherical layer across all four MRI sequences.
With these PC scores as the multivariate response, and the genomic markers as
predictors, we build a layer-wise Bayesian variable selection model. The
estimation occurs sequentially, allowing for borrowing of information between
layers through prior specification.

Constructing one global PDF for each tumor disregards structural information
about tissue characteristics of the tumor. This would prohibit us from
analyzing any potential associations of genomic markers with local (regional)
aspects of the tumor. We propose to evaluate associations between genomic
signatures of the tumor through large (or subtle) changes in the overall shape
of the voxel-intensity PDFs from biologically motivated tumor
compartments/layers. We address this by constructing $\tau$ PDFs corresponding
to $\tau$ tumor layers. Tumors for different subjects will have distinct
characteristics of these layers, i.e., not all tumors might have structural
information about all of the layers. However, tumor spheroids have been used
as in vitro models of solid tumors (Breslin and O’Driscoll,, 2013; Costa et
al.,, 2016). Hence, we divide the tumor region into $\tau$ spherical shells to
mimic the tumor growth process. The number of layers the tumor is divided into
remains the same across subjects. We effectively have $\tau$ voxel-intensity
PDFs as $x\mapsto\big{(}f_{1}(x),\ldots,f_{\tau}(x)\big{)}$, a multivariate
functional data object with each component as a PDF. The space of PDFs is a
non-linear, infinite-dimensional manifold, which poses significant challenges
in their analysis. We use the Riemannian-geometric framework based on the
Fisher-Rao metric on the space of PDFs (Section 2) that (i) leads to a
tractable representation of each PDF, and (ii) engenders dimension reduction
of the functional objects, in a manner compatible with the requirements of
downstream statistical analyses. This process yields a set of imaging
phenotypes based on PDFs corresponding to each spherical layer in each MRI
sequence. The natural dependence between the component functions is
incorporated through a novel Bayesian prior formulation (Section 3).

###  1.2 Statistical Framework

We build sequential regression models for each spherical layer starting from
the inner-most sphere and moving outwards towards the edge of the tumor. The
sets of imaging phenotypes from a spherical layer across multiple MRI
sequences are then used as the multivariate response in the regression
corresponding to that spherical layer, with the genomic markers as the
corresponding predictors. This specific structured form of the imaging
phenotypes as well as assessing their interaction with genomic covariates
leads to the following challenges:

  1. 1.

While handling high-resolution imaging data, the number of parameters to
estimate is usually much higher than the number of subjects, leading to the
curse of dimensionality ($p\gg n$).

  2. 2.

Splitting the tumor voxels into concentric spherical layers induces a natural
sequential ordering imitating tumor growth, which needs to be incorporated
into the model formulation. Additionally, correlation between genomic markers
needs to be accounted for as well.

  3. 3.

High-dimensionality and sequential structure in imaging data, and correlation
structure in genomic markers, create severe computational hurdles warranting
scalable approaches.

In this paper, we build a novel statistical framework in a multiple-
multivariate regression set-up, which addresses the above mentioned challenges
incorporating the aforementioned dependence structures within and between the
responses and covariates. Specifically, we leverage the hierarchical model
specification to incorporate the sequential nature of the spherical layer-
based imaging phenotypes as well as the correlation between the genomic
markers. We propose a Bayesian variable selection approach that allows us to
incorporate prior structural information in a natural and principled manner
within a decision-theoretic framework, and deals with the issue of high-
dimensionality. Specifically, we use continuous structured spike-and-slab
priors (George and McCulloch,, 1997; Ishwaran et al.,, 2005; Andersen et al.,,
2014) within the regression model for a specific layer to (i) blend prior
information about the correlation between the genomic markers, and (ii) borrow
information about the selected genomic markers from the estimation at one
spherical layer and propagate the information through the prior of the model
to the subsequent spherical layer(s). As opposed to existing approaches that
incorporate apriori structural information (Li and Zhang,, 2010; Vannucci and
Stingo,, 2010), our prior formulation incorporates structural information not
only for the dependence between the predictors, but also for the dependence
between the multivariate response and across the spherical layers.

Finally, the complex structure of the model warrants scalable estimation
approaches. Toward this end, as an alternative to computationally expensive
Markov chain Monte Carlo (MCMC) sampling techniques, we employ an
Expectation–Maximization (EM)-based estimation strategy (Ročková and George,,
2014), where we only search for the posterior mode to iteratively estimate the
posterior selection probabilities. Our EM-based estimation is an
generalization of Ročková and George, (2014) to a multiple-multivariate
regression setting, and provides a scalable alternative to sampling
strategies; it also facilitates model selection. While developed in the
context of this specific imaging-genomics case study, our sequential modeling
framework is general, has independent methodological value and could also be
used in other applications, which exhibit a natural ordering, e.g.
multivariate data collected over time, space, or other axes.

The rest of this paper is organized as follows. In Section 2, we describe the
radiological imaging data and the construction of the imaging phenotypes. We
present our statistical modeling framework in Section 3, including model
specification, details about the EM-based estimation (Section 4), and model
selection (Section 4.2). In Section 5, we include results from simulation
studies under different scenarios to assess the performance of the proposed
modeling approach. In Section 6, we provide details of the genomic data,
present results of our radiogenomic analysis of LGGs, and shed light upon some
important biological findings (Section 6). Finally, in Section 7, we conclude
with a discussion and some directions for future work. The indices to
sections, figures, tables, and equations in the supplementary material are
preceded by ‘S’ (e.g. Section S1).

##  2 Data Characteristics and Imaging Phenotypes

In this section, we describe the MRI data characteristics and present the
explicit quantification of density-based imaging phenotypes. MRI scans are a
rich source of imaging data as they provide a wide range of imaging contrasts.
Consequently, different tissue characteristics in the brain are highlighted
differently by the various MRI sequences. Primary MRI sequences include (i)
native (T1), (ii) post-contrast T1-weighted (T1Gd), (iii) T2-weighted (T2),
and (iv) T2 fluid attenuated inversion recovery (FLAIR). The utility of these
imaging sequences in providing complementary information can be seen in Figure
S1, where the tumor sub-regions appear with varying contrasts. MRI scans
typically include data on all four imaging sequences and segmentation labels
that generate a mask indicating the tumor and non-tumor regions in the MRI
scan (top-left panel in Figure 2). We specifically consider MRI scans of 65
LGG subjects obtained from The Cancer Imaging Archive (TCIA) (Clark et al.,,
2013). Further details about the data are provided in Section 6. Each MRI scan
is structured as a 3D array where the third axis corresponds to the different
axial slices of the image. Hence, we have four 3D arrays, each corresponding
to one of the four imaging sequences, and an additional 3D array for the
accompanying segmentation mask. All four MRI scans, and the segmentation mask,
have a voxel-to-voxel correspondence which allows for clear specification of
the tumor region. In Figure S1, we show an example of the axial slice from a
brain MRI for a subject with LGG for each of the four imaging sequences. The
region inside the red boundary overlaid on these images indicates the
segmented tumor region. To integrate MRI data into our modeling, we first
construct imaging phenotypes (or radiomic phenotypes) based on the voxel
intensities of the MRI scans.

Construction of PDF-based imaging phenotypes. Consider MRI scans for $n$ LGG
subjects across all four imaging sequences and the corresponding tumor
segmentation masks. Let $N_{i}$ denote the number of voxels in the tumor
region for each subject $i$. Let $\mbox{\boldmath{$\mathbf{c}$}}_{i}$ denote
array indices of the mid-point for the minimal bounding box of the tumor
region, and $\tau\in\mathbb{N}$ denote the number of spherical shells we
divide the tumor region into. In Algorithm S1 of Section S1, we describe
construction of kernel density estimates $f^{(t,m)}_{i}$ corresponding to
spherical shell $t\in\\{1,\ldots,\tau\\}$ from MRI sequence $m\in
M=\\{T1,T1Gd,T2,FLAIR\\}$ for subject $i$. We divide the tumor region into
$\tau$ spherical shells based on the total number of voxels in the tumor
region, so that we have the same number of spherical shells for all subjects.
While the innermost layer is a sphere, we refer to all of the $\tau$ layers as
spherical shells or layers hereafter. Hence, for each subject, we have the
tumor region divided into $\tau$ spherical shells across all four imaging
sequences. This is represented in panel A of Figure 2, with spherical layers
from four MRI sequences and their corresponding PDFs. Via Algorithm S1, we
obtain kernel density estimates $f^{(t,m)}_{i}$ based on the voxels from MRI
sequence $m$ in spherical shell $D^{(t)}_{i}$ for subject $i$. Using a
Riemannian-geometric framework for analyzing PDFs (Srivastava et al.,, 2007),
we implement a principal component analysis (PCA) at each layer
$t=1,\ldots,\tau$, which generates layer-specific PC scores (as shown in
Figure 2). The geometric framework for PDFs is necessary in order to ensure
that the PCA is compatible with the nonlinear structure of the space of the
voxel-intensity PDFs, and has been successfully used in several works (Saha et
al.,, 2016; Mohammed et al.,, 2021). For ease of exposition, we omit the
mathematical details here and describe them in Section S1.2.

Layer-specific PCA of PDFs thus maps each PDF to a Euclidean vector of PC
scores, which efficiently captures the variability in the voxel-intensity
PDFs. The vector of PC scores is representative of tumor heterogeneity, and
constitutes the imaging phenotype, which is the response in our Bayesian
model. For each imaging sequence $m$ and tumor layer $t$, we perform PCA using
the sample of PDFs $f^{(t,m)}_{1},\ldots,f^{(t,m)}_{n}$ to obtain the PC
scores $Y^{(t,m)}$. The number of PCs included is decided such that the
included PCs cumulatively explain 99% of the total variance. Note that unlike
the case for traditional functional data on linear Hilbert spaces, the
corresponding principal directions of variability for PDFs are not orthogonal
on the space of PDFs—they are orthogonal in the linear tangent space at the
mean. Thus, one needs to be careful when interpreting the PC scores as
capturing magnitude associated with PDF eigenfunctions; details in Section S1
elucidate on this. Although we consider the PDF-based PC scores as response in
our modeling framework, the methods developed in Section 3 are broadly
applicable to multivariate Euclidean responses.

##  3 Models and Methods

Following the above quantification, our data structure is of the following
format. Let
$\mbox{\boldmath{$\mathbf{y}$}}_{i}^{(t,m)}=(y_{i1}^{(t,m)},\ldots,y_{ip^{(t,m)}}^{(t,m)})$
represent a row vector of PC scores from imaging sequence $m\in M$ and layer
$t\in T=\\{1,\ldots,\tau\\}$ for subject $i$. Here, $p^{(t,m)}$ represents the
number of PCs included for layer $t$ and imaging sequence $m$. Let
$\mbox{\boldmath{$\mathbf{y}$}}_{i}^{(t)}=(\mbox{\boldmath{$\mathbf{y}$}}_{i}^{(t,T1)},\ldots,\mbox{\boldmath{$\mathbf{y}$}}_{i}^{(t,FLAIR)})$
represent the PC scores concatenated across all four imaging sequences and
define
$Y^{(t)}=({\mbox{\boldmath{$\mathbf{y}$}}_{1}^{(t)}}^{\top},\ldots,{\mbox{\boldmath{$\mathbf{y}$}}_{n}^{(t)}}^{\top})^{\top}\in\mathbb{R}^{n\times
p^{(t)}}$ where $p^{(t)}=\sum_{m\in M}p^{(t,m)}$ is the total number of PCs
included across all four imaging sequences. A pictorial representation of this
augmentation is provided in Figure 3.

Figure 3: Pictorial representation of $Y^{(t)}$ by augmenting PC scores across
imaging sequences.

Let $X\in\mathbb{R}^{n\times g}$ correspond to the predictors from $g$ genomic
markers for $n$ subjects. For each $t\in T$, we model the PC scores $Y^{(t)}$
using the genomic markers $X$ as

$\displaystyle Y^{(t)}$ $\displaystyle\sim$ $\displaystyle
N_{n,p^{(t)}}(XB^{(t)},I_{n}\otimes\Delta^{(t)}),\text{ or equivalently},$
$\displaystyle vec(Y^{(t)})$ $\displaystyle\sim$ $\displaystyle
N_{np^{(t)}}(vec(XB^{(t)}),I_{n}\otimes\Delta^{(t)}),$ (1)

where $B^{(t)}=[B^{(t,T1)},\ldots,B^{(t,FLAIR)}]\in\mathbb{R}^{g\times
p^{(t)}}$ and
$B^{(t,m)}=[\mbox{\boldmath$\beta$}_{1}^{(t,m)};\ldots;\mbox{\boldmath$\beta$}_{p^{(t,m)}}^{(t,m)}]\in\mathbb{R}^{g\times
p^{(t,m)}}$ for all $m\in M$. Here,
$\mbox{\boldmath$\beta$}_{j}^{(t,m)}=(\beta_{1j}^{(t,m)},\ldots,\beta_{gj}^{(t,m)})^{\top}$
for all $j=1,\ldots,p^{(t,m)}$ and $m\in M$, and $vec(\cdot)$ is the row-wise
vectorization of a matrix. The coefficient $\beta_{kj}^{(t,m)}$ corresponds to
the association of genomic marker $k$ and principal component $j$ from the
imaging sequence $m$ in the layer $t$. Also, $I_{n}$ is the identity matrix of
dimension $n$, and $\otimes$ denotes the Kronecker product. Note that
$\Delta^{(t)}\in\mathbb{R}^{p^{(t)}\times p^{(t)}}$ is a covariance matrix
that captures the covariance between the PCs across imaging sequences, which
naturally exists as these sequences are imaging the same region/tumor. In
other words, although the PC scores corresponding to a specific imaging
sequence are uncorrelated (by construction), cross-sequence correlation could
still exist (e.g. between PC scores from T1 and PC scores from T2); this
correlation is captured through $\Delta^{(t)}$. We can succinctly rewrite
$vec(XB^{(t)})=Z\mbox{\boldmath$\beta$}^{(t)}$, where $Z=X\otimes
I_{p^{(t)}}\in\mathbb{R}^{np^{(t)}\times gp^{(t)}}$ (Gupta and Nagar,, 2018)
and $\mbox{\boldmath$\beta$}^{(t)}=vec(B^{(t)})\in\mathbb{R}^{gp^{(t)}}$. Let
$\mbox{\boldmath{$\mathbf{y}$}}^{(t)}=vec(Y^{(t)})$ and
$\Theta^{(t)}=I_{n}\otimes\Delta^{(t)}$; then, Equation (1) is equivalent to
$\mbox{\boldmath{$\mathbf{y}$}}^{(t)}\sim
N_{np^{(t)}}(Z\mbox{\boldmath$\beta$}^{(t)},\Theta^{(t)})$.

Structured variable selection prior. Our aim is to identify genomic markers
associated with the PC scores, which translates to identifying the nonzero
coefficients $B^{(t)}$ in the model in Equation (1). The genomic markers could
themselves be correlated or have an inherent dependence structure (e.g.
similar biological function or signaling pathway). We incorporate this
information within the variable selection framework by proposing a
generalization of structured continuous spike-and-slab priors as follows
(George and McCulloch,, 1997; Andersen et al.,, 2014):

$\displaystyle\beta_{kj}^{(t,m)}\Big{|}\zeta_{k}^{(t,m)},{\nu^{2}_{kj}}^{(t,m)}$
$\displaystyle\stackrel{{\scriptstyle ind}}{{\sim}}$ $\displaystyle
N\big{(}0,\big{[}(1-\zeta_{k}^{(t,m)})v_{0}+\zeta_{k}^{(t,m)}v_{1}\big{]}{\nu^{2}_{kj}}^{(t,m)}\big{)}\leavevmode\nobreak\
\leavevmode\nobreak\ \forall\leavevmode\nobreak\ k,j,t,m,$
$\displaystyle\zeta_{k}^{(t,m)}\Big{|}\lambda_{k}^{(t,m)}$
$\displaystyle\stackrel{{\scriptstyle ind}}{{\sim}}$ $\displaystyle
Ber(\Phi(\lambda_{k}^{(t,m)}))\leavevmode\nobreak\ \leavevmode\nobreak\
\forall\leavevmode\nobreak\ k,t,m,$ (2)
$\displaystyle\mbox{\boldmath$\lambda$}^{(t,m)}$
$\displaystyle\stackrel{{\scriptstyle ind}}{{\sim}}$ $\displaystyle
N(\mbox{\boldmath$\mu$}^{(t,m)},\Lambda)\leavevmode\nobreak\
\leavevmode\nobreak\ \forall\leavevmode\nobreak\ t,m,$

where
${\eta^{2}_{kj}}^{(t,m)}:=\big{[}(1-\zeta_{k}^{(t,m)})v_{0}+\zeta_{k}^{(t,m)}v_{1}\big{]}{\nu^{2}_{kj}}^{(t,m)}$
is the hypervariance of $\beta_{kj}^{(t,m)}$ with $0<v_{0}\ll v_{1}$. The
indicator $\zeta_{k}^{(t,m)}$ takes values $1$ and $0$ with probability
$\Phi(\lambda_{k}^{(t,m)})$ and $(1-\Phi(\lambda_{k}^{(t,m)}))$, respectively.
Here $\Phi(\cdot)$ is the cumulative distribution function of a standard
normal distribution. Note that any suitable injective function could be chosen
instead of $\Phi$. We use a Gamma prior on ${\nu^{-2}_{kj}}^{(t,m)}$ with
$a_{1}$ and $a_{2}$ as the shape and rate parameters, respectively. If
$\zeta_{k}^{(t,m)}=1$, the hypervariance for $\beta_{kj}^{(t,m)}$ is dictated
by the prior on $v_{1}{\nu^{-2}_{kj}}^{(t,m)}$; when $\zeta_{k}^{(t,m)}=0$,
the hypervariance $v_{0}{\nu^{2}_{kj}}^{(t,m)}$ is small in magnitude allowing
for shrinkage of the coefficient $\beta_{kj}^{(t,m)}$ toward zero. For this
purpose, $v_{0}$ is chosen to be a much smaller positive number than $v_{1}$,
encouraging the exclusion of much smaller nonzero effects. We use an inverse-
Wishart (IW) prior on the covariance matrix $\Delta^{(t)}$ as
$\Delta^{(t)}\stackrel{{\scriptstyle iid}}{{\sim}}IW(\delta,\Psi)$, where
$\delta$ is the degrees of freedom and $\Psi$ is the scale matrix.

The indicator $\zeta_{k}^{(t,m)}$ corresponds to the genomic marker $k$ across
all PC scores $j$. That is, at each level $t$ and imaging modality $m$,
genomic marker $k$ is assumed to be correlated with the entire voxel intensity
PDF, and not just some of its features. This structure is built-in with the
assumption that we are interested in associations of genomic markers with
imaging sequences as a whole, and not just associations with specific PCs
within an imaging sequence.

Borrowing strength using $\mbox{\boldmath$\mu$}^{(t,m)}$ and $\Lambda$. In
essence, our model breaks down the components of a usual multiple-multivariate
regression, into a sequence of such separate models for each layer $t\in T$.
The dependence structure between the responses and hence, the corresponding
regression coefficients is captured through the hyperparameters in our model
specification. For $t=1$, we choose a fairly vague prior by considering
$\mbox{\boldmath$\mu$}^{(1,m)}=\mbox{\boldmath{$\mathbf{0}$}}$. However, for
the model at layer $t>1$, we incorporate information from the estimation for
models at the previous layers $1,\ldots,(t-1)$. We do so through a sequential
posterior updating of the layer-specific parameters. Let
$\hat{\mbox{\boldmath$\lambda$}}^{(s,m)}$ be the posterior estimate of
$\mbox{\boldmath$\lambda$}^{(s,m)}$ for $s<t$; then, some potential choices
for $\mbox{\boldmath$\mu$}^{(t,m)}$ include

  1. 1.

$\mbox{\boldmath$\mu$}^{(t,m)}=\alpha\hat{\mbox{\boldmath$\lambda$}}^{(t-1,m)}+(1-\alpha)\mbox{\boldmath{$\mathbf{0}$}}$,
where $\alpha\in[0,1]$;

  2. 2.

$\mbox{\boldmath$\mu$}^{(t,m)}=\sum\limits_{s=1}^{t-1}\alpha_{s}\hat{\mbox{\boldmath$\lambda$}}^{(t-s,m)}+\alpha_{t}\mbox{\boldmath{$\mathbf{0}$}}$,
such that $\sum\limits_{s=1}^{t}\alpha_{s}=1$;

  3. 3.

$\mbox{\boldmath$\mu$}^{(t,m)}=\sum\limits_{s=1}^{t-1}\alpha_{s}\hat{\mbox{\boldmath$\lambda$}}^{(t-s,m)}$
for $\alpha_{s}\in\mathbb{R}$;

  4. 4.

$\mbox{\boldmath$\mu$}^{(t,m)}=\alpha\hat{\mbox{\boldmath$\lambda$}}^{(t-1,m)}_{+}$,
where $\hat{\lambda}^{(t-1,m)}_{k,+}=\max(\hat{\lambda}^{(t-1,m)}_{k},0)$ for
$k=1,\ldots,g$ and $\alpha\in[0,1]$.

The first choice centers the prior on $\mbox{\boldmath$\lambda$}^{(t,m)}$
based only on the estimates from the previous layer $(t-1)$. The second choice
centers the prior in the convex hull created by estimates from all previous
layers $1,\ldots,(t-1)$. The third choice induces an auto-regressive structure
on the estimates from all previous layers $1,\ldots,(t-1)$. The fourth choice
places higher prior probability of selection on those columns of $X$ which
were selected during the estimation for the previous layer $(t-1)$. That is,
predictor $k$ will have a higher chance of selection ($\zeta_{k}^{(t,m)}=1$)
in layer $t$ for imaging sequence $m$ when $\lambda_{k}^{(t,m)}>0$. However,
$\lambda_{k}^{(t,m)}>0$ has higher probability of occurrence if
$\mu_{k}^{(t,m)}>0$. The construction of $\mbox{\boldmath$\mu$}^{(t,m)}$ in
the fourth choice above ensures $\mu_{k}^{(t,m)}\geq 0$. Note that this choice
does not penalize the selection of the predictor $k$ in case it was not
selected for the previous layer $(t-1)$. In Figure 4, for different choices of
$\alpha$, we plot the marginal prior selection probability
$\Phi(\lambda^{(t,m)}_{k})$ (evaluated at
$\lambda^{(t,m)}_{k}=\mu^{(t,m)}_{k}=\alpha\hat{\lambda}^{(t-1,m)}_{k,+}$) of
gene $k$ for imaging sequence $m$ at layer $t$, versus the marginal posterior
selection probability $\Phi(\hat{\lambda}^{(t-1,m)}_{k})$ at layer $(t-1)$
when $\Lambda_{kk}=1$. For example, if $\hat{\lambda}^{(t-1,m)}_{k}=1$, then
predictor $k$ has marginal posterior selection probability
$\Phi(\hat{\lambda}^{(t-1,m)}_{k})=0.84$ for imaging sequence $m$ at layer
$(t-1)$. In this case, the marginal prior selection probability
$\Phi(\lambda^{(t,m)}_{k})$ at layer $t$ is $0.60,0.69$ and $0.84$ when
$\alpha=0.25,0.5$ and $1$, respectively. The marginal prior selection
probability $\Phi(\lambda^{(t,m)}_{k})=0.5$ at layer $t$ when
$\hat{\lambda}^{(t-1,m)}_{k}\leq 0$.

Figure 4: Marginal prior selection probability $\Phi(\lambda^{(t,m)}_{k})$
(evaluated at
$\lambda^{(t,m)}_{k}=\mu^{(t,m)}_{k}=\alpha\hat{\lambda}^{(t-1,m)}_{k,+}$)
versus the marginal posterior selection probability
$\Phi(\hat{\lambda}^{(t-1,m)}_{k})$ when $\Lambda_{kk}=1$.

We focus on the fourth case for the simulations and case study. Nevertheless,
our estimation process proceeds in a sequential manner in $t$, where the
estimation for model at layer $t$ borrows information from the model for layer
$(t-1)$, and possibly other previous layers. In the prior specification in
Equation (2), the parameter $\lambda_{k}^{(t,m)}$ serves as a structural
parameter as its prior can be used to incorporate any dependence in the
columns of $X$. The correlation structure between the genomic markers is
incorporated into the prior for
$\mbox{\boldmath$\lambda$}^{(t,m)}=(\lambda_{1}^{(t,m)},\ldots,\lambda_{g}^{(t,m)})^{\top}$
through its covariance specification $\Lambda$. If a certain genomic marker
$k$ is associated with the PC scores from imaging sequence $m$ and layer $t$,
then we assume that other positively correlated genomic markers might also
have a similar association. This assumption is driven by the fact that voxels
from the boundaries of the spherical layers are contiguous objects with
potentially similar characteristics. This can be incorporated through a graph
Laplacian or other choice of network structure-based covariance matrix
$\Lambda$. If $\mbox{\boldmath$\mu$}^{(t,m)}=\mbox{\boldmath{$\mathbf{0}$}}$
and $\Lambda=I_{g}$, then the joint prior on $\zeta^{(t,m)}_{k}$ and
$\mbox{\boldmath$\lambda$}_{(t,m)}$ reduces to
$\zeta_{k}^{(t,m)}\stackrel{{\scriptstyle
ind}}{{\sim}}Ber(\upsilon^{(t,m)}_{k})$ with
$\upsilon^{(t,m)}_{k}\stackrel{{\scriptstyle iid}}{{\sim}}U(0,1)$.

Role of $\mbox{\boldmath$\mu$}^{(t,m)}$ and $\Lambda$. While
$\mbox{\boldmath$\mu$}^{(t,m)}$ incorporates the layered dependence structure
between responses for the selection of regression coefficients, $\Lambda$
incorporates the correlation between the columns of $X$, and hence places
higher prior probability of joint selection for any two correlated columns $k$
and $k^{\prime}$. On the other hand, $\mbox{\boldmath$\mu$}^{(t,m)}$ induces
prior dependence for the same column across various layers $t$. In other
words, $\mbox{\boldmath$\mu$}^{(t,m)}_{k}$ centers the prior probability of
selection for column $k$ in layer $t$, based on its selection in previous
layers. As a result, $\mbox{\boldmath$\mu$}^{(t,m)}$ and $\Lambda$ contribute
complementary information to the model and jointly contribute towards the
selection of the regression coefficients.

##  4 Estimation and Model Selection

The full posterior distribution corresponding to the parameters of the model
in Equation (1) with the prior structure considered in Equation (2) is given
in Section S2.1. Due to the multiple-multivariate regression setup, we have
$g\times p^{(t)}$ coefficients ($\beta_{kj}$) and/or $|M|\times g$ posterior
inclusion probabilities to sequentially estimate for each model at layer $t$.
A full MCMC-based estimation can be employed, but the computation time would
be additive with each additional spherical layer included, as the estimation
is required to be sequential based on the model formulation. As a faster
alternative to MCMC, we consider the EM-based strategy to variable selection
as proposed by Ročková and George, (2014). This EM-based estimation approach
is a deterministic alternative to commonly used stochastic search algorithms,
which exploits the prior formulation to rapidly find posterior modes. We
emphasize that the EM-based estimation presented in this work is an adaptation
of Ročková and George, (2014) to the multiple-multivariate regression setting.
We discuss sequential estimation using EM and describe the model selection
procedure below.

###  4.1 Sequential EM Updating

Our implementation of the EM-based estimation approach indirectly maximizes
the parameter posterior given by
$\pi(\beta_{kj}^{(t,m)},\mbox{\boldmath$\lambda$}^{(t,m)},{\nu^{-2}_{kj}}^{(t,m)},\Delta^{(t)}|\mbox{\boldmath{$\mathbf{y}$}}^{(t)})$,
by proceeding iteratively over the complete-data log posterior
$\log\pi(\beta_{kj}^{(t,m)},\zeta_{k}^{(t,m)},\mbox{\boldmath$\lambda$}^{(t,m)},{\nu^{-2}_{kj}}^{(t,m)},\Delta^{(t)}|\mbox{\boldmath{$\mathbf{y}$}}^{(t)})$,
where $\zeta_{k}^{(t,m)}$ are treated as missing data. At each iteration, the
complete-data log posterior is replaced by its conditional expectation given
the observed data and the current parameter estimates; this forms the E-Step.
The expected complete-data log posterior is then maximized with respect to
(w.r.t.)
$\mbox{\boldmath$\chi$}^{(t)}=(\beta_{kj}^{(t,m)},\mbox{\boldmath$\lambda$}^{(t,m)},{\nu^{-2}_{kj}}^{(t,m)},\Delta^{(t)})$
for all $m,j,k$, which constitutes the M-Step. Iterating between the E-Step
and M-Step generates a sequence of parameter estimates, which converge
monotonically to the local maximum of
$\pi(\mbox{\boldmath$\chi$}^{(t)}|\mbox{\boldmath{$\mathbf{y}$}}^{(t)})$,
under appropriate regularity conditions.

The expression for the complete-data log posterior is given in Section S2.1.
For the $u$th iteration of the algorithm, the conditional expectation of the
complete-data log posterior is given by

$\displaystyle
Q\Big{(}\beta_{kj}^{(t,m)},\mbox{\boldmath$\lambda$}^{(t,m)},{\nu^{-2}_{kj}}^{(t,m)},\Delta^{(t)}|\beta_{kj}^{(t,m),(u)},\mbox{\boldmath$\lambda$}^{(t,m),(u)},{\nu^{-2}_{kj}}^{(t,m),(u)},\Delta^{(t),(u)}\Big{)}=$
$\displaystyle:=E_{\mbox{\boldmath$\zeta$}^{(t)}|\mbox{\boldmath$\chi$}^{(t),(u)}}\Big{[}\log\big{(}\pi(\beta_{kj}^{(t,m)},\zeta_{k}^{(t,m)},\mbox{\boldmath$\lambda$}^{(t,m)},{\nu^{-2}_{kj}}^{(t,m)},\Delta^{(t)}|\mbox{\boldmath{$\mathbf{y}$}}^{(t)})\big{)}\Big{]},$
(3)

where $E_{\mbox{\boldmath$\zeta$}^{(t)}|\mbox{\boldmath$\chi$}^{(t),(u)}}$ is
a conditional expectation with
$\mbox{\boldmath$\chi$}^{(t),(u)}=(\beta_{kj}^{(t,m),(u)},\mbox{\boldmath$\lambda$}^{(t,m),(u)},{\nu^{-2}_{kj}}^{(t,m),(u)},\Delta^{(t),(u)})$
for all $m,j,k$ and $\mbox{\boldmath$\zeta$}^{(t)}=(\zeta_{k}^{(t,m)})_{m,k}$.
For the $(u+1)$-th iteration in the EM algorithm, we (i) find the expectation
in Equation (3) to obtain the $Q$-function, and (ii) maximize the $Q$-function
w.r.t. $\mbox{\boldmath$\chi$}^{(t)}$ to obtain new estimates
$\mbox{\boldmath$\chi$}^{(t),(u+1)}$.

####  4.1.1 The E-Step

As the conditional expectation is computed w.r.t. $\zeta^{(t,m)}_{k}$, we note
that the only two terms in the complete-data log posterior involving
$\zeta^{(t,m)}_{k}$ are the terms arising from the prior for
$\beta^{(t,m)}_{j}$ and the prior for $\zeta^{(t,m)}_{k}$. After some algebra
(see Section S2.2), we see that these two terms are given by

$-\sum\limits_{m\in
M}\sum\limits_{k=1}^{g}\sum\limits_{j=1}^{p^{(t,m)}}\frac{1}{2}\log({\nu^{2}_{kj}}^{(t,m)})+\frac{(\beta_{kj}^{(t,m)})^{2}}{2{\nu^{2}_{kj}}^{(t,m)}}E_{\mbox{\boldmath$\zeta$}^{(t)}\big{|}\mbox{\boldmath$\chi$}^{(t),(u)}}\Big{[}\frac{1}{(1-\zeta_{k}^{(t,m)})v_{0}+\zeta_{k}^{(t,m)}v_{1}}\Big{]},\
\text{and}$ $\sum\limits_{m\in
M}\sum\limits_{k=1}^{g}\log(1-\Phi(\lambda_{k}^{(t,m)}))+E_{\mbox{\boldmath$\zeta$}^{(t)}\big{|}\mbox{\boldmath$\chi$}^{(t),(u)}}[\zeta_{k}^{(t,m)}]\log\frac{\Phi(\lambda_{k}^{(t,m)})}{1-\Phi(\lambda_{k}^{(t,m)})}.$

To complete the E-Step, the two expectations we need to compute are given by

$w_{k}^{(t,m)}:=E_{\mbox{\boldmath$\zeta$}^{(t)}\big{|}\mbox{\boldmath$\chi$}^{(t),(u)}}[\zeta_{k}^{(t,m)}],\text{
and
}d_{k}^{(t,m)}:=E_{\mbox{\boldmath$\zeta$}^{(t)}\big{|}\mbox{\boldmath$\chi$}^{(t),(u)}}\Big{[}\frac{1}{(1-\zeta_{k}^{(t,m)})v_{0}+\zeta_{k}^{(t,m)}v_{1}}\Big{]}.$

Note then that
$w_{k}^{(t,m)}=P(\zeta_{k}^{(t,m)}=1\big{|}\mbox{\boldmath$\chi$}^{(t),(u)})=\frac{a_{k}^{(t,m)}}{a_{k}^{(t,m)}+b_{k}^{(t,m)}}$,
where $a_{k}^{(t,m)}$ and $b_{k}^{(t,m)}$ are defined in Equations (S1) and
(S2), respectively. To compute $d_{k}^{(t,m)}$, we note than when
$\zeta_{k}^{(t,m)}\in\\{0,1\\}$, we have
$((1-\zeta_{k}^{(t,m)})v_{0}+\zeta_{k}^{(t,m)}v_{1})^{-1}=\frac{1-\zeta_{k}^{(t,m)}}{v_{0}}+\frac{\zeta_{k}^{(t,m)}}{v_{1}}$,
almost surely. Hence,
$d_{k}^{(t,m)}=\frac{1-w_{k}^{(t,m)}}{v_{0}}+\frac{w_{k}^{(t,m)}}{v_{1}}$.

####  4.1.2 The M-Step

Once we compute $w_{k}^{(t,m)}$ and $d_{k}^{(t,m)}$ in the E-step, we focus on
maximizing the expected log posterior, i.e., the $Q$-function, w.r.t. other
parameters. The $Q$-function, which is a function of
$\mbox{\boldmath$\chi$}^{(t)}=(\beta_{kj}^{(t,m)},\mbox{\boldmath$\lambda$}^{(t,m)},{\nu^{-2}_{kj}}^{(t,m)},\Delta^{(t)})$
for all $m,j,k$, can be split into three terms as

$\displaystyle
Q(\mbox{\boldmath$\chi$}^{(t)}|\mbox{\boldmath$\chi$}^{(t),(u)})$
$\displaystyle=$ $\displaystyle
C+Q_{1}(\beta_{kj}^{(t,m)},{\nu^{-2}_{kj}}^{(t,m)},\Delta^{(t)}|\beta_{kj}^{(t,m),(u)},{\nu^{-2}_{kj}}^{(t,m),(u)},\Delta^{(t),(u)})+Q_{2}(\mbox{\boldmath$\lambda$}^{(t,m)}|\mbox{\boldmath$\lambda$}^{(t,m),(u)}),$

where $Q_{1}$ only depends on $\beta_{kj}^{(t,m)},{\nu^{-2}_{kj}}^{(t,m)}$ and
$\Delta^{(t)}$, and $Q_{2}$ depends on $\mbox{\boldmath$\lambda$}^{(t,m)}$.
The complete forms of $Q_{1}$ and $Q_{2}$ are given in Section S2.3.
Completing the EM-based estimation procedure, we maximize $Q$ w.r.t.
$\mbox{\boldmath$\chi$}^{(t)}$, i.e., we maximize $Q_{1}$ w.r.t.
$\beta_{kj}^{(t,m)},{\nu^{-2}_{kj}}^{(t,m)}$ and $\Delta^{(t)}$, and $Q_{2}$
w.r.t. $\mbox{\boldmath$\lambda$}^{(t,m)}$. The objective functions for all of
the maximization problems are provided in Section S2.4.

Maximizing $Q_{1}$ w.r.t. $\beta_{kj}^{(t,m)}$. Assuming we have the
parameters ${\nu^{2}_{kj}}^{(t,m),(u)},\Delta^{(t),(u)}$ at iteration $u$ and
the value $d_{k}^{(t,m),(u+1)}$ at iteration $(u+1)$, maximizing Equation (S3)
is equivalent to finding the maximum a-posteriori (MAP) estimate of a Bayesian
regression $\mbox{\boldmath{$\mathbf{y}$}}^{(t)}\sim
N_{np^{(t)}}(Z\mbox{\boldmath$\beta$}^{(t)},\Theta^{(t)})$ with the prior
$\mbox{\boldmath$\beta$}^{(t)}\sim
N(\mbox{\boldmath{$\mathbf{0}$}},\mbox{\boldmath$\Gamma$}^{(t)})$. Thus, the
update (MAP estimate) for $\mbox{\boldmath$\beta$}^{(t)}$ is given by

$\mbox{\boldmath$\beta$}^{(t),(u+1)}=\mbox{\boldmath$\Sigma$}^{(t),(u)}Z^{\top}(\Theta^{(t),(u)})^{-1}\mbox{\boldmath{$\mathbf{y}$}}\text{,
where
}\mbox{\boldmath$\Sigma$}^{(t),(u)}=\big{(}Z^{\top}(\Theta^{(t),(u)})^{-1}Z+(\mbox{\boldmath$\Gamma$}^{(t),(u)})^{-1}\big{)}^{-1},$

where $\mbox{\boldmath$\Gamma$}^{(t)}=diag(\nu_{kj}^{(t,m)}/d_{k}^{(t,m)})$
and
$\mbox{\boldmath$\Gamma$}^{(t),(u)}=diag(\nu_{kj}^{(t,m),(u)}/d_{k}^{(t,m),(u+1)})$.

Maximizing $Q_{1}$ w.r.t. ${\nu^{-2}_{kj}}^{(t,m)}$. Each term inside the
product in Equation (S4) is the kernel of a Gamma distribution for
${\nu^{-2}_{kj}}^{(t,m)}$. Assuming we have the parameters
$\mbox{\boldmath$\beta$}^{(t),(u+1)},\Delta^{(t),(u)}$ at iteration $(u+1)$,
the update for ${\nu^{-2}_{kj}}^{(t,m)}$ (maximizer of Equation (S4)) for each
$m,k,j$ is given by

${\nu^{-2}_{kj}}^{(t,m),(u+1)}=\Big{(}a_{1}-\frac{1}{2}\Big{)}\Big{/}\Big{(}a_{2}+\frac{(\beta_{kj}^{(t,m),(u+1)})^{2}d_{k}^{(t,m),(u+1)}}{2}\Big{)},\text{
where }a_{1}>\frac{1}{2}.$

Maximizing $Q_{1}$ w.r.t. $\Delta^{(t)}$. Equation (S5) is the kernel of a
Wishart distribution for $(\Delta^{(t)})^{-1}$; thus, its maximizer is the
mode of the corresponding Wishart distribution. Assuming we have the
parameters $B^{(t),(u+1)}$ at iteration $(u+1)$, the update for
$(\Delta^{(t)})^{-1}$ is given by

$(\Delta^{(t),(u+1)})^{-1}=(n+\delta+p^{(t)}+1)\big{[}\Psi+(Y^{(t)}-XB^{(t)})^{\top}(Y^{(t)}-XB^{(t)})\big{]}^{-1}.$

Maximizing $Q_{2}$ w.r.t. ${\mbox{\boldmath$\lambda$}}^{(t,m)}$. The posterior
for ${\mbox{\boldmath$\lambda$}}^{(t,m)}$ does not arise from a standard
distribution. Thus, we resort to numerical methods, specifically gradient
descent, to find the solution for the optimization in Equation (S6), to
compute the update for ${\mbox{\boldmath$\lambda$}}^{(t,m)}$ at iteration
$(u+1)$. The first derivative of the objective function
$F({\mbox{\boldmath$\lambda$}}^{(t,m)})$ is given by

$\displaystyle\nabla F({\mbox{\boldmath$\lambda$}}^{(t,m)})$ $\displaystyle=$
$\displaystyle({\mbox{\boldmath$\lambda$}^{(t,m)}-\mbox{\boldmath$\mu$}^{(t,m)}})^{\top}\Lambda^{-1}-(s_{1}^{(t,m)},\ldots,s_{g}^{(t,m)}),\text{
where}$ $\displaystyle s_{k}^{(t,m)}$ $\displaystyle=$
$\displaystyle-\frac{(1-w_{k}^{(t,m)})\phi(\lambda_{k}^{(t,m)})}{1-\Phi(\lambda_{k}^{(t,m)})}+\frac{w_{k}^{(t,m)}\phi(\lambda_{k}^{(t,m)})}{\Phi(\lambda_{k}^{(t,m)})}.$

Then, the update at the $(r+1)$-th iteration of the gradient descent algorithm
is given by

${\mbox{\boldmath$\lambda$}}^{(t,m),(u+1,r+1)}={\mbox{\boldmath$\lambda$}}^{(t,m),(u+1,r)}-\kappa\nabla
F({\mbox{\boldmath$\lambda$}}^{(t,m),(u+1,r)}),$

where $\kappa$ is the learning rate. This learning rate can be modified within
each iteration to speed up convergence. The initial value
${\mbox{\boldmath$\lambda$}}^{(t,m),(u+1,0)}$ at iteration $(u+1)$ of the EM
algorithm can be set as the solution of Equation (S6) at iteration $u$. This
solution for ${\mbox{\boldmath$\lambda$}}^{(t,m)}$ is used to assign the prior
mean for ${\mbox{\boldmath$\lambda$}}^{(t+1,m)}$ as
$\mbox{\boldmath$\mu$}^{(t+1,m)}=\alpha\hat{\mbox{\boldmath$\lambda$}}^{(t,m)}_{+}$.

One of the major drawbacks of the EM algorithm is its potential to get trapped
in a local maximum. A possible solution for this is to run the algorithm with
different starting values. In our estimation approach, we consider the
deterministic annealing variant of the EM algorithm (Ueda and Nakano,, 1998),
referred to as the DAEM algorithm, which improves the chance of finding the
global maximum (see Section S2.5 for details). All of the results presented in
this paper are based on the DAEM algorithm. We denote the DAEM algorithm-based
estimates for the parameters in our model by
$\hat{\mbox{\boldmath$\chi$}}^{(t)}$.

###  4.2 Model Selection

Once we have the estimates of the parameters
$\hat{\mbox{\boldmath$\chi$}}^{(t)}$, we can compute the posterior inclusion
probabilities for predictor $k$ in layer $t$ and imaging sequence $m$ as
$\hat{w}^{(t,m)}_{k}=P(\zeta^{(t,m)}_{k}=1|\hat{\mbox{\boldmath$\chi$}}^{(t)})$
using Equation (S8). In terms of variable selection, we set
$\hat{\zeta}^{(t,m)}_{k}=1$ when
$P(\zeta^{(t,m)}_{k}=1|\hat{\mbox{\boldmath$\chi$}}^{(t)})>0.5$, and set
$\hat{\zeta}^{(t,m)}_{k}=0$ otherwise. Note that $v_{0}$ and $v_{1}$ are fixed
hyperparameters and their choice determines the solution set
$\hat{\mbox{\boldmath$\zeta$}}^{(t)}$. The speed of the EM algorithm makes it
feasible to obtain solution sets $\hat{\mbox{\boldmath$\zeta$}}^{(t)}$ for
multiple choices of $v_{0}$ while fixing $v_{1}$. As $v_{0}$ is increased, we
see that more coefficients with smaller magnitude get absorbed by the spike
distribution, leading to sparser models. We consider a grid of values for
$v_{0}$, fix $v_{1}$ to be significantly larger than $v_{0}$, and obtain
solutions for each combination of $v_{0}$ and $v_{1}$. We denote the solution
set corresponding to the choice of $v_{0}$ as
$S(v_{0})=\big{\\{}(m,k,1),\ldots,(m,k,p^{(t,m)})|\zeta^{(t,m)}_{k}=1\leavevmode\nobreak\
\forall\leavevmode\nobreak\ m\in M,k\in\\{1,\ldots,g\\}\text{ and
}t\in\\{1,\ldots,\tau\\}\big{\\}}$.

Once we have identified the set of coefficients to be included in the models,
we compute the Bayesian Information Criterion (BIC) across all layers $t$,
given by
$BIC=\sum_{t=1}^{\tau}-2\big{[}K^{(t)}\log(n)+\log\pi(Y^{(t)}|\hat{\beta}^{(t,m)}_{kj},\widehat{\Delta^{(t)}})\big{]}$,
where $K^{(t)}=\\#\\{\hat{\beta}^{(t,m)}_{kj}\neq 0|m,k,j\\}$ is the number of
nonzero coefficients included in the model for layer $t$. The final model is
based on the inclusion of coefficients, which minimize the $BIC$. We outline
the overall step-by-step procedure from constructing the PC scores to final
model selection in Algorithm 1.

Algorithm 1 Outline of the statistical framework

1:Construct layer-wise kernel density estimates $f_{i}^{(t,m)}$ as described
in Algorithm S1 for all subjects $i=1,\ldots,n$, spherical layers
$t=1,\ldots,\tau$ and MRI sequence $m\in M$.

2:Consider $X\in\mathbb{R}^{n\times g}$ as the gene expression matrix for LGG-
genes.

3:for each layer $t=1,\ldots,\tau$ do

4: for each MRI sequence $m\in M$ do

5: Compute the PC scores $Y^{(t,m)}\in\mathbb{R}^{n\times p^{(t,m)}}$ using
PCA in Algorithm S3.

6: Define
$Y^{(t)}=[Y^{(t,FLAIR)};Y^{(t,T1)};Y^{(t,T1Gd)};Y^{(t,T2)}]\in\mathbb{R}^{n\times
p^{(t)}}$.

7:for each $v_{0}\in$ Set of values to consider do

8: for each layer $t=1,\ldots,\tau$ do

9: Bayesian Modeling:  $\displaystyle Y^{(t)}\sim
N_{n,p^{(t)}}(XB^{(t)},I_{n}\otimes\Delta^{(t)});$
$\displaystyle\beta_{kj}^{(t,m)}\stackrel{{\scriptstyle
ind}}{{\sim}}N\big{(}0,\big{[}(1-\zeta_{k}^{(t,m)})v_{0}+\zeta_{k}^{(t,m)}v_{1}\big{]}{\nu^{2}_{kj}}^{(t,m)}\big{)};$
$\displaystyle\zeta_{k}^{(t,m)}\stackrel{{\scriptstyle
ind}}{{\sim}}Ber(\Phi(\lambda_{k}^{(t,m)}));$
$\displaystyle\mbox{\boldmath$\lambda$}^{(t,m)}\stackrel{{\scriptstyle
ind}}{{\sim}}N(\mbox{\boldmath$\mu$}^{(t,m)},\Lambda);$
$\displaystyle{\nu^{-2}_{kj}}^{(t,m)}\stackrel{{\scriptstyle
iid}}{{\sim}}Gamma(a_{1},a_{2});$
$\displaystyle\Delta^{(t)}\stackrel{{\scriptstyle
iid}}{{\sim}}IW(\delta,\Psi).$

10: Estimation (Section 4) with deterministic annealing:

11:a: E-Step: Compute the E-Step as described in Section 4.1.1.

12:b: M-Step: Sequentially maximize $Q_{1}$ w.r.t.
$\beta_{kj}^{(t,m)},{\nu^{-2}_{kj}}^{(t,m)},\Delta^{(t)}$ and $Q_{2}$ w.r.t.
${\mbox{\boldmath$\lambda$}}^{(t,m)}$.

13:Model Selection: Choose the model corresponding to that $v_{0}$, which
minimizes the BIC as described in Section 4.2.

##  5 Simulation Study

In this section, we demonstrate the effectiveness of our approach through
simulation studies under different scenarios. We do not directly simulate the
MRI images; instead we simulate $Y^{(t)}$ corresponding to the PC scores for
each layer $t$. To generate the matrix
$X=[\mbox{\boldmath{$\mathbf{x}$}}_{1},\ldots,\mbox{\boldmath{$\mathbf{x}$}}_{n}]^{\top}\in\mathbb{R}^{n\times
g}$ of $g$ genomic markers for $n$ subjects, we simulate
$\mbox{\boldmath{$\mathbf{x}$}}_{i}$ independently from
$N_{g}(\mbox{\boldmath{$\mathbf{0}$}},\Sigma^{(x)})$ for $i=1,\ldots,n$, where
$\Sigma^{(x)}$ is the covariance matrix which can be specified to induce
different correlation structures between the columns of $X$. We broadly
consider two cases for our simulation study: the response $Y^{(t)}$ is (i)
generated from the model, and (ii) generated directly based on fixed
coefficients $B^{(t)}$.

Case 1: $Y^{(t)}$ is generated from the model. In this case, we want to assess
the performance of the proposed approach when the response $Y^{(t)}$ is
generated directly from the model in Equation (1), with the prior structure in
Equation (2). The step-by-step procedure to generate the data is described as
an algorithm in Section S3. We construct subcases defined by different choices
of $\Psi$, which incorporates the correlation structure between the columns of
$Y^{(t)}$. For simplicity, we consider $\Psi=\sigma^{2}I_{p}$, where
$p=p^{(t)}$ for all $t$, indicating that the number of columns of $Y^{(t)}$ is
the same for all $t$. Here, $\sigma^{2}$ acts as the noise parameter that
controls the variance in the columns of $Y^{(t)}$. The dependence structure
across the layers $t$ and imaging sequences $m$ is incorporated through the
parameter $\mbox{\boldmath$\mu$}^{(t,m)}$, as described in Section 3.

We use $n=100$ samples, $p=12$ columns in each $Y^{(t)}$ with $p^{(t,m)}=3$
for all $m\in M$ and $t\in T$, $g=20$ genomic markers, and $\tau=3$ layers. We
consider two choices for the covariance in the columns of $X$: (i)
$\Sigma^{(x)}=I_{g}$, indicating no correlation between the genomic markers,
and (ii) $\Sigma^{(x)}$ as a $(2\times 2)$-block matrix, where
$\Sigma^{(x)}_{11}=\big{(}(\Sigma^{(x)}_{11})\big{)}_{1\leq i,j\leq 10}$ with
$(\Sigma^{(x)}_{11})_{ij}=10$ if $i=j$ and $9$ otherwise,
$\Sigma^{(x)}_{12}=\Sigma^{(x)}_{21}=\mbox{\boldmath{$\mathbf{0}$}}$ and
$\Sigma^{(x)}_{22}=I_{10}$. The second choice of $\Sigma^{(x)}$ indicates two
groups of genomic markers, with one group having high between-group
correlation and the other group being uncorrelated. For the noise parameter,
we consider $\sigma^{2}=1,10,20,30$. For each combination of $\sigma^{2}$ and
$\Sigma^{(x)}$, we replicate the simulation 30 times.

Table S1 presents the complete results of EM-based estimation and subsequent
model selection (as described in Section 3) for the simulated data in Case 1.
We include the true positive selection rate (TPR), the average absolute error
for the inclusion probabilities $w^{(t,m)}_{k}$, i.e.,
$E_{w}=\frac{1}{g|M|\tau}\sum_{k,m,t}|\zeta^{(t,m)}_{k}-\hat{w}^{(t,m)}_{k}|$,
and the average squared error for the estimates of $\beta$, i.e.,
$E_{\beta}=\frac{1}{gp|M|\tau}\sum_{k,j,m,t}(\beta^{(t,m)}_{kj}-\hat{\beta}^{(t,m)}_{kj})^{2}$.
In Figure 5, we plot the average TPR across 30 simulations for the four
scenarios based on different choices of $\sigma^{2}$ for the case
$\Sigma^{(x)}=I_{g}$, i.e., assuming no correlation between the genomic
markers. As the magnitude of the noise $\sigma^{2}$ increases, the TPR
decreases across all layers. In this case, we also evaluate the performance of
estimation while ignoring the dependence between the layers, by setting
$\mbox{\boldmath$\mu$}^{(t,m)}=\mbox{\boldmath{$\mathbf{0}$}}$ for all $t$;
this result is shown as model (B) in Figure 5. From the results for layers 2
and 3, we see that including dependence across layers within the prior of
$\mbox{\boldmath$\lambda$}^{(t)}$ improves the performance in terms of the TPR
for all noise levels. Note that these two scenarios perform similarly in layer
1 as we have $\mbox{\boldmath$\mu$}^{(1,m)}=\mbox{\boldmath{$\mathbf{0}$}}$.
The first panel in Table S1 shows these rates in more detail along with
$E_{w}$ and $E_{\beta}$. Since $\Sigma^{(x)}=I_{g}$ and we use $\Lambda$ to
denote the correlation matrix of $X$, i.e., $\Lambda=cor(X)$, we expect
$\Lambda\approx I_{g}$.

|  |  |
---|---|---|---
(a) $\sigma^{2}=1$ | (b) $\sigma^{2}=10$ | (c) $\sigma^{2}=20$ | (d) $\sigma^{2}=30$

Figure 5: Average TPR in each of the three layers for Case 1 when
$\Sigma^{(x)}=I_{g}$. Model (A):
$\mbox{\boldmath$\mu$}^{(t,m)}\neq\mbox{\boldmath{$\mathbf{0}$}};\Lambda=cor(X)\approx
I_{g}$. Model (B):
$\mbox{\boldmath$\mu$}^{(t,m)}=\mbox{\boldmath{$\mathbf{0}$}};\Lambda=cor(X)\approx
I_{g}$.

|  |  |
---|---|---|---
(a) $\sigma^{2}=1$ | (b) $\sigma^{2}=10$ | (c) $\sigma^{2}=20$ | (d) $\sigma^{2}=30$

Figure 6: Average TPR in each of the three layers for Case 1 when
$\Sigma^{(x)}$ is a $(2\times 2)$-block matrix. Model (A):
$\mbox{\boldmath$\mu$}^{(t,m)}\neq\mbox{\boldmath{$\mathbf{0}$}};\Lambda=cor(X)$.
Model (B):
$\mbox{\boldmath$\mu$}^{(t,m)}=\mbox{\boldmath{$\mathbf{0}$}};\Lambda=cor(X)$.
Model (C):
$\mbox{\boldmath$\mu$}^{(t,m)}\neq\mbox{\boldmath{$\mathbf{0}$}};\Lambda=I_{g}$.
Model (D):
$\mbox{\boldmath$\mu$}^{(t,m)}=\mbox{\boldmath{$\mathbf{0}$}};\Lambda=I_{g}$.

Comparative analysis. In Figure 6, we consider the case where $\Sigma^{(x)}$
has a block structure and evaluate performance of the model under four
scenarios: (A) including both $\mbox{\boldmath$\mu$}^{(t,m)}$ and $\Lambda$,
(B) ignoring only $\mbox{\boldmath$\mu$}^{(t,m)}$, (C) ignoring only
$\Lambda$, and (D) ignoring both $\mbox{\boldmath$\mu$}^{(t,m)}$ and
$\Lambda$. Ignoring $\mbox{\boldmath$\mu$}^{(t,m)}$ corresponds to choosing
$\mbox{\boldmath$\mu$}^{(t,m)}=\mbox{\boldmath{$\mathbf{0}$}}$ for all $t$;
ignoring $\Lambda$ implies choosing $\Lambda=I_{g}$. Scenario (A) refers to
incorporating the structure as described in Section 3; Scenario (B) is similar
in spirit to other variable selection approaches (Li and Zhang,, 2010;
Vannucci and Stingo,, 2010) that incorporate dependence in the covariates
through the indicator variable $\zeta$ (while ignoring the layer structure);
Scenarios (C) and (D) correspond to vanilla Bayesian variable selection for
each layer with uniform priors on selection probabilities for $\zeta$ (George
and McCulloch,, 1997). To establish comparison, we use the same data
replications across all four scenarios. We see that incorporating the
structure in $X$ through the hyperparameter $\Lambda=cor(X)$, and borrowing
information from previous layers for the mean $\mbox{\boldmath$\mu$}^{(t,m)}$,
yields better estimation performance in terms of TPR, compared to having no
structure by choosing
$\mbox{\boldmath$\mu$}^{(t,m)}=\mbox{\boldmath{$\mathbf{0}$}}\leavevmode\nobreak\
\forall\leavevmode\nobreak\ t$ or $\Lambda=I_{g}$, or both. This improvement
in performance holds true for all choices of the noise parameter $\sigma^{2}$.
Detailed results are shown in the second panel of Table S1. We have not
included the false positive rates in the results as they were less than
$0.1\%$ on average in all scenarios. We also compared the estimation under a
multivariate-multiple regression with LASSO penalty, and observed high TPRs
but also substantially higher false positive rates (for both Cases 1 and 2).
Results presented in Table S1 report the averages obtained across 30
replications and the corresponding standard deviations. Details of the various
hyperparameter choices in the data generation and model estimation processes
are given in Section S3.

Case 2: $Y^{(t)}$ is generated directly based on fixed $B^{(t)}$. In this
case, we evaluate the performance of our approach in identifying the
associations if they are indeed sequentially dependent across layers. Similar
to Case 1, we evaluate performance of the model under the four scenarios
(A)-(D). Each of these scenarios determines the choices for
$\mbox{\boldmath$\mu$}^{(t,m)}$ and $\Lambda$ as described in Case 1.

We assume $p=p^{(t)}\leavevmode\nobreak\ \forall\leavevmode\nobreak\ t$ and
consider fixed values of $B^{(t)}$ for all layers $t=1,\ldots,\tau=3$. We
further assume that $p^{(t,m)}=3$ for all $m,t$ and $g=20$. The choice of
$B^{(t)}$ is considered so that the genomic markers which are associated with
layer $t$ are based on the associations in layers $1,\ldots,(t-1)$. These
associations are built so that they do not necessarily need to hold across all
imaging sequences. Specifically, we consider $B^{(1)}=\text{block-
diag}(B^{(1)}_{15,9},\mbox{\boldmath{$\mathbf{0}$}})$, $B^{(1)}=\text{block-
diag}(B^{(2)}_{10,6},\mbox{\boldmath{$\mathbf{0}$}})$, and
$B^{(1)}=\text{block-diag}(B^{(3)}_{5,3},\mbox{\boldmath{$\mathbf{0}$}})$,
where $B^{(t)}_{r,s}$ is an $(r\times s)$-matrix whose entries are filled by
randomly sampling $rs$ values from a double-exponential distribution centered
at $0$ with scale parameter $\theta$. The value of $\theta$ acts as the effect
size, and we consider four choices: $\theta=0.7,0.8,0.9,1$. We use
$m=1,\ldots,4$ to index the four MRI sequences. Here $B^{(1)}$ indicates that
the first 15 genomic markers are associated with all three PC scores
corresponding to $m=1,2,3$. As indicated by $B^{(2)}$, the first 10 of these
15 markers are further associated with all three PC scores for $m=1,2$.
Similarly, from $B^{(3)}$ we see that the first five markers are associated
with all three PC scores for $m=1$. We construct $X$ such that the first 15
columns are positively correlated with each other and uncorrelated with the
next five (16-20) columns. The last five columns are also considered to be
positively correlated with each other. That is, we consider $\Sigma^{(x)}$ as
a $(4\times 4)$-block matrix as shown in Figure S6 of Section S3. To generate
$Y^{(t)}$, we first simulate $\Delta^{(t)}\sim
IW(\delta=p^{(t)},\Psi=\sigma^{2}I_{g})$ with $\sigma^{2}=20$, which is one of
the higher choices for the noise parameter in Case 1. We then generate
$vec(Y^{(t)})\sim N_{np^{(t)}}(vec(XB^{(t)}),I_{n}\otimes\Delta^{(t)})$ and
appropriately reshape $vec(Y^{(t)})$ to $Y^{(t)}\in\mathbb{R}^{n\times
p^{(t)}}$. This step-by-step procedure of generating the data under Case 2 is
provided in Section S3.

|  |  |
---|---|---|---
(a) $\theta=1$ | (b) $\theta=0.9$ | (c) $\theta=0.8$ | (d) $\theta=0.7$

Figure 7: Average TPR in each of the three layers for Case 2 when
$\Sigma^{(x)}$ is a $(4\times 4)$-block matrix. Model (A):
$\mbox{\boldmath$\mu$}^{(t,m)}\neq\mbox{\boldmath{$\mathbf{0}$}};\Lambda=cor(X)$.
Model (B):
$\mbox{\boldmath$\mu$}^{(t,m)}=\mbox{\boldmath{$\mathbf{0}$}};\Lambda=cor(X)$.
Model (C):
$\mbox{\boldmath$\mu$}^{(t,m)}\neq\mbox{\boldmath{$\mathbf{0}$}};\Lambda=I_{g}$.
Model (D):
$\mbox{\boldmath$\mu$}^{(t,m)}=\mbox{\boldmath{$\mathbf{0}$}};\Lambda=I_{g}$.

We now discuss the results of the simulation study under Case 2. Similar to
Case 1, we present the TPR, $E_{w}$ and $E_{\beta}$ for each of the $\tau=3$
layers separately in Table S2. For each choice of the effect size $\theta$, we
present results corresponding to four scenarios, which are determined based on
incorporating (or ignoring) the structure through
$\mbox{\boldmath$\mu$}^{(t,m)}$ and $\Lambda$. Figure 7 shows the average TPR
under all four choices of $\theta$. We note that including the structure
through $\mbox{\boldmath$\mu$}^{(t,m)}$ and $\Lambda$ provides higher values
for TPR thereby improving performance, even for smaller effect sizes. From
Table S2, we see that by including this structure we also obtain lower values
for error metrics $E_{w}$ and $E_{\beta}$. We note that including the
dependence structure in the columns of $X$, through $\Lambda$, has a
pronounced influence on the performance compared to ignoring it by considering
$\Lambda=I_{g}$. Including the structure between the layers through
$\mbox{\boldmath$\mu$}^{(t,m)}$ also has an incremental effect in the
performance for all choices of $\theta$. The results for layer 1 are similar
for scenarios with $\mbox{\boldmath$\mu$}^{(t,m)}=0$ and
$\mbox{\boldmath$\mu$}^{(t,m)}\neq 0$ for fixed $\Lambda$ as
$\mbox{\boldmath$\mu$}^{(1,m)}=0$. Table S2 reports the averages obtained
across 30 replications, and are accompanied by the standard deviations; false
positive rates are not reported as they were negligible across all scenarios.
Details of the hyperparameter choices are provided in Section S3.

##  6 Radiogenomic Analysis in Lower Grade Gliomas

We leverage the radiology-based imaging data from TCIA and the genomic data
from The Cancer Genome Atlas (TCGA; www.cancer.gov/tcga) in the context of
LGG, to understand the associations between imaging phenotypes and genomic
markers. In this section, we describe the data acquisition and pre-processing
steps, and present results from our Bayesian modeling approach.

Imaging data. The imaging data acquired from TCIA contain publicly available
pre-operative multi-institutional MRI scans for the TCGA LGG cohort.
Specifically, we obtained MRI scans and the corresponding tumor segmentation
labels for 65 LGG subjects from Bakas et al., 2017a ; Bakas et al., 2017b .
The segmentation labels were constructed using an automated segmentation
method called GLISTRboost (Bakas et al.,, 2015). One of the hurdles posed by
MRI for statistical analysis is that the voxel intensity values are not
comparable either between study visits for the same subject or across
different subjects. Hence, it is required to pre-process these images to
perform intensity value normalization. We employ a biologically motivated
normalization technique using the _R_ package _WhiteStripe_ for intensity
normalization (Shinohara et al.,, 2014). However, we discard imaging data for
two subjects due to issues with segmentation masks and intensity value
normalization, leaving us with a total of 63 LGG subjects with imaging data.

Imaging phenotypes. Gliomas contain various heterogeneous subregions broadly
categorized into three groups: necrosis and non-enhancing core, edema, and
enhancing core. These subregions reflect differences in tumor biology, have
variable histologic and genomic phenotypes, and exhibit highly variable
clinical prognosis (Bakas et al., 2017b, ). Typically, the area of dead cells
in the center is referred to as necrosis or non-enhancing core, and the tissue
swelling caused in the outer region due to accumulation of fluid is referred
to as edema (mayfieldclinic.com/pe-braintumor.htm). Following suit, we divide
the tumor region for each subject into $\tau=3$ spherical shells, i.e., one
inner sphere and two spherical shells. For each subject, we follow the
procedure described in Algorithm S1 to construct the PDFs corresponding to the
three spherical shells. The PDFs estimated in each layer are used to construct
the PC scores $Y^{(t)}$, which act as the imaging phenotypes.

Genomic data. We obtained the genomic data from LinkedOmics (Vasaikar et al.,,
2017), a publicly available portal that includes multi-omics data from
multiple cancer types in TCGA. The genomic data we acquire are normalized
gene-level RNA sequencing data from the Illumina HiSeq system (high-throughput
sequencing) with the expression values on the $\log_{2}$ scale. This data is
obtained for all of the 63 matched LGG subjects across 20086 genes. We focus
our analysis on LGG-genes, which are genes specifically reported as cancer
drivers in LGG (Bailey et al.,, 2018). The cancer driver genes were cataloged
as driver genes based on a PanCancer and PanSoftware analysis (comprising of
all 33 TCGA projects and several computational tools). The set of cancer
driver genes for LGG includes the 24 genes reported in Table S3, whose
expression profile was retrieved from the expression of the full set of genes
from LinkedOmics. We work with the gene expression data corresponding to these
24 LGG-genes as the genomic markers for the 63 subjects.

Prior elicitation. For our statistical analysis of the LGG data, we choose
$\mbox{\boldmath$\mu$}^{(t,m)}=\alpha\hat{\mbox{\boldmath$\lambda$}}^{(t-1,m)}_{+}$,
which prioritizes selection of genes in the previous layer, and $\Lambda$ as
the correlation matrix of $X$ capturing dependence structure between the
genes. The choices of some of the other hyperparameters include
$v_{1}=100,a_{1}=4,a_{2}=5,\alpha=0.5,\Psi=I_{g},\delta=p^{(t)}$. As $v_{1}$
is expected to be very large, we have compensated by placing an informative
prior on $\nu^{-2(t,m)}_{kj}$ with more mass around the value $1$ based on our
choice of $a_{1}$ and $a_{2}$. The choice of $\delta$ and $\Psi$ are made such
that we have a non-informative prior around $\Delta^{(t)}$. Within the
estimation for each model for layer $t$, we choose $10^{-5}$ as the threshold
for convergence. The choice of $\alpha$ is such that we place the prior mean
at the midpoint of zero and $\hat{\mbox{\boldmath$\lambda$}}^{(t-1,m)}$ for
the selected genes. Further details about hyperparameter choices are provided
in Section S4.4. For model selection, we estimate the model on a grid of
values for $v_{0}\in\\{0.001+(i-1)0.001|i=1,\ldots,40\\}$ and compute the BIC
based on each value of $v_{0}$ to choose the best model as the one with lowest
BIC.

We implement the multiple-multivariate Bayesian regression model with the gene
expression of the LGG-genes as predictors and the PC scores from the MRI
images as responses. Figure 8 presents the results of our analysis, where each
panel corresponds to one of the four imaging sequences (FLAIR, T1, T1Gd and
T2). The inner sphere from the tumor region is represented as the inner circle
(layer 1). Similarly, the two spherical shells are represented as concentric
annuli (layers 2 and 3). The circular region in the plot is divided into 24
sectors, with each representing one of the 24 LGG-genes under consideration.
Consequently, each sector defined by a gene is divided into three concentric
blocks corresponding to the three layers. If $\hat{\zeta}^{(t,m)}_{k}=1$ for
gene $k$ in layer $t$ of imaging sequence $m$, the corresponding block is
colored in blue. Based on our estimation and model selection, if
$\hat{\zeta}^{(t,m)}_{k}=1$, the expression of gene $k$ has significant
association with the PC scores corresponding to layer $t$ from the imaging
sequence $m$. For example, in FLAIR, the block in layer 3 corresponding to the
sector for gene ATRX is colored, indicating significant association of the
expression of ATRX with PC scores from layer 3 in this imaging sequence.

---

Figure 8: Associations of layer-wise PC scores from the four imaging sequences
with the gene expression. Each panel corresponds to an imaging sequence, and a
block filled with color indicates association for that gene with one of the
three specific layers ($\tau=3$).

Figure S13 shows plots corresponding to the $Q$-function maximization in the
EM algorithm. In terms of computation time, our algorithm requires about $20$
seconds for estimation and model selection across $\tau=3$ layers with $g=24$
genes and $p^{(1)}=22,p^{(2)}=18,p^{(3)}=17$ PC scores. The analysis was
executed on a computer with an Intel(R) Core(TM) i7-7700 CPU @ 3.60GHz, 3601
Mhz, 4 cores, 8 logical processors with 32 GB RAM; model selection for $40$
values of $v_{0}$ was done in parallel across $4$ processors. Under the same
settings, generating 200 MCMC samples of the model parameters for each layer
sequentially would require about 50 seconds. We demonstrate the utility of our
EM-based estimation compared to a MCMC-based approach. Details of this
comparison under different choices for the number of parameters are provided
in Section S5.

Biological findings and implications. We now describe the biological
significance of the radiogenomic associations identified by our model. These
results are presented separately for each imaging sequence in Figure 8. Most
of the radiogenomic associations appear in FLAIR, which presents most of the
tumor regions in LGG with relative hypointensity (dark signal on MRI), except
for a hyperintense (bright) rim (Weerakkody et al.,, 2020). The other three
sequences usually present a hypointense signal compared to other brain
tissues, with exceptions at certain times. We specifically focus on the
following identified associations:

  1. 1.

The genes IDH1 and IDH2 demonstrate significant associations with all three
layers for the FLAIR sequence. Several recent studies (Yan et al.,, 2009; Leu
et al.,, 2013, 2016) have highlighted the prognostic importance of IDH
mutations in gliomas as well as their distinctive genetic and clinical
characteristics. Mutations that affect IDH1 were identified in more than 70%
of low grade gliomas and in glioblastomas that developed from lower-grade
lesions; most of the other tumors often had mutations affecting the gene IDH2
(Yan et al.,, 2009). These mutations were also shown to be strong predictors
for survival, specifically in the case of LGG (Leu et al.,, 2013). During the
developmental stage of most LGG tumors, one of the earliest genetic events is
the acquisition of IDH1 or IDH2 mutation (Leu et al.,, 2016).

  2. 2.

The PI3K pathway (group of genes), which includes genes such as PTEN, PIK3CA
and PIK3R1, is one of the most deregulated and druggable pathways in human
cancer (Arafeh and Samuels,, 2019). We see associations of the genes PTEN,
PIK3CA and PIK3R1 with the spherical shells in both T2 and FLAIR. The PI3K
pathway controls multiple cellular processes including metabolism, motility,
proliferation, growth and survival (Janku et al.,, 2018). Somatic mutations in
PIK3R1 act as oncogenic driver events in gliomas and are known to increase
signaling through the PI3K pathway; they also promote tumorigenesis of primary
normal human astrocytes in an orthotopic xenograft model (Quayle et al.,,
2012).

  3. 3.

Another significant association across all of the three spherical shells in
the FLAIR imaging sequence, and partially in T2, is with the genes ARID1A and
ARID2. Studies have shown that somatic mutations in the chromatin remodeling
gene ARID1A occur in several tumor types (Jones et al.,, 2012). ARID1A is
known to be a tumor suppressor and inhibits glioma cell proliferation via the
PI3K pathway (Zeng et al.,, 2013). ARID2 is a gene associated with chromatin
organization and has been predicted to be a glioma driver (Ceccarelli et al.,,
2016). Moreover, recent studies have identified ARID2 as a direct target and
functional effector in other cancers such as oral cancer (Wu et al.,, 2020).

  4. 4.

The genes PTPN11 and NIBPL also present significant associations across all
three spherical shells in the FLAIR sequence. Additionally, NIPBL presents
association across all spherical shells in T2. NIPBL is a somatically altered
glioma gene that is known to be a crucial adherin subunit, and is essential
for loading cohesins on chromatin (Ceccarelli et al.,, 2016). PTPN11 is known
to mediate gliomagenesis in both mice and humans (Liu et al.,, 2011).

  5. 5.

NRAS is a member of the RAS oncogene family, which encodes small enzymes
involved in cellular signal transduction (Fiore et al.,, 2016). In our
analysis, NRAS shows significant association with the exterior layers of the
tumor in FLAIR. RAS pathways, when activated, trigger downstream signaling
pathways such as MAPKs and PI3K/AKT, which modulate cell growth and survival
(Schubbert et al.,, 2007). NRAS was also identified as a gene which promoted
oncogenesis in glioma stem cell (Gong et al.,, 2016).

  6. 6.

We see significant associations of the gene CIC with the inner and middle
layers in T1, T1Gd and T2. CIC is a transcriptional repressor that counteracts
activation of genes downstream of receptor tyrosine kinase (RTK)/RAS/ERK
signaling pathways (Bunda et al.,, 2019). Loss of CIC potentiated the
formation and reduced the latency in tumor development in an orthotopic mouse
model of glioma (Yang et al.,, 2017). Tumorigenesis in high-grade gliomas is
attributed to hyperactive RTK/RAS/ERK signaling, while CIC has a tumor
suppressive function (Bunda et al.,, 2019). This indicates an inverse
functional relationship between CIC and the genes in RTK/RAS/ERK pathways
(e.g., NRAS). Interestingly, we see that associations with CIC are picked only
when there are no significant associations with NRAS.

These findings indicate several associations between the LGG-genes, and
related pathways, with imaging phenotypes from the spherical shells. It is
promising that several of these genes have been found to play important roles
in the genesis of glioma. It is not surprising that we see most of the
associations in the FLAIR sequence, as it exhibits the best contrast between
normal brain tissue and presumed infiltrating tumor margins (Grier and
Batchelor,, 2006). It is also the principal imaging sequence for assessment of
LGG growth (Bynevelt et al.,, 2001). These findings encourage a deeper
corroboration, which could potentially indicate intricate characteristics of
tumor growth.

In the analysis presented above, the tumor region was divided into three
spherical shells to conform with the number of natural subregions: necrosis,
edema and enhancing tumor. However, we also deployed our framework with other
choices for the number of spherical shells including $\tau=4,5,$ and $6$. We
plot these results in Figures S10-S12 in Section S4.2 of the supplementary
material. We see that the results are reasonably robust, i.e., most of the
genes associated with the layers from the imaging sequences are consistent
(except for a few additions/deletions) across different choices of $\tau$.
This indicates that our approach is broadly able to identify the underlying
radiogenomic associations between the radiomic phenotypes and gene expression.

##  7 Discussion and Future Work

In this paper, we propose a new framework for integrating radiological imaging
and genomics data that is biologically-informed, mimics the tumor evolution
process, and accounts for (structural) tumor heterogeneity. In the context of
LGG, we identify molecular determinants of this tumor evolution process by
borrowing information sequentially from concentric spherical shells; this
process tries to emulate tumor growth since 3D spheroids are known to better
mimic real tumors (Breslin and O’Driscoll,, 2013). We build layer-wise
sequential Bayesian regression models that consider (i) the PC scores
constructed from the PDFs of the MRI-based tumor voxel intensities as the
responses, which effectively capture tumor heterogeneity, and (ii) the
transcriptomic profiling data in terms of the gene expression of LGG-genes as
predictors. Information about the multiple imaging sequences, multiple PC
scores from a given sequence, and the correlation between the genes was
encoded into the prior structure of the regression model. The estimation is
based on an computationally-efficient EM algorithm, and it identifies genes
which have significant association with tumor heterogeneity in each spherical
shell. We have proposed novel methodology in terms of multiple-multivariate
regression to simultaneously incorporate complex dependence structure between
(i) the genomic covariates, and (ii) the components of the multivariate
response. Our model incorporates a sequential variable selection strategy to
borrow information across layers of the tumor. Furthermore, this sequential
strategy can be deployed to model similar structured data, e.g., any imaging
data that mimics a growth process or data where the sequencing arises due to
time.

An important aspect of our analysis is that we borrow information about the
selection of genes starting from the interior sphere and moving in an outward
direction. This approach captures the process of tumor growth which starts as
a single cancerous cell and grows radially outward to form the entire tumor.
The choice of spherical shapes for the different tumor layers is an
uninformative one, due to the lack of prior knowledge about the tumor’s growth
process. If the tumor’s growth dynamics are available, this new information
can be used to build the layers. This growth process is incorporated through
the parameter $\mbox{\boldmath$\mu$}^{(t,m)}$, necessitating sequential
estimation. However, we are only interested in the posterior inclusion
probabilities for each of the LGG-genes to identify the radiogenomic
associations. Hence, instead of using MCMC sampling techniques, we perform
estimation using an EM-algorithm, which iteratively searches for the posterior
mode of the inclusion probabilities; this, in turn, reduces computational time
of estimation. Information about the selection of genes across spherical
shells is propagated, and many significant associations in the inner sphere
carry their effect forward to subsequent spherical shells. Clearly, this is
not always the case, which is reassuring, as these results are not dominated
by our choice of prioritizing the previously selected genes through
$\mbox{\boldmath$\mu$}^{(t,m)}$. Additionally, dependence structure between
the genes is incorporated using $\Lambda$. We see that gene pairs, such as
NIBPL and ARID2 or NF1 and ZBTB20, have high positive correlation (Figure S9),
and are either jointly selected or jointly not selected in most layers and
imaging sequences. Similarly, the genes ATRX and CIC have negative
correlation, and their selection alternates, i.e., one of them is not selected
while the other one is selected.

Our analysis identifies several associations between image-based tumor
characteristics and expression profiles of LGG-genes. The genes associated
with the imaging phenotypes in the inner spherical layers can act as potential
biomarkers for early events. Monitoring these selected genes from the onset,
or the time of diagnosis, could lead to a better understanding of the
molecular underpinnings as well as provide effective diagnostic options prior
to invasive approaches such as a biopsy. One of the limitations of our data is
that the gene expression might not be derived from the spherical layers
separately during the biopsy. This can be potentially addressed in the future
using more modern technologies that facilitate analyses using spatial
transcriptomics (Burgess,, 2019) \- an area which is still maturing. In the
proposed model setup, the parameter space is high-dimensional due to the
multivariate regression setting. This results in a high computational burden
during estimation while incorporating an extremely large number of genomic
markers within the sequential multivariate regression framework; this issue
needs to be explored further using highly scalable or approximate estimation
algorithms.

Future work. Our work tries to identify radiogenomic associations between
radiomic phenotypes and genomic markers in LGGs by incorporating the
biological structure of the tumor into the modeling framework. However, this
could be further explored in several directions such as identifying the tumor
regions/spherical shells based on known segmentation of the tumor subregions.
This could additionally be informed by the shape of the tumor, which the PDFs
do not capture, as they only assess the heterogeneity in the tumor voxel
values. Also, our model explores linear relationships between the imaging-
based PC scores and the gene expression, which could potentially be extended
to study nonlinear associations (e.g. using semi-parametric additive models
(Scheipl et al.,, 2012; Ni et al.,, 2015)). Our future work will include
integrating these radiogenomic findings into predictive models for clinical
outcomes such as overall survival, time to progression.

## Acknowledgments

All of the authors acknowledge support by the NIH-NCI grant R37-CA214955. SM
was partially supported by Precision Health at The University of Michigan
(U-M). SM and AR were partially supported by CCSG P30 CA046592, U-M
Institutional Research Funds and ACS RSG-16-005-01. SK was partially supported
by NSF grants CCF-1740761, CCF-1839252 and DMS-2015226. SK and KB were
partially supported by NSF DMS-1613054. KB was partially supported by NSF
DMS-2015374. VB was supported by NIH grants R01-CA160736, R21-CA220299, and
P30 CA46592, NSF grant 1463233, and start-up funds from the U-M Rogel Cancer
Center and School of Public Health.

## References

  * Andersen et al., (2014) Andersen, M. R., Winther, O., and Hansen, L. K. (2014).  Bayesian inference for structured spike and slab priors.  In Advances in Neural Information Processing Systems, pages 1745–1753.
  * Arafeh and Samuels, (2019) Arafeh, R. and Samuels, Y. (2019).  PIK3CA in cancer: The past 30 years.  In Seminars in Cancer Biology. Elsevier.
  * Bailey et al., (2018) Bailey, M. H., Tokheim, C., Porta-Pardo, E., Sengupta, S., Bertrand, D., Weerasinghe, A., Colaprico, A., Wendl, M. C., Kim, J., Reardon, B., et al. (2018).  Comprehensive characterization of cancer driver genes and mutations.  Cell, 173(2):371–385.
  * (4) Bakas, S., Akbari, H., Sotiras, A., Bilello, M., Rozycki, M., Kirby, J., Freymann, J., Farahani, K., and Davatzikos, C. (2017a).  Segmentation labels and radiomic features for the pre-operative scans of the TCGA-LGG collection.  https://doi.org/10.7937/K9/TCIA.2017.GJQ7R0EF.
  * (5) Bakas, S., Akbari, H., Sotiras, A., Bilello, M., Rozycki, M., Kirby, J. S., Freymann, J. B., Farahani, K., and Davatzikos, C. (2017b).  Advancing the cancer genome atlas glioma MRI collections with expert segmentation labels and radiomic features.  Nature Scientific Data, 4:170117.
  * Bakas et al., (2015) Bakas, S., Zeng, K., Sotiras, A., Rathore, S., Akbari, H., Gaonkar, B., Rozycki, M., Pati, S., and Davatzikos, C. (2015).  GLISTRboost: combining multimodal MRI segmentation, registration, and biophysical tumor growth modeling with gradient boosting machines for glioma segmentation.  In International Workshop on Brainlesion: Glioma, Multiple Sclerosis, Stroke and Traumatic Brain Injuries, pages 144–155. Springer.
  * Bhattacharyya, (1943) Bhattacharyya, A. (1943).  On a measure of divergence between two statistical populations defined by their probability distributions.  Bulletin of the Calcutta Mathematical Society, 35:99–109.
  * Breslin and O’Driscoll, (2013) Breslin, S. and O’Driscoll, L. (2013).  Three-dimensional cell culture: the missing link in drug discovery.  Drug Discovery Today, 18(5-6):240–249.
  * Bunda et al., (2019) Bunda, S., Heir, P., Metcalf, J., Li, A. S. C., Agnihotri, S., Pusch, S., Yasin, M., Li, M., Burrell, K., Mansouri, S., et al. (2019).  CIC protein instability contributes to tumorigenesis in glioblastoma.  Nature Communications, 10(1):1–17.
  * Burgess, (2019) Burgess, D. J. (2019).  Spatial transcriptomics coming of age.  Nature Reviews Genetics, 20(6):317–317.
  * Bynevelt et al., (2001) Bynevelt, M., Britton, J., Seymour, H., MacSweeney, E., Thomas, N., and Sandhu, K. (2001).  Flair imaging in the follow-up of low-grade gliomas: time to dispense with the dual-echo?  Neuroradiology, 43(2):129–133.
  * Ceccarelli et al., (2016) Ceccarelli, M., Barthel, F. P., Malta, T. M., Sabedot, T. S., Salama, S. R., Murray, B. A., Morozova, O., Newton, Y., Radenbaugh, A., Pagnotta, S. M., et al. (2016).  Molecular profiling reveals biologically discrete subsets and pathways of progression in diffuse glioma.  Cell, 164(3):550–563.
  * Cencov, (1982) Cencov, N. N. (1982).  Statistical decision rules and optimal inference.  American Mathematical Society.
  * Clark et al., (2013) Clark, K., Vendt, B., Smith, K., Freymann, J., Kirby, J., Koppel, P., Moore, S., Phillips, S., Maffitt, D., Pringle, M., et al. (2013).  The cancer imaging archive (tcia): maintaining and operating a public information repository.  Journal of digital imaging, 26(6):1045–1057.
  * Claus et al., (2015) Claus, E. B., Walsh, K. M., Wiencke, J. K., Molinaro, A. M., Wiemels, J. L., Schildkraut, J. M., Bondy, M. L., Berger, M., Jenkins, R., and Wrensch, M. (2015).  Survival and low-grade glioma: the emergence of genetic information.  Neurosurgical Focus, 38(1):E6.
  * Costa et al., (2016) Costa, E. C., Moreira, A. F., de Melo-Diogo, D., Gaspar, V. M., Carvalho, M. P., and Correia, I. J. (2016).  3D tumor spheroids: an overview on the tools and techniques used for their analysis.  Biotechnology Advances, 34(8):1427–1441.
  * Dryden et al., (1998) Dryden, I. L., Mardia, K. V., et al. (1998).  Statistical shape analysis.
  * Elshafeey et al., (2019) Elshafeey, N., Kotrotsou, A., Hassan, A., Elshafei, N., Hassan, I., Ahmed, S., Abrol, S., Agarwal, A., El Salek, K., Bergamaschi, S., et al. (2019).  Multicenter study demonstrates radiomic features derived from magnetic resonance perfusion images identify pseudoprogression in glioblastoma.  Nature Communications, 10(1):1–9.
  * Fiore et al., (2016) Fiore, D., Donnarumma, E., Roscigno, G., Iaboni, M., Russo, V., Affinito, A., Adamo, A., De Martino, F., Quintavalle, C., Romano, G., et al. (2016).  mir-340 predicts glioblastoma survival and modulates key cancer hallmarks through down-regulation of NRAS.  Oncotarget, 7(15):19531.
  * George and McCulloch, (1997) George, E. I. and McCulloch, R. E. (1997).  Approaches for Bayesian variable selection.  Statistica Sinica, pages 339–373.
  * Gevaert et al., (2014) Gevaert, O., Mitchell, L. A., Achrol, A. S., Xu, J., Echegaray, S., Steinberg, G. K., Cheshier, S. H., Napel, S., Zaharchuk, G., and Plevritis, S. K. (2014).  Glioblastoma multiforme: exploratory radiogenomic analysis by using quantitative image features.  Radiology, 273(1):168–174.
  * Gong et al., (2016) Gong, W., Zheng, J., Liu, X., Ma, J., Liu, Y., and Xue, Y. (2016).  Knockdown of NEAT1 restrained the malignant progression of glioma stem cells by activating microRNA let-7e.  Oncotarget, 7(38):62208–62223.
  * Grier and Batchelor, (2006) Grier, J. T. and Batchelor, T. (2006).  Low-grade gliomas in adults.  The Oncologist, 11(6):681–693.
  * Gupta and Nagar, (2018) Gupta, A. K. and Nagar, D. K. (2018).  Matrix variate distributions.  Chapman and Hall/CRC.
  * Ishwaran et al., (2005) Ishwaran, H., Rao, J. S., et al. (2005).  Spike and slab variable selection: frequentist and Bayesian strategies.  The Annals of Statistics, 33(2):730–773.
  * Janku et al., (2018) Janku, F., Yap, T. A., and Meric-Bernstam, F. (2018).  Targeting the PI3K pathway in cancer: are we making headway?  Nature Reviews Clinical Oncology, 15(5):273–291.
  * Jones et al., (2012) Jones, S., Li, M., Parsons, D. W., Zhang, X., Wesseling, J., Kristel, P., Schmidt, M. K., Markowitz, S., Yan, H., Bigner, D., et al. (2012).  Somatic mutations in the chromatin remodeling gene ARID1A occur in several tumor types.  Human Mutation, 33(1):100–103.
  * Just, (2014) Just, N. (2014).  Improving tumour heterogeneity MRI assessment with histograms.  British Journal of Cancer, 111(12):2205–2213.
  * Karcher, (1977) Karcher, H. (1977).  Riemannian center of mass and mollifier smoothing.  Communications on Pure and Applied Mathematics, 30(5):509–541.
  * Kass and Vos, (2011) Kass, R. E. and Vos, P. W. (2011).  Geometrical foundations of asymptotic inference, volume 908.  John Wiley & Sons.
  * Kurtek and Bharath, (2015) Kurtek, S. and Bharath, K. (2015).  Bayesian sensitivity analysis with the Fisher-Rao metric.  Biometrika, 102(3):601–616.
  * Lang, (2012) Lang, S. (2012).  Fundamentals of differential geometry.  Springer Science & Business Media.
  * Leu et al., (2016) Leu, S., von Felten, S., Frank, S., Boulay, J.-L., and Mariani, L. (2016).  IDH mutation is associated with higher risk of malignant transformation in low-grade glioma.  Journal of Neuro-oncology, 127(2):363–372.
  * Leu et al., (2013) Leu, S., von Felten, S., Frank, S., Vassella, E., Vajtai, I., Taylor, E., Schulz, M., Hutter, G., Hench, J., Schucht, P., et al. (2013).  IDH/MGMT-driven molecular classification of low-grade glioma is a strong predictor for long-term survival.  Neuro-oncology, 15(4):469–479.
  * Li and Zhang, (2010) Li, F. and Zhang, N. R. (2010).  Bayesian variable selection in structured high-dimensional covariate spaces with applications in genomics.  Journal of the American statistical association, 105(491):1202–1214.
  * Liu et al., (2011) Liu, K.-W., Feng, H., Bachoo, R., Kazlauskas, A., Smith, E. M., Symes, K., Hamilton, R. L., Nagane, M., Nishikawa, R., Hu, B., et al. (2011).  SHP-2/PTPN11 mediates gliomagenesis driven by PDGFRA and INK4A/ARF aberrations in mice and humans.  The Journal of Clinical Investigation, 121(3):905–917.
  * Lu et al., (2018) Lu, C.-F., Hsu, F.-T., Hsieh, K. L.-C., Kao, Y.-C. J., Cheng, S.-J., Hsu, J. B.-K., Tsai, P.-H., Chen, R.-J., Huang, C.-C., Yen, Y., et al. (2018).  Machine learning–based radiomics for molecular subtyping of gliomas.  Clinical Cancer Research, 24(18):4429–4436.
  * Mazurowski et al., (2017) Mazurowski, M. A., Clark, K., Czarnek, N. M., Shamsesfandabadi, P., Peters, K. B., and Saha, A. (2017).  Radiogenomics of lower-grade glioma: algorithmically-assessed tumor shape is associated with tumor genomic subtypes and patient outcomes in a multi-institutional study with The Cancer Genome Atlas data.  Journal of Neuro-oncology, 133(1):27–35.
  * Mohammed et al., (2021) Mohammed, S., Bharath, K., Kurtek, S., Rao, A., and Baladandayuthapani, V. (2021).  Radiohead: Radiogenomic analysis incorporating tumor heterogeneity in imaging through densities.  arXiv preprint arXiv:2104.00510.
  * Morris and Baladandayuthapani, (2017) Morris, J. S. and Baladandayuthapani, V. (2017).  Statistical contributions to bioinformatics: Design, modelling, structure learning and integration.  Statistical Modelling, 17(4-5):245–289.
  * Ni et al., (2015) Ni, Y., Stingo, F. C., and Baladandayuthapani, V. (2015).  Bayesian nonlinear model selection for gene regulatory networks.  Biometrics, 71(3):585–595.
  * Ohgaki et al., (2004) Ohgaki, H., Dessen, P., Jourde, B., Horstmann, S., Nishikawa, T., Di Patre, P.-L., Burkhard, C., Schüler, D., Probst-Hensch, N. M., Maiorka, P. C., et al. (2004).  Genetic pathways to glioblastoma: a population-based study.  Cancer Research, 64(19):6892–6899.
  * Ohgaki and Kleihues, (2007) Ohgaki, H. and Kleihues, P. (2007).  Genetic pathways to primary and secondary glioblastoma.  The American Journal of Pathology, 170(5):1445–1453.
  * Quayle et al., (2012) Quayle, S. N., Lee, J. Y., Cheung, L. W. T., Ding, L., Wiedemeyer, R., Dewan, R. W., Huang-Hobbs, E., Zhuang, L., Wilson, R. K., Ligon, K. L., et al. (2012).  Somatic mutations of PIK3R1 promote gliomagenesis.  PloS One, 7(11).
  * Rao, (1992) Rao, C. R. (1992).  Information and the accuracy attainable in the estimation of statistical parameters.  In Breakthroughs in Statistics, pages 235–247. Springer.
  * Ročková and George, (2014) Ročková, V. and George, E. I. (2014).  EMVS: The EM approach to Bayesian variable selection.  Journal of the American Statistical Association, 109(506):828–846.
  * Saha et al., (2016) Saha, A., Banerjee, S., Kurtek, S., Narang, S., Lee, J., Rao, G., Martinez, J., Bharath, K., Rao, A. U., and Baladandayuthapani, V. (2016).  DEMARCATE: Density-based magnetic resonance image clustering for assessing tumor heterogeneity in cancer.  NeuroImage: Clinical, 12:132–143.
  * Scheipl et al., (2012) Scheipl, F., Fahrmeir, L., and Kneib, T. (2012).  Spike-and-slab priors for function selection in structured additive regression models.  Journal of the American Statistical Association, 107(500):1518–1532.
  * Schubbert et al., (2007) Schubbert, S., Shannon, K., and Bollag, G. (2007).  Hyperactive ras in developmental disorders and cancer.  Nature Reviews Cancer, 7(4):295–308.
  * Shinohara et al., (2014) Shinohara, R. T., Sweeney, E. M., Goldsmith, J., Shiee, N., Mateen, F. J., Calabresi, P. A., Jarso, S., Pham, D. L., Reich, D. S., and Crainiceanu, C. M. (2014).  Statistical normalization techniques for magnetic resonance imaging.  NeuroImage: Clinical, 6:9–19.
  * Silverman, (1986) Silverman, B. W. (1986).  Density estimation for statistics and data analysis.  CRC Press.
  * Srivastava et al., (2007) Srivastava, A., Jermyn, I. H., and Joshi, S. H. (2007).  Riemannian analysis of probability density functions with applications in vision.  In IEEE Conference on Computer Vision and Pattern Recognition, pages 1–8.
  * Swanson et al., (2003) Swanson, K. R., Bridge, C., Murray, J., and Alvord Jr, E. C. (2003).  Virtual and real brain tumors: using mathematical modeling to quantify glioma growth and invasion.  Journal of the Neurological Sciences, 216(1):1–10.
  * Ueda and Nakano, (1998) Ueda, N. and Nakano, R. (1998).  Deterministic annealing EM algorithm.  Neural Networks, 11(2):271–282.
  * Vannucci and Stingo, (2010) Vannucci, M. and Stingo, F. C. (2010).  Bayesian models for variable selection that incorporate biological information.  Bayesian Statistics, 9:1–20.
  * Vasaikar et al., (2017) Vasaikar, S. V., Straub, P., Wang, J., and Zhang, B. (2017).  LinkedOmics: analyzing multi-omics data within and across 32 cancer types.  Nucleic Aids Research, 46(D1):D956–D963.
  * Wang et al., (2013) Wang, W., Baladandayuthapani, V., Morris, J. S., Broom, B. M., Manyam, G., and Do, K.-A. (2013).  iBAG: integrative Bayesian analysis of high-dimensional multiplatform genomics data.  Bioinformatics, 29(2):149–159.
  * Weerakkody et al., (2020) Weerakkody, Y., Gaillard, F., et al. (2020).  Anaplastic astrocytoma: radiology reference article.  www.radiopaedia.org/articles/anaplastic-astrocytoma.  Accessed: October 6, 2020.
  * Wu et al., (2020) Wu, M., Duan, Q., Liu, X., Zhang, P., Fu, Y., Zhang, Z., Liu, L., Cheng, J., and Jiang, H. (2020).  MiR-155-5p promotes oral cancer progression by targeting chromatin remodeling gene ARID2.  Biomedicine & Pharmacotherapy, 122:109696.
  * Yan et al., (2009) Yan, H., Parsons, D. W., Jin, G., McLendon, R., Rasheed, B. A., Yuan, W., Kos, I., Batinic-Haberle, I., Jones, S., Riggins, G. J., et al. (2009).  IDH1 and IDH2 mutations in gliomas.  New England Journal of Medicine, 360(8):765–773.
  * Yang et al., (2020) Yang, H., Baladandayuthapani, V., Rao, A. U., and Morris, J. S. (2020).  Quantile function on scalar regression analysis for distributional data.  Journal of the American Statistical Association, 115(529):90–106.
  * Yang et al., (2017) Yang, R., Chen, L. H., Hansen, L. J., Carpenter, A. B., Moure, C. J., Liu, H., Pirozzi, C. J., Diplas, B. H., Waitkus, M. S., Greer, P. K., et al. (2017).  CIC loss promotes gliomagenesis via aberrant neural stem cell proliferation and differentiation.  Cancer Research, 77(22):6097–6108.
  * Zeng et al., (2013) Zeng, Y., Liu, Z., Yang, J., Liu, Y., Huo, L., Li, Z., Lan, S., Wu, J., Chen, X., Yang, K., et al. (2013).  ARID1A is a tumour suppressor and inhibits glioma cell proliferation via the PI3K pathway.  Head & Neck Oncology, 5(1):6.

Supplement for “Tumor Radiogenomics with Bayesian Layered Variable Selection”

We provide several details for the methods and results from the main paper
including (i) the construction of layer-wise radiomic phenotypes and the
computation of the principal component scores from a sample of probability
density functions (Section S1), (ii) the full posterior of the model, and its
logarithm, as well as some computational details related to the Expectation-
Maximization (EM) algorithm used for estimation (Section S2), (iii) details
related to the simulation study, (iv) details related to the gene expression
(Section S3), and some additional results for the The Cancer Genome Atlas
lower grade glioma data (Section S4), and (v) comparison of computational
burden for the EM-based estimation and model selection versus a Markov Chain
Monte Carlo sampling-based approach.

##  S1 Computation of Imaging Phenotypes

Figure S1 shows examples of axial slices of the whole brain magnetic resonance
imaging (MRI) for an lower grade glioma (LGG) patient, with the tumor region
highlighted using a red boundary. The four images correspond to the T1, T1Gd,
T2 and FLAIR imaging sequences, respectively.

|  |  |
---|---|---|---
(a) T1 | (b) T1Gd | (c) T2 | (d) FLAIR

Figure S1: Axial slice of a brain MRI for a subject with LGG for the four
imaging sequences: (a) T1, (b) T1Gd, (c) T2, and (d) FLAIR. The region inside
the red boundary corresponds to the segmented tumor region.

###  S1.1 Construction of Layer-wise Radiomic Phenotypes

The construction of kernel density estimates $f^{(t,m)}_{i}$ corresponding to
spherical shell $t\in\\{1,\ldots,\tau\\}$ from MRI sequence $m\in
M=\\{T1,T1Gd,T2,FLAIR\\}$ for subject $i$ are described in Algorithm S1. These
kernel density estimates are computed using the _density_ function in R using
the default bandwidth selection method, also commonly referred to as
Silverman’s rule, which is an optimal bandwidth for estimating Gaussian
densities (Silverman,, 1986).

Algorithm S1 Construction of layer-wise kernel density estimates

1:Input: MRI scans from $M$ sequences and corresponding tumor segmentation for
$n$ subjects.

2:Output: PDFs corresponding to $\tau$ layers from $M$ imaging sequences for
$n$ subjects.

3:for each subject $i=1,\ldots,n$ do

4: Define $N_{i}$ as the tumor volume (number of tumor voxels).

5: for each MRI sequence $m\in M$ do

6: for each layer $t=1,\ldots,\tau$ do

7: if $t=1$ then

8: Define $S^{(1)}_{i}$ as the sphere with
$\mbox{\boldmath{$\mathbf{c}$}}_{i}$ as its center and radius
$r^{(1)}_{i}=\inf\limits_{r\in\mathbb{R}^{+}}r$ such that $r$ is the radius of
all spheres centered at $\mbox{\boldmath{$\mathbf{c}$}}_{i}$ and containing at
least $N_{i}/\tau$ voxels.
---

9: Construct $f^{(1,m)}_{i}$ as the kernel density estimate based on the
voxels in sequence $m$ belonging to the sphere $S^{(1)}_{i}$.
---

10: if $t>1$ then

11: Define the spherical shell $D^{(t)}_{i}=S^{(t)}_{i}\backslash
S^{(t-1)}_{i}$, where $S^{(t)}_{i}$ is the sphere with center
$\mbox{\boldmath{$\mathbf{c}$}}_{i}$ and radius
$r^{(t)}_{i}=\inf\limits_{r\in\mathbb{R}^{+}}r$ such that $r$ is the radius of
all spheres centered at $\mbox{\boldmath{$\mathbf{c}$}}_{i}$ and containing at
least $tN_{i}/\tau$ voxels.
---

12: Construct $f^{(t,m)}_{i}$ as the kernel density estimate based on the
voxels in sequence $m$ belonging to the spherical shell $D^{(t)}_{i}$
---
.

###  S1.2 Mapping Probability Density Functions to PC Scores

The kernel density estimates are probability density functions (PDFs) and
belong to the Banach manifold of all PDFs. We focus on PDFs on the domain
$[0,1]$, however the following description can be applied to other domains
with simple adjustments. Let
$\mbox{$\mathcal{F}$}=\big{\\{}f:[0,1]\rightarrow\mathbb{R}_{>0}\big{|}\int_{0}^{1}f(x)dx=1\big{\\}}$.
We use the Fisher-Rao (F-R) Riemannian metric (Rao,, 1992; Kass and Vos,,
2011; Srivastava et al.,, 2007) on $\mathcal{F}$ to make it a Riemannian
manifold and facilitate computation. This metric has nice statistical
properties such as invariance to bijective and smooth transformations of the
PDF domain (Cencov,, 1982); it is also closely related to the Fisher
information matrix. However, the metric is difficult to employ for computation
so we skip the specific formula here for brevity. To build some basic
statistical tools, e.g., the distance between a pair of densities, Karcher
mean for a sample of densities, and PCA on the space of densities, we work on
a transformed space of PDFs using a square-root transformation
(Bhattacharyya,, 1943). Consequently, the space of PDFs becomes the positive
orthant of the unit sphere in $\mathbb{L}^{2}$, whose geometry is well-known,
and the F-R metric flattens to the standard $\mathbb{L}^{2}$ metric.

##### Distances Between PDFs.

The square-root transformation (SRT) of a PDF is defined as a function
$h=+\sqrt{f}$ (we omit the $+$ sign hereafter for notational convenience). The
inverse mapping is unique and is simply given by $f=h^{2}$ (Kurtek and
Bharath,, 2015). The space of SRTs is given by
$\mathcal{H}=\\{h:[0,1]\rightarrow\mathbb{R}_{>0}|\int_{0}^{1}h(x)^{2}dx=1\\}$,
which represents the positive orthant of a unit Hilbert sphere (Lang,, 2012).
The $\mathbb{L}^{2}$ Riemannian metric on $\mathcal{H}$ is defined as
$\langle\langle\delta h_{1},\delta h_{2}\rangle\rangle=\int_{0}^{1}\delta
h_{1}(t)\delta h_{2}(t)dt$, where $\delta h_{1},\delta h_{2}\in
T_{h}(\mathcal{H})$ and $T_{h}(\mathcal{H})=\\{\delta
h:[0,1]\rightarrow\mathbb{R}|\int_{0}^{1}h(x)\delta h(x)dx=0\\}$. The geodesic
paths and lengths can now be analytically computed due to the Riemannian
geometry of $\mathcal{H}$ equipped with the $\mathbb{L}^{2}$ metric
($\mathcal{H}$ is a Riemannian manifold and the $\mathbb{L}^{2}$ metric on
$\mathcal{H}$ corresponds to the F-R metric on $\mathcal{F}$). The geodesic
distance between $h_{1},h_{2}\in\mathcal{H}$ is simply given by
$d(h_{1},h_{2})=\theta=\cos^{-1}\big{(}\int_{0}^{1}h_{1}(x)h_{2}(x)dx\big{)}$.
Under this setup, the geodesic distance between two densities
$f_{1},f_{2}\in\mbox{$\mathcal{F}$}$, represented by their SRTs
$h_{1},h_{2}\in\mbox{$\mathcal{H}$}$ is defined as
$d(f_{1},f_{2})=d(h_{1},h_{2})$.

Algorithm S2 Sample Karcher mean of densities

1:$\bar{h}_{0}$ (initial estimate for the Karcher mean) $\leftarrow$ any one
of the densities in the sample OR the extrinsic average. Set $j\leftarrow 0$
and $\epsilon_{1},\epsilon_{2}>0$ small.

2:For $i=1,\ldots,n$ compute $u_{i}=\exp^{-1}_{\bar{h}_{j}}(h_{i})$.

3:Compute the average direction in the tangent space
$\bar{u}=\frac{1}{n}\sum_{i=1}^{n}u_{i}$.

4:if $||\bar{u}||_{L^{2}}<\epsilon_{1}$ then

5: return $\bar{h}_{j}$ as the Karcher mean.

6:else

7: $\bar{h}_{j+1}=\exp_{\bar{h}_{j}}(\epsilon_{2}\bar{u})$.

8: Set $j\leftarrow j+1$.

9: Return to step $2$.

##### Karcher Mean for a Sample of PDFs.

The geometry of the space of SRTs can be used to define an average PDF of a
sample of PDFs, which allows us to efficiently summarize and visualize the
sample. The average PDF can be computed using a generalized version of a mean
on a metric space called the Karcher mean (Karcher,, 1977; Dryden et al.,,
1998). Suppose we have $n$ pdfs $f_{1},\ldots,f_{n}$ and their corresponding
SRTs $h_{1},\ldots,h_{n}$. The sample Karcher mean $\bar{h}\in\mathcal{H}$ is
defined as the minimizer of the Karcher variance
$\rho(\bar{h})\propto\sum\limits_{i=1}^{n}d(h_{i},\bar{h})^{2}$, i.e.,
$\bar{h}=\text{argmin}_{h\in\mathcal{H}}\leavevmode\nobreak\ \rho(h)$.
Algorithm S2 presents a gradient-based approach to compute the Karcher mean on
$\mathcal{H}$. The Karcher mean of the sample PDFs is an intrinsic average
that is computed directly on $\mathcal{H}$ (or equivalently $\mathcal{F}$).
For the computation of the Karcher mean, we require tools from differential
geometry called the exponential and inverse-exponential maps. The exponential
map at a point $h_{1}\in\mathcal{H}$, denoted by
$\exp:T_{h_{1}}(\mathcal{H})\mapsto\mathcal{H}$, is defined as
$\exp_{h_{1}}(\delta h)=\cos(\|\delta h\|)h_{1}+\sin(\|\delta h\|)(\delta
h/\|\delta h\|)$, where $\|\delta h\|=\big{(}\int_{0}^{1}\delta
h(x)^{2}dx\big{)}^{1/2}$. The inverse-exponential map, denoted by
$\exp^{-1}_{h_{1}}:\mathcal{H}\mapsto T_{h_{1}}(\mathcal{H})$, is given by
$\exp^{-1}_{h_{1}}(h_{2})=(\theta/\sin(\theta))(h_{2}-h_{1}\cos(\theta))$.

##### Principal Component Analysis (PCA) for a Sample of PDFs.

Under standard settings, visualizing variability on the space of PDFs is not
straightforward. PCA is an effective method to explore the variability in the
PDFs through the primary modes of variation in the data. Note that the tangent
space is a (Euclidean) vector space, hence PCA can be implemented there, as in
standard multivariate problems. Algorithm S3 describes the computation of PCA
on the space generated by SRTs $h_{1},\ldots,h_{n}$ corresponding to the
sample of PDFs $f_{1},\ldots,f_{n}$, assuming that they have been sampled
using $m$ points on the domain $[0,1]$.

Algorithm S3 PCA on $T_{\bar{h}}(\mathcal{H})$

1:Compute the Karcher mean of $h_{1},\ldots,h_{n}$ as $\bar{h}$ using
Algorithm S2.

2:for $i=1,\ldots,n$ do

3: Compute projections ($v_{i}=\exp^{-1}_{\bar{h}}(h_{i})$) of $h_{i}$ onto
$T_{\bar{h}}(\mathcal{H})$.

4:Evaluate sample covariance matrix
$K=\frac{1}{n-1}\sum\limits_{i=1}^{n}v_{i}v_{i}^{\top}\in\mathbb{R}^{m\times
m}$.

5:Compute the SVD of $K=U\Sigma U^{\top}$.

The first $r$ columns of $U$ (denoted by $\tilde{U}\in\mathbb{R}^{m\times r}$)
span the $r$-dimensional principal subspace. We can compute the principal
component (PC) scores as $Y=V\tilde{U}$, where
$V^{\top}=[v_{1}\leavevmode\nobreak\ v_{2}\leavevmode\nobreak\
\ldots\leavevmode\nobreak\ v_{n}]\in\mathbb{R}^{m\times n}$. The PC scores in
$Y_{i}$ act as Euclidean coordinates corresponding to the density $f_{i}$, and
can be used for downstream modelling. In our case, for each imaging sequence
$m$ and each tumor layer $t$, we perform PCA with the sample of PDFs
$f^{(t,m)}_{1},\ldots,f^{(t,m)}_{n}$ to obtain the PC scores
$Y^{(t,m)}_{1},\ldots,Y^{(t,m)}_{n}$.

##  S2 Estimation Details

###  S2.1 Full Posterior and Log Posterior

The posterior distribution,
$\pi(\beta_{kj}^{(t,m)},\zeta_{k}^{(t,m)},\mbox{\boldmath$\lambda$}^{(t,m)},{\nu^{-2}_{kj}}^{(t,m)},{\theta^{-2}}^{(t,m)},\Delta^{(t)}|\mbox{\boldmath{$\mathbf{y}$}}^{(t)},X)$,
corresponding to the model in Equation (1), with the prior structure in
Equation (2), is proportional to

$\displaystyle|\Delta^{(t)}|^{-\frac{n}{2}}\exp\Big{(}-\frac{1}{2}(\mbox{\boldmath{$\mathbf{y}$}}^{(t)}-Z\mbox{\boldmath$\beta$}^{(t)})^{\top}(\Theta^{(t)})^{-1}(\mbox{\boldmath{$\mathbf{y}$}}^{(t)}-Z\mbox{\boldmath$\beta$}^{(t)})\Big{)}\times|\Delta^{(t)}|^{-\frac{\delta+p^{(t)}+1}{2}}\exp\Big{(}-\frac{1}{2}tr(\Psi(\Delta^{(t)})^{-1}\Big{)}\times$
$\displaystyle\times\prod\limits_{m\in
M}\prod\limits_{j=1}^{p^{(t,m)}}\prod\limits_{k=1}^{g}\Big{[}({\eta^{2}_{kj}}^{(t,m)})^{-1/2}\exp\Big{(}-\frac{(\beta_{kj}^{(t,m)})^{2}}{2{\eta^{2}_{kj}}^{(t,m)}}\Big{)}\Big{]}\times\prod\limits_{m\in
M}\prod\limits_{k=1}^{g}(1-\Phi(\lambda_{k}^{(t,m)}))^{(1-\zeta_{k}^{(t,m)})}\Phi(\lambda_{k}^{(t,m)})^{\zeta_{k}^{(t,m)}}\times$
$\displaystyle\times\prod\limits_{m\in
M}\exp\Big{(}-\frac{({\mbox{\boldmath$\lambda$}^{(t,m)}-\mbox{\boldmath$\mu$}^{(t,m)}})^{\top}\Lambda^{-1}(\mbox{\boldmath$\lambda$}^{(t,m)}-\mbox{\boldmath$\mu$}^{(t,m)})}{2}\Big{)}\times\prod\limits_{m\in
M}\prod\limits_{j=1}^{p^{(t,m)}}\prod\limits_{k=1}^{g}({\nu^{-2}_{kj}}^{(t,m)})^{a_{1}-1}\exp(-a_{2}{\nu^{-2}_{kj}}^{(t,m)}).$

The complete-data log posterior as described for the EM-based estimation is
given by

$\displaystyle\log\big{(}\pi(\beta_{kj}^{(t,m)},\zeta_{k}^{(t,m)},\mbox{\boldmath$\lambda$}^{(t,m)},{\nu^{-2}_{kj}}^{(t,m)},\Delta^{(t)}|\mbox{\boldmath{$\mathbf{y}$}}^{(t)},X)\big{)}=c-\frac{1}{2}(\mbox{\boldmath{$\mathbf{y}$}}^{(t)}-Z\mbox{\boldmath$\beta$}^{(t)})^{\top}(\Theta^{(t)})^{-1}(\mbox{\boldmath{$\mathbf{y}$}}^{(t)}-Z\mbox{\boldmath$\beta$}^{(t)})+$
$\displaystyle+$ $\displaystyle\sum\limits_{m\in
M}\sum\limits_{j=1}^{p^{(t,m)}}\sum\limits_{k=1}^{g}\Big{[}-\frac{1}{2}\log({\eta^{2}_{kj}}^{(t,m)})-\frac{(\beta_{kj}^{(t,m)})^{2}}{2{\eta^{2}_{kj}}^{(t,m)}}\Big{]}+$
$\displaystyle+$ $\displaystyle\sum\limits_{m\in
M}\sum\limits_{k=1}^{g}{(1-\zeta_{k}^{(t,m)})}\log(1-\Phi(\lambda_{k}^{(t,m)}))+{\zeta_{k}^{(t,m)}}\log\Phi(\lambda_{k}^{(t,m)})-$
$\displaystyle-$ $\displaystyle\sum\limits_{m\in
M}\frac{({\mbox{\boldmath$\lambda$}^{(t,m)}-\mbox{\boldmath$\mu$}^{(t,m)}})^{\top}\Lambda^{-1}(\mbox{\boldmath$\lambda$}^{(t,m)}-\mbox{\boldmath$\mu$}^{(t,m)})}{2}+\sum\limits_{m\in
M}\sum\limits_{j=1}^{p^{(t,m)}}\sum\limits_{k=1}^{g}-(a_{1}-1)\log({\nu^{2}_{kj}}^{(t,m)})-a_{2}{\nu^{-2}_{kj}}^{(t,m)}-$
$\displaystyle-$
$\displaystyle\frac{n+\delta+p^{(t)}+1}{2}\log|\Delta^{(t)}|-\frac{1}{2}tr(\Psi(\Delta^{(t)})^{-1}).$

###  S2.2 E-Step Computation

We have two terms which involve $\zeta_{k}^{(t,m)}$ and hence the expectation
$E_{\mbox{\boldmath$\zeta$}^{(t)}|\mbox{\boldmath$\chi$}^{(t),(u)}}$. The
first one is the term for the prior on $\beta_{kj}^{(t,m)}$, which is given by

$\displaystyle-\sum\limits_{m\in
M}\sum\limits_{k=1}^{g}\sum\limits_{j=1}^{p^{(t,m)}}E_{\mbox{\boldmath$\zeta$}^{(t)}\big{|}\mbox{\boldmath$\chi$}^{(t),(u)}}\Big{[}\frac{1}{2}\log({\eta^{2}_{kj}}^{(t,m)})+\frac{(\beta_{kj}^{(t,m)})^{2}}{2{\eta^{2}_{kj}}^{(t,m)}}\Big{]}=$
$\displaystyle=$ $\displaystyle-\sum\limits_{m\in
M}\sum\limits_{k=1}^{g}\sum\limits_{j=1}^{p^{(t,m)}}E_{\mbox{\boldmath$\zeta$}^{(t)}\big{|}\mbox{\boldmath$\chi$}^{(t),(u)}}\Big{[}\frac{1}{2}\log(\big{[}(1-\zeta_{k}^{(t,m)})v_{0}+\zeta_{k}^{(t,m)}v_{1}\big{]})+\frac{1}{2}\log({\nu^{2}_{kj}}^{(t,m)})+\frac{(\beta_{kj}^{(t,m)})^{2}}{2{\eta^{2}_{kj}}^{(t,m)}}\Big{]}=$
$\displaystyle=$ $\displaystyle c-\sum\limits_{m\in
M}\sum\limits_{k=1}^{g}\sum\limits_{j=1}^{p^{(t,m)}}\frac{1}{2}\log({\nu^{2}_{kj}}^{(t,m)})+\frac{(\beta_{kj}^{(t,m)})^{2}}{2{\nu^{2}_{kj}}^{(t,m)}}E_{\mbox{\boldmath$\zeta$}^{(t)}\big{|}\mbox{\boldmath$\chi$}^{(t),(u)}}\Big{[}\frac{1}{(1-\zeta_{k}^{(t,m)})v_{0}+\zeta_{k}^{(t,m)}v_{1}}\Big{]}.$

Note that we do not need to compute
$E_{\mbox{\boldmath$\zeta$}^{(t)}\big{|}\mbox{\boldmath$\chi$}^{(t),(u)}}\Big{[}\log(\big{[}(1-\zeta_{k}^{(t,m)})v_{0}+\zeta_{k}^{(t,m)}v_{1}\big{]})\Big{]}$
as it is a scalar and does not involve any other parameters. The second term
arises from the prior for $\zeta_{k}^{(t,m)}$ which is given by

$\displaystyle\sum\limits_{m\in
M}\sum\limits_{k=1}^{g}E_{\mbox{\boldmath$\zeta$}^{(t)}\big{|}\mbox{\boldmath$\chi$}^{(t),(u)}}[(1-\zeta_{k}^{(t,m)})]\log(1-\Phi(\lambda_{k}^{(t,m)}))+E_{\mbox{\boldmath$\zeta$}^{(t)}\big{|}\mbox{\boldmath$\chi$}^{(t),(u)}}[\zeta_{k}^{(t,m)}]\log\Phi(\lambda_{k}^{(t,m)})=$
$\displaystyle=$ $\displaystyle\sum\limits_{m\in
M}\sum\limits_{k=1}^{g}\log(1-\Phi(\lambda_{k}^{(t,m)}))+E_{\mbox{\boldmath$\zeta$}^{(t)}\big{|}\mbox{\boldmath$\chi$}^{(t),(u)}}[\zeta_{k}^{(t,m)}]\log\frac{\Phi(\lambda_{k}^{(t,m)})}{1-\Phi(\lambda_{k}^{(t,m)})}.$

Additionally, we have

$\displaystyle a_{k}^{(t,m)}$ $\displaystyle=$ $\displaystyle
P(\zeta_{k}^{(t,m)}=1\big{|}\lambda_{k}^{(t,m),(u)})\prod\limits_{j=1}^{p^{(t,m)}}\pi(\beta_{kj}^{(t,m),(u)}\big{|}{\nu^{2}_{kj}}^{(t,m),(u)},\zeta_{k}^{(t,m)}=1)$
(S1) $\displaystyle=$
$\displaystyle\Phi(\lambda_{k}^{(t,m),(u)})\prod\limits_{j=1}^{p^{(t,m)}}\phi(\beta_{kj}^{(t,m),(u)}\big{|}0,v_{1}{\nu^{2}_{kj}}^{(t,m),(u)}),\text{
and }$ $\displaystyle b_{k}^{(t,m)}$ $\displaystyle=$ $\displaystyle
P(\zeta_{k}^{(t,m)}=0\big{|}\lambda_{k}^{(t,m),(u)})\prod\limits_{j=1}^{p^{(t,m)}}\pi(\beta_{kj}^{(t,m),(u)}\big{|}{\nu^{2}_{kj}}^{(t,m),(u)},\zeta_{k}^{(t,m)}=0)$
(S2) $\displaystyle=$
$\displaystyle\big{(}1-\Phi(\lambda_{k}^{(t,m),(u)})\big{)}\prod\limits_{j=1}^{p^{(t,m)}}\phi(\beta_{kj}^{(t,m),(u)}\big{|}0,v_{0}{\nu^{2}_{kj}}^{(t,m),(u)}).$

Here, $\phi(.|\theta,\sigma^{2})$ denotes the PDF of a normal distribution
with mean $\theta$ and variance $\sigma^{2}$, and $\Phi(.)$ denotes the
cumulative distribution function of a standard normal distribution.

###  S2.3 $Q$-Function

The $Q$-function is separable in the parameters as $Q_{1}$ and $Q_{2}$ are
given by

$\displaystyle
Q_{1}(\beta_{kj}^{(t,m)},{\nu^{-2}_{kj}}^{(t,m)},\Delta^{(t)}|\beta_{kj}^{(t,m),(u)},{\nu^{-2}_{kj}}^{(t,m),(u)},\Delta^{(t),(u)})=$
$\displaystyle=$
$\displaystyle-\frac{1}{2}(\mbox{\boldmath{$\mathbf{y}$}}^{(t)}-Z\mbox{\boldmath$\beta$}^{(t)})^{\top}(\Theta^{(t)})^{-1}(\mbox{\boldmath{$\mathbf{y}$}}^{(t)}-Z\mbox{\boldmath$\beta$}^{(t)})-$
$\displaystyle-$ $\displaystyle\sum\limits_{m\in
M}\sum\limits_{k=1}^{g}\sum\limits_{j=1}^{p^{(t,m)}}\frac{1}{2}\log({\nu^{2}_{kj}}^{(t,m)})+\frac{(\beta_{kj}^{(t,m)})^{2}}{2{\nu^{2}_{kj}}^{(t,m)}}d_{k}^{(t,m)}+$
$\displaystyle+$ $\displaystyle\sum\limits_{m\in
M}\sum\limits_{j=1}^{p^{(t,m)}}\sum\limits_{k=1}^{g}-(a_{1}-1)\log({\nu^{2}_{kj}}^{(t,m)})-a_{2}{\nu^{-2}_{kj}}^{(t,m)}-$
$\displaystyle-$
$\displaystyle\frac{n+\delta+p^{(t)}+1}{2}\log|\Delta^{(t)}|-\frac{1}{2}tr(\Psi(\Delta^{(t)})^{-1}),\text{
and}$ $\displaystyle
Q_{2}(\mbox{\boldmath$\lambda$}^{(t,m)}|\mbox{\boldmath$\lambda$}^{(t,m),(u)})=\sum\limits_{m\in
M}\sum\limits_{k=1}^{g}\log(1-\Phi(\lambda_{k}^{(t,m)}))+E_{\mbox{\boldmath$\zeta$}^{(t)}\big{|}\mbox{\boldmath$\chi$}^{(t),(u)}}[\zeta_{k}^{(t,m)}]\log\frac{\Phi(\lambda_{k}^{(t,m)})}{1-\Phi(\lambda_{k}^{(t,m)})}-$
$\displaystyle-$ $\displaystyle\sum\limits_{m\in
M}\frac{({\mbox{\boldmath$\lambda$}^{(t,m)}-\mbox{\boldmath$\mu$}^{(t,m)}})^{\top}\Lambda^{-1}(\mbox{\boldmath$\lambda$}^{(t,m)}-\mbox{\boldmath$\mu$}^{(t,m)})}{2}.$

###  S2.4 M-Step Computation

##### Maximizing $Q_{1}$ w.r.t. $\beta_{kj}^{(t,m)}$.

This involves maximizing

$\displaystyle\max\limits_{\mbox{\boldmath$\beta$}^{(t)}}-\frac{1}{2}(\mbox{\boldmath{$\mathbf{y}$}}^{(t)}-Z\mbox{\boldmath$\beta$}^{(t)})^{\top}(\Theta^{(t)})^{-1}(\mbox{\boldmath{$\mathbf{y}$}}^{(t)}-Z\mbox{\boldmath$\beta$}^{(t)})-\sum\limits_{m\in
M}\sum\limits_{k=1}^{g}\sum\limits_{j=1}^{p^{(t,m)}}\frac{(\beta_{kj}^{(t,m)})^{2}}{2{\nu^{2}_{kj}}^{(t,m)}}d_{k}^{(t,m)}=$
(S3) $\displaystyle=$
$\displaystyle\max\limits_{\mbox{\boldmath$\beta$}^{(t)}}-\frac{1}{2}(\mbox{\boldmath{$\mathbf{y}$}}^{(t)}-Z\mbox{\boldmath$\beta$}^{(t)})^{\top}(\Theta^{(t)})^{-1}(\mbox{\boldmath{$\mathbf{y}$}}^{(t)}-Z\mbox{\boldmath$\beta$}^{(t)})-\frac{1}{2}{\mbox{\boldmath$\beta$}^{(t)}}^{\top}(\mbox{\boldmath$\Gamma$}^{(t)})^{-1}\mbox{\boldmath$\beta$}^{(t)},$

where $\mbox{\boldmath$\Gamma$}^{(t)}=diag(\nu_{kj}^{(t,m)}/d_{k}^{(t,m)})$.

##### Maximizing $Q_{1}$ w.r.t. ${\nu^{-2}_{kj}}^{(t,m)}$.

This involves maximizing

$\max\limits_{{\nu^{-2}_{kj}}^{(t,m)}}-\sum\limits_{m\in
M}\sum\limits_{k=1}^{g}\sum\limits_{j=1}^{p^{(t,m)}}\bigg{[}(a_{1}-\frac{1}{2})\log({\nu^{2}_{kj}}^{(t,m)})+(a_{2}+\frac{(\beta_{kj}^{(t,m)})^{2}d_{k}^{(t,m)}}{2}){\nu^{-2}_{kj}}^{(t,m)}\bigg{]},$

which is the same as maximizing

$\max\limits_{{\nu^{-2}_{kj}}^{(t,m)}}\prod\limits_{m\in
M}\prod\limits_{k=1}^{g}\prod\limits_{j=1}^{p^{(t,m)}}({\nu^{-2}_{kj}}^{(t,m)})^{(a_{1}-\frac{1}{2})}\exp\bigg{(}-\Big{(}a_{2}+\frac{(\beta_{kj}^{(t,m)})^{2}d_{k}^{(t,m)}}{2}\Big{)}{\nu^{-2}_{kj}}^{(t,m)}\bigg{)}.$
(S4)

##### Maximizing $Q_{1}$ w.r.t. $\Delta^{(t)}$.

This involves maximizing

$\displaystyle\max\limits_{(\Delta^{(t)})^{-1}}-\frac{1}{2}(\mbox{\boldmath{$\mathbf{y}$}}^{(t)}-Z\mbox{\boldmath$\beta$}^{(t)})^{\top}(\Theta^{(t)})^{-1}(\mbox{\boldmath{$\mathbf{y}$}}^{(t)}-Z\mbox{\boldmath$\beta$}^{(t)})-\frac{n+\delta+p^{(t)}+1}{2}\log|\Delta^{(t)}|-\frac{1}{2}tr(\Psi(\Delta^{(t)})^{-1})$
$\displaystyle=$
$\displaystyle\max\limits_{(\Delta^{(t)})^{-1}}-\frac{n+\delta+p^{(t)}+1}{2}\log|\Delta^{(t)}|-\frac{1}{2}tr(\big{[}\Psi+(Y^{(t)}-XB^{(t)})^{\top}(Y^{(t)}-XB^{(t)})\big{]}(\Delta^{(t)})^{-1}),$

which is the same as maximizing

$\max\limits_{(\Delta^{(t)})^{-1}}|(\Delta^{(t)})^{-1}|^{\frac{n+\delta+p^{(t)}+1}{2}}\exp\bigg{(}-\frac{1}{2}tr\big{(}\big{[}\Psi+(Y^{(t)}-XB^{(t)})^{\top}(Y^{(t)}-XB^{(t)})\big{]}(\Delta^{(t)})^{-1}\big{)}\bigg{)}.$
(S5)

Equation (S5) is the kernel of a Wishart distribution for
$(\Delta^{(t)})^{-1}$; hence, its maximizer is the mode of the corresponding
Wishart distribution.

##### Maximizing $Q_{2}$ w.r.t. ${\mbox{\boldmath$\lambda$}}^{(t,m)}$.

The terms in $Q_{2}$ involving ${\mbox{\boldmath$\lambda$}}^{(t,m)}$ are given
by

$\displaystyle\max\limits_{{\mbox{\boldmath$\lambda$}}^{(t,m)}}\sum\limits_{m\in
M}\sum\limits_{k=1}^{g}(1-w_{k}^{(t,m)})\log(1-\Phi(\lambda_{k}^{(t,m)}))+w_{k}^{(t,m)}\log\Phi(\lambda_{k}^{(t,m)})$
$\displaystyle+$ $\displaystyle\sum\limits_{m\in
M}-\frac{({\mbox{\boldmath$\lambda$}^{(t,m)}-\mbox{\boldmath$\mu$}^{(t,m)}})^{\top}\Lambda^{-1}(\mbox{\boldmath$\lambda$}^{(t,m)}-\mbox{\boldmath$\mu$}^{(t,m)})}{2}.$

These terms are separable in $m$ and thus let us consider the maximization
problem for each $m$, which can be termed as a minimization problem given by

$\displaystyle\min\limits_{{\mbox{\boldmath$\lambda$}}^{(t,m)}}F({\mbox{\boldmath$\lambda$}}^{(t,m)})=\min\limits_{{\mbox{\boldmath$\lambda$}}^{(t,m)}}-\sum\limits_{k=1}^{g}(1-w_{k}^{(t,m)})\log(1-\Phi(\lambda_{k}^{(t,m)}))+w_{k}^{(t,m)}\log\Phi(\lambda_{k}^{(t,m)})$
(S6) $\displaystyle+$
$\displaystyle\frac{({\mbox{\boldmath$\lambda$}^{(t,m)}-\mbox{\boldmath$\mu$}^{(t,m)}})^{\top}\Lambda^{-1}(\mbox{\boldmath$\lambda$}^{(t,m)}-\mbox{\boldmath$\mu$}^{(t,m)})}{2}.$

###  S2.5 Deterministic Annealing to Mitigate Multi-modality

In our estimation approach, we consider the deterministic annealing variant of
the EM algorithm (Ueda and Nakano,, 1998), referred to as the DAEM algorithm,
which improves the chance of finding the global maximum. Using the DAEM
algorithm, we aim to maximize a tempered version of the complete-data log
posterior, which is the same as computing the E-step as

$\tilde{Q}(\mbox{\boldmath$\chi$}^{(t)}|\mbox{\boldmath$\chi$}^{(t),(u)})=\int\log\big{(}\pi(\mbox{\boldmath$\chi$}^{(t)},\zeta_{k}^{(t,m)}|\mbox{\boldmath{$\mathbf{y}$}}^{(t)})\big{)}\frac{\pi(\mbox{\boldmath$\chi$}^{(t),(u)},\zeta_{k}^{(t,m)}|\mbox{\boldmath{$\mathbf{y}$}}^{(t)})^{q}}{\int\pi(\mbox{\boldmath$\chi$}^{(t),(u)},\zeta_{k}^{(t,m)}|\mbox{\boldmath{$\mathbf{y}$}}^{(t)})^{q}d\mbox{\boldmath$\zeta$}^{(t)}}d\mbox{\boldmath$\zeta$}^{(t)},$
(S7)

where $0<q\leq 1$ and $1/q$ corresponds to a temperature parameter which
determines the degree of separation between the multiple modes of
$\tilde{Q}(\mbox{\boldmath$\chi$}^{(t)}|\mbox{\boldmath$\chi$}^{(t),(u)})$.
Large values of temperature smooth the local modes of
$\tilde{Q}(\mbox{\boldmath$\chi$}^{(t)}|\mbox{\boldmath$\chi$}^{(t),(u)})$;
the modes begin to reappear as the temperature decreases, i.e., as $q$ gets
close to $1$. This can be incorporated in the iterative process by increasing
the value of $q$ in each iteration of the DAEM algorithm (bounded above by
$1$), which decreases the sensitivity to initialization. For our analyses, we
start with $q=0.01$ and increase it by 10% with each iteration until $q=1$.
This change in the E-step does not affect the M-step. The computation of the
E-step in Equation (S7) is straightforward: the two terms requiring evaluation
of expectation remain the same as described in Section 4.1.1. However, their
computations change as follows:

$\displaystyle\tilde{w}_{k}^{(t,m)}$ $\displaystyle=$ $\displaystyle
E_{\mbox{\boldmath$\zeta$}^{(t)}\big{|}\mbox{\boldmath$\chi$}^{(t),(u)}}[\zeta_{k}^{(t,m)}]=\frac{{a_{k}^{(t,m)}}^{q}}{{a_{k}^{(t,m)}}^{q}+{b_{k}^{(t,m)}}^{q}}\text{,
and}$ (S8) $\displaystyle\tilde{d}_{k}^{(t,m)}$ $\displaystyle=$
$\displaystyle
E_{\mbox{\boldmath$\zeta$}^{(t)}\big{|}\mbox{\boldmath$\chi$}^{(t),(u)}}\Big{[}\frac{1}{(1-\zeta_{k}^{(t,m)})v_{0}+\zeta_{k}^{(t,m)}v_{1}}\Big{]}=\frac{1-\tilde{w}_{k}^{(t,m)}}{v_{0}}+\frac{\tilde{w}_{k}^{(t,m)}}{v_{1}},$

where $a_{k}^{(t,m)}$ and $b_{k}^{(t,m)}$ are given in Equations (S1) and
(S2), respectively. Other details about the limiting behavior w.r.t. the
temperature, and the superiority of DAEM over the regular EM approach in the
context of variable selection, are discussed in Ročková and George, (2014).

##  S3 Details of Simulation Study

###  S3.1 Case 1

Algorithm S4 describes the procedure to generate the simulated data under Case
1. To generate data from the model, we make the following choices for the
hyperparameters:
$\delta=p^{(t)}=12,\Psi=\sigma^{2}I_{g},a_{1}=5,a_{2}=5,v_{0}=0.01,v_{1}=1$
and $\alpha=0.8$. We use the same choices for all cases of
$\sigma^{2}=1,10,20,30$.

For estimation, we choose
$a_{1}=4,a_{2}=5,\alpha=0.5,\Psi=I_{g},\delta=p^{(t)}$ and set $v_{1}$ as the
smallest power of 10 greater than the maximum magnitude of the entries in the
ordinary least square estimate of $B^{(t)}$. For model selection, we first
perform estimation on a grid of values for $v_{0}$. From the estimation based
on each value of $v_{0}$, we compute the Bayesian information criterion (BIC)
and choose the best model based on the lowest BIC value. In Table S1, we show
the average true positive rates (TPR), $E_{w}$ and $E_{\beta}$ after model
selection based on 30 replications. We include results for both of the choices
for $\Sigma^{(x)}$ as described in the main manuscript. Figures S2 and S5 show
the plots of $E_{w}$ across layers for various choices of $\sigma^{2}$ when
$\Sigma^{(x)}=I_{g}$ and the $(2\times 2)$-block matrix structure,
respectively. Moreover, we also include plots for $E_{\beta}$ in Figures S3
and S4 for the cases when $\Sigma^{(x)}=I_{g}$ and the $(2\times 2)$-block
matrix structure, respectively. We considered the grid of values for $v_{0}$
as $\\{0.001+(v-1)0.001|v=1,\ldots,10\\}$.

Algorithm S4 Simulation from the model

1:Specify the choices for
$n=100,p^{(t)}=12,p^{(t,m)}=3,m=4,\tau=3,g=20,\hat{\lambda}_{k}^{(t,m)}=0\leavevmode\nobreak\
\forall k,m,t$.

2:Specify the hyperparameters $\delta,\Psi,a_{1},a_{2},v_{0},v_{1},\alpha$ and
$\sigma^{2}\in\\{1,10,20,30\\}$.

3:Generate
$X=[\mbox{\boldmath{$\mathbf{x}$}}_{1},\ldots,\mbox{\boldmath{$\mathbf{x}$}}_{n}]^{\top}\in\mathbb{R}^{n\times
g}$, where $\mbox{\boldmath{$\mathbf{x}$}}_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}N_{g}(\mbox{\boldmath{$\mathbf{0}$}},\Sigma^{(x)})$ for
$i=1,\ldots,n$.

4:$\Lambda\leftarrow\Sigma^{(x)}$.

5:for $t=1,\ldots,\tau$ do

6: Generate $\Delta^{(t)}\sim IW(\delta,\Psi)$.

7: Generate ${\nu^{-2}_{kj}}^{(t,m)}\stackrel{{\scriptstyle
iid}}{{\sim}}Gamma(a_{1},a_{2})$ for all $j,k,t,m$.

8: if $t=1$ then

9: return $\mbox{\boldmath$\mu$}^{(t,m)}=\mbox{\boldmath{$\mathbf{0}$}}$.

10: else

11:
$\mbox{\boldmath$\mu$}^{(t,m)}=\alpha\hat{\mbox{\boldmath$\lambda$}}^{(t-1,m)}_{+}$,
where $\hat{\lambda}^{(t-1,m)}_{k,+}=\max(\hat{\lambda}^{(t-1,m)}_{k},0)$ for
$k=1,\ldots,g$ and $\alpha\in[0,1]$.

12: Generate $\mbox{\boldmath$\lambda$}^{(t,m)}\stackrel{{\scriptstyle
ind}}{{\sim}}N(\mbox{\boldmath$\mu$}^{(t,m)},\Lambda)$ for all $t,m$.

13: Generate $\zeta_{k}^{(t,m)}\stackrel{{\scriptstyle
ind}}{{\sim}}Ber(\Phi(\lambda_{k}^{(t,m)}))$ for all $k,t,m$.

14: Generate $\beta_{kj}^{(t,m)}\stackrel{{\scriptstyle
ind}}{{\sim}}N\big{(}0,\big{[}(1-\zeta_{k}^{(t,m)})v_{0}+\zeta_{k}^{(t,m)}v_{1}\big{]}{\nu^{2}_{kj}}^{(t,m)}\big{)}$
for all $k,j,t,m$.

15: Generate $vec(Y^{(t)})\sim
N_{np^{(t)}}(vec(XB^{(t)}),I_{n}\otimes\Delta^{(t)})$ and appropriately
reshape $Y^{(t)}\in\mathbb{R}^{n\times p^{(t)}}$.

|  |  |
---|---|---|---
(a) $\sigma^{2}=1$ | (b) $\sigma^{2}=10$ | (c) $\sigma^{2}=20$ | (d) $\sigma^{2}=30$

Figure S2: Average values of $E_{w}$ in each of the three layers for Case 1:
$\Sigma^{(x)}=I_{g}$. Model (A):
$\mbox{\boldmath$\mu$}^{(t,m)}\neq\mbox{\boldmath{$\mathbf{0}$}};\Lambda=cor(X)\approx
I_{g}$. Model (B):
$\mbox{\boldmath$\mu$}^{(t,m)}=\mbox{\boldmath{$\mathbf{0}$}};\Lambda=cor(X)\approx
I_{g}$.

|  |  |
---|---|---|---
(a) $\sigma^{2}=1$ | (b) $\sigma^{2}=10$ | (c) $\sigma^{2}=20$ | (d) $\sigma^{2}=30$

Figure S3: Average values of $E_{\beta}$ in each of the three layers for Case
1: $\Sigma^{(x)}=I_{g}$. Model (A):
$\mbox{\boldmath$\mu$}^{(t,m)}\neq\mbox{\boldmath{$\mathbf{0}$}};\Lambda=cor(X)\approx
I_{g}$. Model (B):
$\mbox{\boldmath$\mu$}^{(t,m)}=\mbox{\boldmath{$\mathbf{0}$}};\Lambda=cor(X)\approx
I_{g}$.

|  |  |
---|---|---|---
(a) $\sigma^{2}=1$ | (b) $\sigma^{2}=10$ | (c) $\sigma^{2}=20$ | (d) $\sigma^{2}=30$

Figure S4: Average values of $E_{\beta}$ in each of the three layers for Case
1: $\Sigma^{(x)}$ is a $(2\times 2)$-block matrix. Model (A):
$\mbox{\boldmath$\mu$}^{(t,m)}\neq\mbox{\boldmath{$\mathbf{0}$}};\Lambda=cor(X)$.
Model (B):
$\mbox{\boldmath$\mu$}^{(t,m)}=\mbox{\boldmath{$\mathbf{0}$}};\Lambda=cor(X)$.
Model (C):
$\mbox{\boldmath$\mu$}^{(t,m)}\neq\mbox{\boldmath{$\mathbf{0}$}};\Lambda=I_{g}$.
Model (D):
$\mbox{\boldmath$\mu$}^{(t,m)}=\mbox{\boldmath{$\mathbf{0}$}};\Lambda=I_{g}$.

|  |  |
---|---|---|---
(a) $\sigma^{2}=1$ | (b) $\sigma^{2}=10$ | (c) $\sigma^{2}=20$ | (d) $\sigma^{2}=30$

Figure S5: Average values of $E_{w}$ in each of the three layers for Case 1:
$\Sigma^{(x)}$ is a $(2\times 2)$-block matrix. Model (A):
$\mbox{\boldmath$\mu$}^{(t,m)}\neq\mbox{\boldmath{$\mathbf{0}$}};\Lambda=cor(X)$.
Model (B):
$\mbox{\boldmath$\mu$}^{(t,m)}=\mbox{\boldmath{$\mathbf{0}$}};\Lambda=cor(X)$.
Model (C):
$\mbox{\boldmath$\mu$}^{(t,m)}\neq\mbox{\boldmath{$\mathbf{0}$}};\Lambda=I_{g}$.
Model (D):
$\mbox{\boldmath$\mu$}^{(t,m)}=\mbox{\boldmath{$\mathbf{0}$}};\Lambda=I_{g}$.

|  |  |  | Layer 1 |  | Layer 2 |  | Layer 3
---|---|---|---|---|---|---|---|---
$\Sigma^{(x)}$ | $\mbox{\boldmath$\mu$}^{(t,m)},\Lambda$ | $\sigma^{2}$ |  | TPR | $E_{w}$ | $E_{\beta}$ |  | TPR | $E_{w}$ | $E_{\beta}$ |  | TPR | $E_{w}$ | $E_{\beta}$
|  |  |  |  |  |  |  |  |  |  |  |  |  |
Identity Matrix |  $\mbox{\boldmath$\mu$}^{(t,m)}\neq\mbox{\boldmath{$\mathbf{0}$}}$ $\Lambda=cor(X)$ | 1 |  | 0.901 (0.048) | 0.051 (0.024) | 0.018 (0.006) |  | 0.904 (0.039) | 0.057 (0.024) | 0.018 (0.006) |  | 0.901 (0.038) | 0.067 (0.027) | 0.018 (0.008)
10 |  | 0.804 (0.074) | 0.101 (0.04) | 0.036 (0.012) |  | 0.787 (0.075) | 0.125 (0.047) | 0.046 (0.016) |  | 0.818 (0.067) | 0.122 (0.047) | 0.046 (0.019)
20 |  | 0.746 (0.084) | 0.126 (0.045) | 0.055 (0.015) |  | 0.732 (0.065) | 0.157 (0.044) | 0.065 (0.02) |  | 0.711 (0.075) | 0.192 (0.056) | 0.079 (0.025)
30 |  | 0.669 (0.089) | 0.168 (0.056) | 0.077 (0.024) |  | 0.735 (0.057) | 0.154 (0.036) | 0.075 (0.017) |  | 0.711 (0.07) | 0.186 (0.045) | 0.089 (0.021)
|  |  |  |  |  |  |  |  |  |  |  |  |
$\mbox{\boldmath$\mu$}^{(t,m)}=\mbox{\boldmath{$\mathbf{0}$}}$ $\Lambda=cor(X)$ $\approx I_{g}$ | 1 |  | 0.901 (0.048) | 0.051 (0.024) | 0.018 (0.006) |  | 0.901 (0.04) | 0.059 (0.025) | 0.018 (0.006) |  | 0.895 (0.038) | 0.071 (0.027) | 0.019 (0.008)
10 |  | 0.804 (0.074) | 0.101 (0.04) | 0.036 (0.012) |  | 0.782 (0.076) | 0.128 (0.047) | 0.047 (0.017) |  | 0.808 (0.067) | 0.13 (0.048) | 0.049 (0.019)
20 |  | 0.746 (0.084) | 0.126 (0.045) | 0.055 (0.015) |  | 0.728 (0.064) | 0.16 (0.045) | 0.067 (0.02) |  | 0.707 (0.076) | 0.195 (0.056) | 0.081 (0.026)
30 |  | 0.669 (0.089) | 0.168 (0.056) | 0.077 (0.024) |  | 0.733 (0.056) | 0.155 (0.035) | 0.076 (0.016) |  | 0.71 (0.071) | 0.187 (0.045) | 0.09 (0.021)
|  |  |  |  |  |  |  |  |  |  |  |  |  |
Block Matrix |  $\mbox{\boldmath$\mu$}^{(t,m)}\neq\mbox{\boldmath{$\mathbf{0}$}}$ $\Lambda=cor(X)$ | 1 |  | 0.918 (0.048) | 0.04 (0.024) | 0.015 (0.007) |  | 0.922 (0.036) | 0.045 (0.021) | 0.015 (0.006) |  | 0.915 (0.043) | 0.057 (0.028) | 0.017 (0.008)
10 |  | 0.819 (0.062) | 0.091 (0.036) | 0.035 (0.013) |  | 0.81 (0.066) | 0.112 (0.041) | 0.041 (0.015) |  | 0.801 (0.071) | 0.13 (0.049) | 0.046 (0.015)
20 |  | 0.737 (0.08) | 0.132 (0.039) | 0.053 (0.016) |  | 0.734 (0.073) | 0.155 (0.042) | 0.064 (0.018) |  | 0.704 (0.066) | 0.191 (0.041) | 0.076 (0.018)
30 |  | 0.73 (0.065) | 0.134 (0.042) | 0.061 (0.017) |  | 0.699 (0.066) | 0.181 (0.044) | 0.083 (0.025) |  | 0.724 (0.06) | 0.18 (0.041) | 0.085 (0.023)
|  |  |  |  |  |  |  |  |  |  |  |  |
$\mbox{\boldmath$\mu$}^{(t,m)}=\mbox{\boldmath{$\mathbf{0}$}}$ $\Lambda=cor(X)$ | 1 |  | 0.918 (0.048) | 0.04 (0.024) | 0.015 (0.007) |  | 0.92 (0.037) | 0.046 (0.021) | 0.016 (0.007) |  | 0.91 (0.044) | 0.06 (0.028) | 0.019 (0.009)