# Structured Bayesian variable selection for
multiple related response variables and high-dimensional predictors

Zhi Zhao  Oslo Centre for Biostatistics and Epidemiology, University of Oslo,
Norway Marco Banterle  Department of Medical Statistics, London School of
Hygiene & Tropical Medicine, UK  Alex Lewin# Department of Medical Statistics,
London School of Hygiene & Tropical Medicine, UK  Manuela Zucknick# Oslo
Centre for Biostatistics and Epidemiology, University of Oslo, Norway

###### Abstract

It is becoming increasingly common to study the complex association between
multiple outcomes and high-dimensional predictors in biomedicine. However, the
related multiple outcomes and related high-dimensional predictors require
flexible and efficient joint statistical models. We propose a multivariate
structured Bayesian variable selection model to identify sparse predictors
associated with multiple outcomes. The approach uses a known structure prior
between the multiple outcomes and predictors via a Markov random field (MRF)
prior for the latent indicator variables of the coefficient matrix of a sparse
seemingly unrelated regression. The prior knowledge in the MRF prior has the
potential to generate new structure knowledge. In addition, we employ random
effects to capture sample heterogeneity of grouped data. The proposed approach
is validated by simulation studies and a large-scale pharmacogenomic study
which includes a pharmacological profiling and multi-omics data from _in
vitro_ anti-cancer drug screening.

Keywords: Bayesian multivariate regression; seemingly unrelated regression;
structured spike-and-slab prior; Markov random field prior, random effects

##  1 Introduction

Multivariate regression models involve multiple response variables that have
been used in multi-traits and genomic data analysis. For example, expression
quantitative trait loci (eQTL) analysis and analysis of large-scale drug
screening datasets have multiple correlated response variables and high-
dimensional predictors. A fundamental task is to study the
structures/associations in the integrative datasets to improve the model
prediction performance.

In this work, we analyze data from the Genomics of Drug Sensitivity in Cancer
(GDSC) database [15, 48]. The large-scale in-vitro cancer drug screens
produced a large amount of drug sensitivity data which are correlated, since
multiple drugs have a similar mechanism or the common targeted gene/pathway.
Meanwhile the genomic information of cancer tissue cell lines was measured
which is expected to guide the personalized cancer therapy [15, 5, 36, 16],
for example genome-wide measurements of mRNA expression, DNA copy numbers and
DNA mutations. These omics characterizations reflect different aspects’
information of the same system or cofunctionality of multiple genes
within/between different characterizations [27]. Therefore the analysis of
joint associations between the correlated multiple phenotypes and high-
dimensional molecular features poses both theoretical and computational
challenges.

Bayesian modelling allows us the flexibility to specify the relationships
between the complex data sets. However, classic Markov chain Monte Carlo
(MCMC) sampling becomes infeasible to search the whole huge model space when
jointly analyzing multiple phenotypes and high-dimensional genomic
information. To overcome the infeasibility of sampling the full enumerated
model space with MCMC samplers in high dimensional settings, [9] proposed an
Evolutionary Stochastic Search (ESS) algorithm based on Evolutionary Monte
Carlo. [39] extended the ESS to the computationally efficient hierarchical
evolutionary stochastic search (HESS) for multi-response regressions by
assuming independent residuals. To model dependent residuals, the seemingly
unrelated regression (SUR) provides a general framework for multi-response
regressions [2, 20, 13]. [37] showed that the inverse Wishart prior on the
covariance matrix of residuals is powerful to do simultaneous analysis.
Furthermore the hyper-inverse Wishart prior has been used to learn a sparse
graph structure for the covariance matrix of high-dimensional variables [10,
47, 6], thus performing covariance selection.

Besides imposing different structures on the covariances via the choice of
prior, it is necessary to impose variable selection priors on the high-
dimensional predictors and connect the variable selection to the covariance
selection. [24] mapped each genetic marker to multiple phenotypes (i.e., eQTL)
by using the spike-and-slab prior and hyper predictor-effect prior. [8]
proposed the hotspot prior for variable selection in multivariate regressions,
which induces both predictor and response effects. [28] used the Markov random
field (MRF) prior to encourage the joint selection of the same predictor
across several correlated response variables.

To address the large computational expense to learn the covariance structure
of multiple response variables, [20] adapted the SUR to the Bayesian framework
and used MCMC algorithms for the analytically intractable posterior inference.
[2] used a Direct Monte Carlo approach to solve the SUR models efficiently. In
high-dimensional situations, [3] presented a Bayesian sparse SUR (SSUR) model
to learn a sparse structure of the covariance matrix of the high-dimensional
responses, perform variable section and implement the ESS algorithm for
computational efficiency by assuming a underlying decomposable graph between
response variables.

In this article, we propose a structured Bayesian variable selection approach
based on the SSUR model, which employs the MRF prior on the latent indicator
variable to emphasize a joint structure between the multiple response
variables and high-dimensional predictors. This joint structure represents
known associations, for example from the literature or previous studies. For
example, Figure 1 shows two groups of the response variables (i.e., drugs) and
their corresponding two predictor groups (i.e., related target gene pathways)
in the omics data sets. The MRF prior is able to address the joint structure
by including all edges in Figure 1 that reflect the relationships between
response variables and predictors. In addition, we use random effects to
capture the sample heterogeneity, since the data set may be collected from
multiple sample groups, for example, drug sensitivity and genomic information
measured from different cancer tissues.

Figure 1: Illustration of the drug groups and omics path (adapted from [41]).

The rest of the paper is organized as follows. In Section 2, we introduce the
SUR model briefly and propose the MRF prior on the latent indicator variable.
Then the SSUR model and SUR model with the MRF prior and random effects are
presented. Section 3 shows the performance of the SSUR model and our approach
(i.e., SSUR model with MRF prior) with respect to the structure recovery and
prediction by the simulated data sets. In Section 4, we perform data analysis
from the Genomics of Drug Sensitivity in Cancer project [15]. Finally, we
discuss the main findings and conclude the paper in Section 5.

##  2 Methodology

###  2.1 Sparse SUR (SSUR) model

The multivariate regression has a response matrix $\mathbf{Y}_{n\times m}$
from $n$ samples and $m$ response variables. All response variables correspond
to the same $p$ predictors which are measured on the $n$ samples, so that the
predictor matrix is $\mathbf{X}_{n\times p}$. We first assume correlated
response variables but independent samples. The SUR model is defined as

$\displaystyle\begin{split}\mathbf{Y}&=\mathbf{X}\bm{B}+\mathbf{U},\\\
\text{vec}\\{\mathbf{U}\\}&\sim\mathcal{N}(\mathbf{0},\
\Psi\otimes\mathbb{I}_{n}),\end{split}$ (1)

where the residuals have dependent columns with covariance $\Psi$ and
independent rows, and $\text{vec}\\{\cdot\\}$ is to vectorize a matrix by
column.

In the Bayesian framework, it is difficult to sample the posteriors of the
regression coefficients from the model (1). In order to improve computation,
[49] reparametrized the SUR model and proposed a direct Monte Carlo procedure
for Bayesian analysis. They used a reparametrization proposed by [12] with
Jeffrey’s priors. [3] use the same transformation but with more general
priors. If the covariance matrix is given an inverse Wishart prior
$\Psi\sim\mathcal{IW}(\nu,\tau\mathbb{I}_{m})$, model (1) can be rewritten as

$\displaystyle\begin{split}\bm{y}_{j}&=\mathbf{X}\bm{\beta}_{j}+\sum_{l<j}\bm{u}_{l}\rho_{jl}+\bm{\epsilon}_{j},\text{\
\ }j=1,\cdots,m,\\\ \bm{\epsilon}_{j}&\sim\mathcal{N}(\mathbf{0},\
\sigma_{j}^{2}\mathbb{I}_{n}),\end{split}$ (2)

where $\bm{u}_{l}=\bm{y}_{l}-\mathbf{X}\bm{\beta}_{l}$. The reparametrized
parameters $(\sigma_{j}^{2},\rho_{jl})$ have priors

$\displaystyle\begin{split}\sigma_{j}^{2}&\sim\mathcal{IG}\left(\frac{\nu-m+2j-1}{2},\frac{\tau}{2}\right),\\\
\rho_{jl}|\sigma_{j}^{2}&\sim\mathcal{N}\left(0,\frac{\sigma_{j}^{2}}{\tau}\right),j>l,\end{split}$
(3)

where $v$ is fixed and $\tau\sim\mathcal{G}amma(a_{\tau},b_{\tau})$. [2] shows
that the conditional prior of regression coefficients after such
reparametrization is invariant when changing the order of response variables.

The transformation (3) makes parallel computing possible for the likelihood,
which benefits high-dimensional response variables especially. If there are
only a few associated predictors, but $p\gg n$, i.e., high-dimensional
predictors, we use a latent indicator matrix
$\mathbf{\Gamma}=\\{\gamma_{kj}\\}$ to show the variable selection. If
$\gamma_{kj}=1$, then $\beta_{kj}\neq 0$ and the $k$-th predictor is regarded
as an associated predictor to the $j$th response variable; otherwise
$\beta_{kj}=0$. The independent spike-and-slab priors on coefficients
$\\{\beta_{kj}:j=1,\cdots,m;k=1,\cdots,p\\}$ are able to find a small subset
of predictors that explains the variability of $\mathbf{Y}$:

$\beta_{kj}|\gamma_{kj},w\sim\gamma_{kj}\mathcal{N}(0,\
w)+(1-\gamma_{kj})\delta_{0}(\beta_{kj}),\\\ $ (4)

where $w\sim\mathcal{IG}(a_{w},b_{w})$ and $\delta_{0}(\cdot)$ is the Dirac
delta function.

An alternative of the inverse Wishart prior on the covariance $\Psi$ is a
hyper-inverse Wishart prior

$\Psi\sim\mathcal{HIW}_{\mathcal{G}}(\nu,\tau\mathbb{I}_{m}).\\\ $ (5)

It assumes an underlying decomposable graph $\mathcal{G}$ between response
variables. [17] propose a junction tree sampler for a decomposable graph
$\mathcal{G}$ which can be efficiently estimated. A sparse graph may be of
interest, so [3] specify an $Bernoulli(\eta)$ prior on each edge of the graph.
Then towards all edges of $\mathcal{G}$, a Binomial prior is on the
cardinality edge-set

$|\mathcal{E}|\sim Binomial(m(m-1)/2,\eta),$ (6)

where $\eta\sim\mathcal{B}eta(a_{\eta},b_{\eta})$ controls the sparsity of the
graph. The new reparameterization is similar to (3) but with respect to the
response variables of each clique of $\mathcal{G}$ rather than all response
variables.

Given a graph with fewer edges via (6), i.e., small $\eta$, a sparse SUR
(SSUR) model proposed by [3] achieves sparse covariance and variable selection
via a hotspot prior on the latent indicator variable.

Prior specification

For the hyperparameters $\nu$ and $\tau$ in the hyper-inverse Wishart prior
$\Psi\sim\mathcal{HIW}_{\mathcal{G}}(\nu,\tau\mathbb{I}_{m})$, we take
$\nu=m+2$, $a_{\tau}=0.1$ and $b_{\tau}=10$ by default. $\nu=m+2$ is the
smallest integer degree and makes the reparametrized parameter
$\sigma_{j}^{2}$ with the proper distribution. $\nu\gg m+2$ may enforce a weak
correlation between response variables (see Supplementary S3 for sensitivity
analysis). The model is not sensitive to the hyperparameters of $\tau$ and
$\eta$ (see Supplementary S3 for sensitivity analysis).

In order to specify the hyperparameters (i.e., $a_{w}$ and $b_{w}$) of $w$ in
(4) properly, we need to be sure sampling $w$ larger than
$\mathbb{E}[\beta_{kj}]+1.96\sqrt{\mathbb{V}\text{ar}[\beta_{kj}]}$, where
$\mathbb{E}[\cdot]$ and $\mathbb{V}\text{ar}[\cdot]$ are the mean and variance
of the prior of $\beta_{kj}$. Therefore the drawn $w$ provides enough large
value as the variance of $\beta_{kj}$ almostly, so that the prior (4) of
$\beta_{kj}$ is able to cover the scale of true effects. Since $w$ has
posterior conditional

$w|a_{w},b_{w},\mathbf{\Gamma},\bm{B}\sim\mathcal{IG}\left(a_{w}+\frac{1}{2}\sum_{kj}\gamma_{kj},b_{w}+\frac{1}{2}\sum_{k,j}\beta_{kj}^{2}\right),$
(7)

we can choose the optimal hyperparameters $a_{w}$ and $b_{w}$ according to
this posterior conditional with 5% quantile larger than
$\mathbb{E}[\beta_{kj}]+1.96\sqrt{\mathbb{V}\text{ar}[\beta_{kj}]}$.

For convenience, the hyperparameters $a_{w}$ and $b_{w}$ can be chosen with
similar scales as $\frac{1}{2}\sum_{kj}\gamma_{kj}$ and
$\frac{1}{2}\sum_{kj}\beta_{0,tj}^{2}$, respectively. For example, we assume a
model sparsity (i.e., proportion of nonzero coefficients)
$r_{sparsity}\in(0,1)$ and can either estimate the effect
$\mathbb{E}[\beta_{kj}]$ and its variation $\mathbb{V}\text{ar}[\beta_{kj}]$
from previous studies or expert knowledge about typical effect size in similar
studies. Then we specify

$\displaystyle a_{w}$ $\displaystyle=const_{a}\cdot 1/2\cdot m\cdot p\cdot
r_{sparsity},$ $\displaystyle b_{w}$ $\displaystyle=const_{b}\cdot 1/2\cdot
m\cdot p\cdot r_{sparsity}\cdot(\mathbb{E}[\beta_{kj}])^{2},$

where $const_{a}$ and $const_{b}$ are chosen to be sure that 5% quantile of
the posterior conditional (7) larger than
$\mathbb{E}[\beta_{kj}]+1.96\sqrt{\mathbb{V}\text{ar}[\beta_{kj}]}$.

###  2.2 SSUR model with MRF prior

In Figure 1, there are known relationships between the drug responses and
genomic predictors from the literature. That means we may have prior knowledge
on the regression coefficients in the multivariate regression model (1).
Instead of independent or hotspot priors [39, 29, 3] on indicator variable
$\gamma_{kj}$ of $\beta_{kj}$, we address the structures between the response
variables and predictors via an MRF prior on the latent indicator vector
$\bm{\gamma}=\text{vec}\\{\bm{\Gamma}\\}$

$f(\bm{\gamma}|d,e,G)\propto\exp\\{d\mathbbm{1}^{\top}\bm{\gamma}+e\bm{\gamma}^{\top}G\bm{\gamma}\\},$
(8)

where $d$ controls the overall model sparsity, $e$ determines the strength of
the structure relationships between responses and predictors, and $G$ is an
$mp\times mp$ matrix of the structure prior knowledge. To illustrate the idea,
we consider a simple case with four response variables ($\bm{y}_{1}$,
$\bm{y}_{2}$ and $\bm{y}_{3}$) and four predictors ($\bm{x}_{1}$,
$\bm{x}_{2}$, $\bm{x}_{3}$ and $\bm{x}_{4}$). When predictors $\bm{x}_{1}$ and
$\bm{x}_{2}$ are associated with responses $\bm{y}_{1}$ and $\bm{y}_{2}$, and
$\bm{x}_{3}$ and $\bm{x}_{4}$ are associated with $\bm{y}_{3}$, then $G$ is
the $12\times 12$ matrix given by

$G=\blockarray{ccccccccccccc}&\gamma_{11}\gamma_{21}\gamma_{31}\gamma_{41}\gamma_{12}\gamma_{22}\gamma_{32}\gamma_{42}\gamma_{13}\gamma_{23}\gamma_{33}\gamma_{43}\\\
\block{c(cccccccccccc)}\gamma_{11}010011000000\\\ \gamma_{21}100011000000\\\
\gamma_{31}000000000000\\\ \gamma_{41}000000000000\\\
\gamma_{12}110001000000\\\ \gamma_{22}110010000000\\\
\gamma_{32}000000000000\\\ \gamma_{42}000000000000\\\
\gamma_{13}000000000000\\\ \gamma_{23}000000000000\\\
\gamma_{33}000000000001\\\ \gamma_{43}000000000010\\\ \ .$

Elements of $G$ can be any real number which indicates weights for the prior
relationships between response variables and predictors. If only knowing the
relationships between response variables or relationships between predictors,
then $G=G_{y}\otimes G_{x}-\mathbbm{I}_{mp}$. For example, if we know
predictors $\bm{x}_{1}$ and $\bm{x}_{2}$ are related, $\bm{x}_{3}$ and
$\bm{x}_{4}$ are related, but do not know any related information to the three
response variables, then $G=G_{y}\otimes G_{x}-\mathbbm{I}_{12}$ where
$G_{y}=\mathbbm{I}_{3}$ and

$G_{x}=\begin{pmatrix}1&1&0&0\\\ 1&1&0&0\\\ 0&0&1&1\\\ 0&0&1&1\\\
\end{pmatrix}.$

Prior specification

It is known that the hyperparameter $e$ in (8) displays the phase transition
behaviour [30, 44] which results in an increased number of selected predictors
(i.e., $\beta_{kj}\neq 0$) with a small change of $e$ given a value $(w,d)$.
After choosing an appropriate logit${}^{-1}d$ according to the prior
assumption about model sparsity, the hyperparameter $e$ is optimized by a grid
search strategy suggested by [28]. If one assumes the model sparsity lies in
the range $(c_{1},c_{2})$, then one can specify the hyperparameter $d$ subject
to logit${}^{-1}d=c_{1}$ and search hyperparameter $e$ from 0 to
$-\frac{dmp}{2c_{2}}(\mathbbm{1}^{\top}G\mathbbm{1})^{-1}$ (see Supplementary
S1 for more details). This can be tailored to avoid searching in the space of
too dense models which would be computationally slow.

###  2.3 SSUR model with MRF prior and random effects

The methods in Section 2.1 and 2.2 both assume independent samples. However,
the large-scale drug screens may have cell line samples from different cancer
tissue types. We use random effects to address the heterogeneity of multiple
sample groups.

Let $\mathbf{Z}_{n\times T}$ be indicator variables representing the $T$
heterogeneous sample groups. We define a multivariate linear regression which
specifies a spike-and-slab prior (4), hyper-inverse Wishart prior (5), MRF
prior (8) and random effects:

$\displaystyle\begin{split}\mathbf{Y}&=\mathbf{Z}\bm{B}_{0}+\mathbf{X}\bm{B}+\mathbf{U},\\\
\beta_{0,tj}|w_{0}&\sim\mathcal{N}(0,\ w_{0}),\\\
\beta_{kj}|\gamma_{kj},w&\sim\gamma_{kj}\mathcal{N}(0,\
w)+(1-\gamma_{kj})\delta_{0}(\beta_{kj}),\\\
w_{0}&\sim\mathcal{IG}(a_{w_{0}},b_{w_{0}}),\\\
w&\sim\mathcal{IG}(a_{w},b_{w}),\\\
\bm{\gamma}|d,e,G&\propto\exp\\{d\mathbbm{1}^{\top}\bm{\gamma}+e\bm{\gamma}^{\top}G\bm{\gamma}\\},\\\
\text{vec}\\{\mathbf{U}\\}&\sim\mathcal{N}(\mathbf{0},\
\Psi\otimes\mathbb{I}_{n}),\\\
\Psi&\sim\mathcal{HIW}_{\mathcal{G}}(\nu,\tau\mathbb{I}_{m}),\\\
\tau&\sim\mathcal{G}amma(a_{\tau},b_{\tau}),\\\ \end{split}$ (9)

where the random effects
$\bm{B}_{0}=\\{\beta_{0,tj}:t=1,\cdots,T;j=1,\cdots,m\\}$, $\bm{B}_{0}$,
$\bm{B}$ and $\mathbf{U}$ are mutually independent.

Let us look into details of the random effects. For the $i$th subject and the
$j$th response variable, we have

$y_{ij}=x_{i\cdot}^{\top}\bm{\beta}_{j}+z_{i\cdot}^{\top}\bm{\beta}_{0,j}+u_{ij}.$

Then the covariance between the $j$th and $j^{\prime}$th responses is

$\mathbb{C}\text{ov}[y_{ij},y_{ij^{\prime}}|u_{\cdot\cdot}]=\mathbb{C}\text{ov}[x_{i\cdot}^{\top}\bm{\beta}_{j}+z_{i\cdot}^{\top}\bm{\beta}_{0,j},x_{i\cdot}^{\top}\bm{\beta}_{j^{\prime}}+z_{i\cdot}^{\top}\bm{\beta}_{0,j^{\prime}}]=0.$

Although the priors on the coefficients $\bm{B}_{0}$ and $\bm{B}$ don’t
specify any correlations between different responses for the same subject, the
inverse Wishart or hyper-Inverse Wishart prior on $\Psi$
($=\mathbb{C}\text{ov}[u_{i\cdot}]$) models their variances. After a similar
reparametrization to (3), the relationships between response variables reflect
the transformed parameter $\bm{\rho}$.

The covariance between the $i$th and $i^{\prime}$th subjects for the $j$th
response is

$\displaystyle\mathbb{C}\text{ov}[y_{ij},y_{i^{\prime}j}|u_{\cdot\cdot}]$
$\displaystyle=\mathbb{C}\text{ov}[x_{i\cdot}^{\top}\bm{\beta}_{j}+z_{i\cdot}^{\top}\bm{\beta}_{0,j},x_{i^{\prime}\cdot}^{\top}\bm{\beta}_{j}+z_{i^{\prime}\cdot}^{\top}\bm{\beta}_{0,j}]$
$\displaystyle=wx_{i\cdot}^{\top}x_{i^{\prime}\cdot}+w_{0}z_{i\cdot}^{\top}z_{i^{\prime}\cdot}$
$\displaystyle=\begin{cases}wx_{i\cdot}^{\top}x_{i^{\prime}\cdot},&\text{if
}i\text{th}\text{ and }i^{\prime}\text{th subjects belong to different
groups},\\\ wx_{i\cdot}^{\top}x_{i^{\prime}\cdot}+w_{0},&\text{if
}i\text{th}\text{ and }i^{\prime}\text{th subjects belong to the same
group}.\end{cases}$

Prior specification

Here the hyperparameter $w_{0}$ determines the difference of covariances
between subjects within the same group or not. In the pharmacogenomic datasets
from multiple cancer tissues, the tissue effects to predict drug sensitivity
are usually stronger than the genomic effects. Thus it is better to specify
larger hyperparameter $w_{0}$ than $w$. Similar to specify the hyperparameter
$a_{w}$, $a_{w_{0}}=const_{a_{0}}\cdot 1/2\cdot m\cdot T$. $const_{a_{0}}$ and
$b_{w_{0}}$ are chosen to satisfy the 0.95-quantile of the posterior
$\mathcal{IG}(a_{w_{0}}+1/2mT,b_{w_{0}})$ larger than
$1.96\mathbb{V}ar[\beta_{0,tj}]=0$, where the variance of the random effect
can be from previous studies or expert knowledge about typical effect size in
similar studies.

###  2.4 Posterior computation

Similar to [3], the joint posterior distribution of model (2) with inverse
Wishart and MRF priors is

$\displaystyle
f(\bm{B},\mathbf{\Gamma},w,\bm{\rho},\bm{\sigma}^{2},\tau|\mathbf{Y},\mathbf{X})$
$\displaystyle=$ $\displaystyle
f(\mathbf{Y}|\mathbf{X},\bm{B},\bm{\rho},\bm{\sigma}^{2})f(\bm{B}|\mathbf{\Gamma},w)f(\mathbf{\Gamma}|G,d,e)f(w)f(\bm{\rho}|\bm{\sigma}^{2},\tau)f(\bm{\sigma}^{2}|\tau)f(\tau)$
(9) $\displaystyle=$
$\displaystyle\prod_{j}f(\bm{y}_{j}|\mathbf{X},\bm{B},\bm{\rho},\bm{\sigma}^{2})\prod_{k,j}f(\beta_{kj}|\gamma_{kj},w)f(\bm{\gamma}|G,d,e)f(w)\prod_{j,l<j}f(\rho_{jl}|\sigma_{j}^{2},\tau)\prod_{j}f(\sigma_{j}^{2}|\tau)f(\tau),$

where $\bm{\rho}$ and $\bm{\sigma}^{2}$ are vectors of $\\{\rho_{jl}\\}$ and
$\\{\sigma_{j}^{2}\\}$, respectively. Since
$\bm{y}_{j}|\mathbf{X},\bm{B},\bm{\rho},\bm{\sigma}^{2}$ is a normal
distribution with mean
$\mathbf{X}\bm{\beta}_{j}+\sum_{l<j}\bm{u}_{l}\rho_{jl}$ and variance
$\sigma_{j}^{2}\mathbb{I}_{n}$, it is not complicated to obtain the full
conditionals of regression coefficients, $w$, $\sigma_{j}^{2}$, $\rho_{jl}$
and $\tau$. The posterior of the indicator variable $\mathcal{\Gamma}$ can be
estimated by the Metropolis-Hastings sampler.

The inference of the SUR model with hyper-inverse Wishart and MRF priors is
similar to (9), but each cumulative product in (9) with respect to the
response index $j$ or $l$ will be multiple products with respect to the clique
index of graph $\mathcal{G}$ and the response index of each clique. If there
are random effects for group factors as (9), then the joint posterior
distribution (9) will include parameters $\bm{B}_{0}$ and $w_{0}$, i.e.,

$\displaystyle
f(\bm{B}_{0},\bm{B},\mathbf{\Gamma},w_{0},w,\bm{\rho},\bm{\sigma}^{2},\tau|\mathbf{Y},\mathbf{X},\mathbf{Z})$
$\displaystyle=$ $\displaystyle
f(\mathbf{Y}|\mathbf{X},\mathbf{Z},\bm{B}_{0},\bm{B},\bm{\rho},\bm{\sigma}^{2})f(\bm{B}_{0}|w_{0})f(w_{0})$
$\displaystyle\times
f(\bm{B}|\mathbf{\Gamma},w)f(\mathbf{\Gamma}|G,d,e)f(w)f(\bm{\rho}|\bm{\sigma}^{2},\tau)f(\bm{\sigma}^{2}|\tau)f(\tau)$

We implement Gibbs samplers to obtain posterior estimates of $\bm{B}$,
$\bm{\rho}$ and $\bm{\sigma}^{2}$, and update the latent indicator
$\bm{\Gamma}$ via a Metropolis-Hastings sampler with parallel tempering the
same as [3]. The Thompson sampling [42] is used as the proposal for each
latent indicator $\gamma_{kj}$. The hyperparameter $\tau$ is updated via a
random walk Metropolis sampler as [3] proposed.

To overcome the prohibitive computational time in high-dimensional settings,
the ESS algorithm [9, 39] is used to update the posteriors. For each iteration
of MCMC, after sampling the latent indicator variables $\bm{\Gamma}$, we first
update most other hyperparameters, and then update the reparametrized
parameters $\bm{\sigma}^{2}$ and $\bm{\rho}$, and finally update the
regression coefficients.

  1. (i)

sampling latent indicator variables $\bm{\Gamma}$ using the Metropolis-
Hastings sampler;

  2. (ii)

sampling hyperparameter $\tau$ using the random walk Metropolis sampler;

  3. (iii)

sampling hyperparameter $w$ (and $w_{0}$) using the Gibbs sampler;

  4. (iv)

sampling the graph $\mathcal{G}$ from the junction tree sampler (essentially
Metropolis-Hastings sampling);

  5. (v)

sampling reparametrized $\bm{\sigma}^{2}$ and parameters $\bm{\rho}$ from
their posteriors (Gibbs sampler);

  6. (vi)

sampling coefficients $\bm{B}$ (and $\bm{B}_{0}$) from their posteriors (Gibbs
sampler).

At each iteration, the ESS algorithm implements a local move to add/delete and
swap the latent indicators within each chain, and then a global move is to
exchange and crossover the latent indicators between any two parallel tempered
chains. The temperature across all response variables is adapted based on the
acceptance rate of the global exchange operator. The ESS algorithm with
parallel tempering is potentially effective to search a high-dimensional model
space with multiple modes [9].

###  2.5 Model performance evaluation

To evaluate the performance of the proposed methods, the structure recovery
and prediction performance are considered. The structures recovery includes
the estimation of the latent indicator matrix $\mathbf{\Gamma}$ which captures
the relationships between the multiple response variables and high-dimensional
predictors, and the estimation of the graph $\mathcal{G}$ which addresses the
relationships between response variables. The prediction accuracy of the
Bayesian models can be estimated based on posterior predictive distribution by
the leave-one-out cross-validation (LOO) and widely applicable information
criterion (WAIC) [46], and based on the median probability model (MPM) [4] by
the mean squared error (MSE) and mean squared prediction error (MSPE).

[46] proposed an efficient computation for the Bayesian LOO estimate of out-
of-sample predictive fit

$\text{elpd}_{\text{loo}}=\sum_{j=1}^{m}\sum_{i=1}^{n}\log
f(y_{ij}|\bm{y}_{(-i)j})\approx\sum_{j=1}^{m}\sum_{i=1}^{n}\frac{1}{\frac{1}{N}\sum_{t=1}^{N}\frac{1}{f(y_{ij}|\bm{\vartheta}^{(t)})}},$

where $\bm{y}_{(-i)j}$ is the observation vector of the $j$th response
variable except the $i$th observation, $N$ is the length of an MCMC chain and
$\bm{\vartheta}^{(t)}$ is the estimates of the $t$th iteration for all related
parameters. The WAIC is estimated by

$\widehat{\text{elpd}}_{\text{waic}}=\widehat{\text{elpd}}_{\text{loo}}-\sum_{j=1}^{m}\sum_{i=1}^{n}\mathbb{V}\text{ar}_{t=1}^{N}[\log
f(y_{ij}|\bm{\vartheta}^{(t)})],$

where the second term above is used as a measure of the complexity of the
model.

For future prediction, a single model may be required that can be used by
practitioners or evaluated in different situations. [4] suggested the median
probability model (MPM) which is often the optimal predictive model. So we
have the MPM for each coefficient $\beta_{kj}$:

$\hat{\beta}_{kj,MPM}=\begin{cases}\frac{\sum_{t=1}^{N}\beta_{kj}^{(t)}}{\sum_{t=1}^{N}\gamma_{kj}^{(t)}},&\text{if
}\frac{\sum_{t=1}^{N}\gamma_{kj}^{(t)}}{N}>0.5,\\\
0,&\text{otherwise},\end{cases}$

where $\gamma_{kj}^{(t)}$ is the estimate of the $t$th MCMC iteration for the
latent indicator variable of $\beta_{kj}$. After obtaining
$\hat{\bm{B}}_{MPM}=\\{\hat{\beta}_{kj,MPM}\\}$, the MSE and MSPE are

MSE
$\displaystyle=\frac{1}{mn}\|\mathbf{Y}-\mathbf{X}\hat{\bm{B}}_{MPM}\|^{2},$
MSPE
$\displaystyle=\frac{1}{mn^{\prime}}\|\mathbf{Y}^{*}-\mathbf{X}^{*}\hat{\bm{B}}_{MPM}\|^{2},$

where $\mathbf{Y}_{n\times m}$ and $\mathbf{X}_{n\times p}$ were used to
estimate $\hat{\bm{B}}_{MPM}$, and $\mathbf{Y}_{n^{\prime}\times m}^{*}$ and
$\mathbf{X}_{n^{\prime}\times p}^{*}$ are new data.

##  3 Simulation study

In this section, we evaluate the structure recovery and prediction performance
of the proposed methods described in Section 2 on two simulation settings, one
with independent samples for the SSUR model with MRF prior and the other with
grouped or heterogeneous samples for the SSUR model with MRF prior and random
effects. The first method will be compared with the SSUR model with hotspot
prior.

###  3.1 Simulation scenario

To construct a complex structure between the response variables and
predictors, a network is designed as as shown in Figure 2(a). It assumes three
groups of the 20 response variables and six groups of 120 relevant predictors
among 300 total predictors. The first group responses (i.e.,
$\bm{y}_{1},\cdots,\bm{y}_{5}$) are related to four groups predictors, i.e.,
$\bm{x}_{1}\sim\bm{x}_{5}$, $\bm{x}_{30}\sim\bm{x}_{50}$,
$\bm{x}_{51}\sim\bm{x}_{60}$ and $\bm{x}_{110}\sim\bm{x}_{120}$. The second
group responses (i.e., $\bm{y}_{6},\cdots,\bm{y}_{12}$) are also related to
four groups predictors, i.e., $\bm{x}_{10}\sim\bm{x}_{20}$,
$\bm{x}_{51}\sim\bm{x}_{60}$, $\bm{x}_{70}\sim\bm{x}_{90}$ and
$\bm{x}_{110}\sim\bm{x}_{120}$. The third group responses (i.e.,
$\bm{y}_{13},\cdots,\bm{y}_{20}$) are related to three groups predictors,
i.e., $\bm{x}_{30}\sim\bm{x}_{50}$, $\bm{x}_{370}\sim\bm{x}_{90}$ and
$\bm{x}_{110}\sim\bm{x}_{120}$. Corresponding to this network (Figure 2(a)), a
sparse latent indicator matrix $\bm{\Gamma}$ (Figure 2(b)) reflects the same
relationships between the response variables and predictors in the
multivariate regression model (1). In addition, a decomposable graph
$\mathcal{G}$ (Figure 2(c)) is added to reflect the extra relationships
between the response variables. The graph $\mathcal{G}$ has six blocks which
denote six possible groups among the 20 response variables. That means all
given 300 predictors are not able to explain all relationships between
response variables. In the SUR model (1), the information $\mathcal{G}$ is
included in residuals and expected to be recovered by our methods.

Figure 2:  True relationships within and between response variables and
predictors. (a) Network between multiple response variables and predictors
(left panel); (b) sparse $\bm{\Gamma}$ for the associations between the
response variables and predictors corresponding to (a); (c) additional
structure $\mathcal{G}$ between response variables besides the three groups
structure in (a). Black indicates a true relation between or within response
variables and predictors.

All data are generated based on a multivariate linear regression model

$\mathbf{Y}=\mathbbm{1}\bm{\alpha}^{\top}+(\mathbf{X}\bm{B})_{\mathbf{\Gamma}}+\mathbf{U}.$

The input data $\mathbf{X}=\\{x_{ik}\\}$ ($i=1,\cdots,250;k=1,\cdots,300$) are
simulated independently from the standard normal distribution
$\mathcal{N}(0,1)$. The intercepts $\bm{\alpha}=\\{\alpha_{j}\\}$ and
coefficients $\bm{B}=\\{\beta_{kj}\\}$ ($k=1,\cdots,300;j=1,\cdots,20$) are
simulated independently from the normal distribution
$\mathcal{N}(0.5,0.2^{2})$ but truncated by the latent indicator variable
$\mathbf{\Gamma}=\\{\gamma_{kj}\\}$. The noise matrix $\mathbf{U}$ is
simulated based on a standard multivariate normal distribution and a G-Wishart
distribution [33]. Algorithm 1 shows the detailed steps for the simulation
scenario.

To assess the proposed SSUR model with MRF prior and random effects, Algorithm
2 is designed which is similar to Algorithm 1. Algorithm 2 includes group
indicators $\mathbf{Z}=\\{z_{it}\\}$ ($i=1,\cdots,250;t=1,\cdots,4$) with
$(z_{i1},z_{i2},z_{i3},z_{i4})$ independently from a multinomial distribution
$mulitinomial(0.1,0.2,0.3,0.4)$ and random effects
$\bm{B}_{0}=\\{\beta_{0,tj}\\}$ ($t=1,\cdots,4;j=1,\cdots,20$) independently
from a normal distribution $\mathcal{N}(0,2^{2})$. The response data are
generated from the linear mixed model

$\mathbf{Y}=(\mathbf{X}\bm{B})_{\mathbf{\Gamma}}+\mathbf{Z}\bm{B}_{0}+\mathbf{U}.$

We run a MCMC sampler with 5 chains, 500 000 iterations and a burn-in of 300
000 iterations. Both simulation algorithms generate validation datasets
independently with the same sample size to evaluate the performance of the
proposed methods.

Algorithm 1 Simulation steps without random effects

1: Design a decomposable graph $\mathcal{G}$, dim$(\mathcal{G})=m\times m$

2: Design a sparse matrix $\bm{\Gamma}$, dim$(\Gamma)=p\times m$

3: Simulate $x_{ik}\sim\mathcal{N}(0,1)$, $i=1,\cdots,n$ and $k=1,\cdots,p$

4: Simulate $\alpha_{j}$, $\beta_{kj}\sim\mathcal{N}(0.5,0.2^{2})$,
$k=1,\cdots,p$ and $j=1,\cdots,m$

5: Simulate $\tilde{u}_{ij}\sim\mathcal{N}(0,1)$, $i=1,\cdots,n$ and
$j=1,\cdots,m$

6: Simulate $P\sim\mathcal{W}_{\mathcal{G}}(3,M)$ (G-Wishart distribution)
where diagonals of $M$ are 1 and off-diagonals are 0.9, dim$(P)=m\times m$

7: Use Cholesky decomposition $\text{chol}(P^{-1})$ to get
$\mathbf{U}=\tilde{\mathbf{U}}\cdot\text{chol}(P^{-1})$

8: Generate
$\mathbf{Y}=\mathbbm{1}\bm{\alpha}^{\top}+(\mathbf{X}\bm{B})_{\mathbf{\Gamma}}+\mathbf{U}$

Algorithm 2 Simulation steps with random effects

Design a decomposable graph $\mathcal{G}$, dim$(\mathcal{G})=m\times m$

2: Design a sparse matrix $\bm{\Gamma}$, dim$(\Gamma)=p\times m$

Simulate $x_{ik}\sim\mathcal{N}(0,1)$, $i=1,\cdots,n$ and $k=1,\cdots,p$

4: Simulate $\beta_{kj}\sim\mathcal{N}(0.5,0.2^{2})$, $k=1,\cdots,p$ and
$j=1,\cdots,m$

Simulate $\bm{z}_{i}\sim multinomial(0.1,0.2,0.3,0.4)$, $i=1,\cdots,n$

6: Simulate $\beta_{0,tj}\sim\mathcal{N}(0,2^{2})$, $t=1,\cdots,T$ and
$j=1,\cdots,m$

Simulate $\tilde{u}_{ij}\sim\mathcal{N}(0,1)$, $i=1,\cdots,n$ and
$j=1,\cdots,m$

8: Simulate $P\sim\mathcal{W}_{\mathcal{G}}(3,M)$ (G-Wishart distribution)
where diagonals of $M$ are 1 and off-diagonals are 0.9, dim$(P)=m\times m$

Use Cholesky decomposition $\text{chol}(P^{-1})$ to get
$\mathbf{U}=\tilde{\mathbf{U}}\cdot\text{chol}(P^{-1})$

10: Generate
$\mathbf{Y}=(\mathbf{X}\bm{B})_{\mathbf{\Gamma}}+\mathbf{ZB}_{0}+\mathbf{U}$

###  3.2 Comparison of SSUR models with hotspot and MRF priors

We first compare our approach, the SSUR model with MRF prior, to the SSUR
model with hotspot prior based on the simulation Algorithm 1. For the
convenience of comparison, the simulated $\mathbf{Y}$ and $\mathbf{X}$ are
standardized before fitting the two models. Our approach uses the network in
Figure 2 as prior information for the MRF prior, which is implemented as the
illustration in Section 2.2. Throughout this article, we refer to a predictor
as being selected or identified, if the corresponding latent indicator
variable has posterior mean larger than 0.5. Figure 3 shows that the SSUR
model with hotspot prior and our approach both have good recovery for the
coefficients matrix $\bm{\Gamma}$ and relationship of response variables
(i.e., $\mathcal{G}$). Table 1 displays the accuracy of variable selection and
prediction performance. The SSUR model with MRF prior has slightly higher
accuracy and sensitivity for variable selection than the SSUR model with
hotspot prior, i.e., slightly better identification of truly associated
predictors. The better variable selection is also more possible to lead the
better MSPE (0.271) based on validation data.

Figure 3:  Estimators by the SSUR model with hotspot prior (panels (a) and (b)) and the SSUR model with MRF prior (panels (c) and (d))  Table 1:  Accuracy of variable selection and prediction performance by the SSUR models with hotspot priors and MRF prior  prior for $\mathbf{\Gamma}$ | accuracy | sensitivity | specificity | MSE | MSPE
---|---|---|---|---|---
hotspot | 0.985 | 0.930 | 0.996 | 0.101 | 0.309
MRF | 0.989 | 0.996 | 0.988 | 0.102 | 0.271

###  3.3 Sensitivity analysis of the SSUR model with MRF prior

Since the MRF prior may be strongly informative, sensitivity analysis is done
by using different $G$ for the MRF prior. Different prior information of the
network is provided by deleting partial true edge potentials uniformly or
blockly or adding noise edge potentials111Deleting edge potentials uniformly,
for example deleting 1% means that the
$1+(1-1/|\mathcal{E}|)/1\%\cdot\\{0:(1\%\cdot|\mathcal{E}|)\\}$th edges of the
whole true edge list are deleted. Deleting 1% blocks of edges means that the
last 1% blocks of edges are deleted. The true edge list includes the edges of
each pattern (i.e., association block) together. Adding 1% noise edges means
that there are $1\%\cdot mp(mp-1)/2$ wrong edges are included randomly.. The
three cases are as follows.

  * •

Case 1: delete 1%, 10%, 50% and 90% edges uniformly from the full informative
$G$, which always remain partial true associations across all patterns in
$\bm{\Gamma}$.

  * •

Case 2: delete 1%, 10%, 50%, 90% and 100% blocks of edges from the informative
$G$ which gradually provide less until no association knowledge in
$\bm{\Gamma}$.

  * •

Case 3: add 0.1%, 0.5%, 1%2221% noise edges are already more than the total
true edges which are only $\sim$ 0.3% of all possible edges. noise edges into
the full informative $G$.

Table 2 shows that our approach can identify truly associated predictors very
well even after deleting true edges uniformly (i.e., case 1) or adding noise
edges (i.e., case 2) in the MRF prior. This reveals that even without knowing
90% true edges, almost all truly associated predictors for the corresponding
response variables can still be identified if given 10% true edges across most
patterns. It means that if a little information about the true associated
predictors is given, more true associated predictors that are related to the
given predictors can be identified. Case 3 reveals that if the provided
information is wrong, our approach can still perform good variable selection.
However, in case 2, the variable selection becomes slightly worse without
providing complete patterns (i.e., deleting blocks of edges). Figure 4
indicates the information of deleted blocks cannot be recovered fully. But
even without any prior knowledge on the associations (i.e., deleting 100%
blocks of edges in case 2) a sparser $\mathbf{\Gamma}$ is estimated, which
includes most true patterns. Supplementary S2 shows that all cases have good
structure recovery of the relationships between response variables (i.e.,
$\hat{\mathcal{G}}$).

In terms of prediction performance, our approach has almost the same MSE and
MSPE after deleting true edges uniformly or adding noise edges in the MRF
prior (i.e., case 1). Although deleting 90% blocks of edges in case 2, our
approach still has quite good MSE (0.260). However, the predictions based on
validation data become worse after deleting more blocks of edges in case 2,
i.e., worse MSPE. Comparing with the SSUR model with hotspot prior in Table
tab:paramEst, the prediction performance of our approach is still comparable
without prior knowledge of the 90% true associated predictors (i.e., deleting
90% blocks of edges in case 2).

Table 2:  Sensitivity analysis of the SSUR model with different MRF priors

accuracy sensitivity specificity
$\widehat{\text{elpd}_{\text{loo}}}^{{\dagger}}$
$\widehat{\text{elpd}_{\text{waic}}}^{{\dagger}}$ MSE MSPE delete edges
uniformly 1% 0.989 0.995 0.988 -6685 -6733 0.101 0.267 10% 0.989 0.997 0.988
-6688 -6744 0.106 0.274 50% 0.988 0.994 0.987 -6691 -6734 0.108 0.275 90%
0.990 0.995 0.989 -6710 -6756 0.110 0.283 delete blocks of edgess 1% 0.989
0.996 0.988 -6667 -6729 0.101 0.271 10% 0.987 0.976 0.989 -6678 -6739 0.156
0.343 50% 0.976 0.905 0.990 -6654 -6718 0.246 0.473 90%‡ 0.977 0.883 0.996
-6658 -6717 0.260 0.477 100%‡ 0.967 0.803 0.999 -6689 -6743 0.408 0.611 add
noise edges 0.1% 0.989 0.997 0.988 -6689 -6744 0.106 0.273 0.5% 0.989 0.993
0.988 -6699 -6745 0.101 0.273 1% 0.989 0.997 0.988 -6714 -6762 0.102 0.266
†The comparison of the elpd values should be considered with caution, since
the number of edges influences the posterior predictive directly and further
changes the elpd. ‡The hyperparameters in the MRF prior are chosen more
properly for deleting 90% blocks of edges ($d=-0.1$, $e=5$) and for deleting
100% blocks of edges ($d=-0.5$). All others cases have $d=-2$ and $e=0.1$.

Figure 4:  Sensitivity analysis case 2: delete blocks of edges

###  3.4 Results and discussion of the SSUR model with MRF prior and random
effects

In Algorithm 2, $T=4$ indicator variables that reflect sample groups are
included to assess the performance of the SSUR model with MRF prior and random
effects in Section 2.3. As described in Section2.1 and 2.3, the
hyperparameters of the shrinkage coefficients and random effects can be
specified appropriately for the model without data standardization. Here we
choose $const_{a_{w_{0}}}=3/2$, $b_{w_{0}}=300$, $const_{a_{w}}=1/2$,
$const_{b_{w}}=2$ and the model sparsity $r_{sparse}=20\%$ to obtain
$a_{w_{0}}=60$, $a_{w}=300$ and $b_{w}=300$. Figure 5 and Table 3 both show
good recovery of the true associated predictors $\hat{\mathbf{\Gamma}}$ and
relationships between response variables. The estimated random effects of
group variables have larger error (i.e.,
$\frac{1}{\sqrt{mT}}\|\hat{\bm{B}}_{0}-\bm{B}_{0}\|_{\ell_{2}}^{2}$) than the
estimated effects of other predictors (i.e.,
$\frac{1}{\sqrt{mp}}\|\hat{\bm{B}}_{MPM}-\bm{B}\|_{\ell_{2}}^{2}$). This is
because each true random effect is simulated from a normal distribution with
larger standard deviation 2 than other effects with standard deviation 0.2.

Figure 5:  Posterior mean of $\mathbf{\Gamma}$ and $\mathcal{G}$ by the SSUR model with MRF prior and random effects based on the simulated data from Algorithm 2 Table 3:  Structure recovery and prediction by the SSUR model with MRF prior and random effects  | accuracy | sensitivity | specificity | $\frac{1}{\sqrt{mT}}\|\hat{\bm{B}}_{0}-\bm{B}_{0}\|_{\ell_{2}}$ | $\frac{1}{\sqrt{mp}}\|\hat{\bm{B}}_{MPM}-\bm{B}\|_{\ell_{2}}$
---|---|---|---|---|---
$\mathcal{G}$ | 0.935 | 0.871 | 1.000 | - | -
Coefficients | 0.988 | 0.992 | 0.988 | 1.315 | 0.514

##  4 Applications

To validate the methodology that we proposed in practice, a subset of the
Genomics of Drug Sensitivity in Cancer (GDSC) dataset [48] is used to identify
the relationship between multiple cancer drugs and high-dimensional
predictors. The preprocessed data sets by [15] are from
ftp://ftp.sanger.ac.uk/pub4/cancerrxgene/releases/release-5.0/. We selected
three groups of cancer drugs: four MEK inhibitors (i.e., RDEA119, PD-0325901,
CI-1040 and AZD6244), two Bcr-Abl tyrosine kinase inhibitors (i.e., Nilotinib
and Axitinib) and one generic chemotherapy agent (i.e., Methotrexate). The
seven drugs were tested on 498 selected cell lines from 13 cancer tissue
types. The drug sensitivity of the cell lines was summarized by the log(IC50)
from in vitro drug does response experiments. In order to (initially) explore
the relationship between the three groups of drugs and their genomic
information, we preselected 343 gene expression features (continuous and log-
transformed), 426 copy number variation features (counts) and 68 mutated genes
(binary).

Figure 6: Illustration of the relationship between drugs and a group of
related genes which is used for the MRF prior. The left panel is for the Bcr-
Abl fusion gene and the corresponding related genes. The right panel is for
all drugs and gene TP35 as one example with features representing all three
data sources. The rectangles indicate drugs, solid circles indicate genes and
dashed circles indicate that elements inside are related. The edges between
drugs and dashed circled genes indicate that those genes corresponding to the
multiple drugs or individual drugs are related. The names with suffix “.GEX”,
“.CNV” and “.MUT” are features of expression, copy number variation and
mutation, respectively.

Figure 6 illustrates the prior information used to construct the graph for the
MRF prior. The four MEK inhibitors were developed in order to reduce the
activity of the MAPK/ERK pathway, and the two Bcr-Abl tyrosine kinase
inhibitors were developed in order to inhibit Bcr-Abl tyrosine kinase activity
and proliferation of Bcr-Abl expressing cells. The left panel of Figure 6
illustrates different characteristics of BCR and ABL genes corresponding to
drugs Nilotinib and Axitinib are correlated, so the model is likely to select
these genes together. In addition, different characteristics of the same gene
are correlated and are likely to be selected simultaneously as well.

Based on Figure 6, we construct edge potentials for the MRF prior in the
proposed model (9) in Section 2.3. To specify the hyperparameters for the
random effects for the 13 cancer tissue types. We assume each cancer tissue
tissue with prior mean 0 and prior standard deviation 3, which results in the
hyperparameters $w_{0}\sim\mathcal{IG}(54.6,400)$ according to the prior
specification in Section 2.3. Similarly, we assume each genomic predictor with
prior mean -0.1 and prior standard deviation 2 and crude proportion of
associated genomic predictors as 0.1, which results in the hyperparameters
$w\sim\mathcal{IG}(11.718,70.308)$ according to the prior specification in
Section 2.1. The two hyperparameters in the MRF prior are chosen as $d=-2.7$
and $e=0.1$ according to the prior specification in Section 2.2. We run the
MCMC algorithm with 25 chains, 800 000 iterations and a burn-in of 500 000
iterations.

The SSUR model with MRF prior and random effects have good MCMC diagnostics
(see Supplementary S4 for more details). Figure 7 shows the estimated network
representation between the seven drugs and three omics characteristics. There
are in total 72 of 837 features selected. Drugs Methotrexate, Nilotinib,
Axitinib, RDEA119, PD-0325901, CI-1040 and AZD6244 are found to be associated
with 13, 15, 11, 60, 57, 57 and 59 features, respectively.

Figure 7: Estimated network between the seven drugs and selected genes based
on thresholds 0.5.

Drugs RDEA119, PD-0325901, CI-1040 and AZD6244 are clustered together in
Figure 7 which confirms our prior expectations, since these drugs share the
same target pathway, i.e., targeting the MAPK/ERK pathway. In fact, genes
FGFR1, FGFR3, FGFR3.CNV, FGFR3.MUT, PDGFRA, PDGFRA.CNV, PDGFRA.MUT, PDGFRB,
PDGFRB.CNV, AKT2 and MAP2K4 which are important genes in the MAPK/ERK pathway
are all identified here.

The two Bcr-Abl tyrosine kinase inhibitors are clustered in Figure 7. Only
drug nilotinib identifies BCR-ABL mutation with a strong effect (i.e.,
regression coefficient=-2.854 in Figure 9). Based on the “Drug Target Common
2.0” platform [45], three genes (i.e., ALK.MUT, EXT1, PIM1) of the 14 selected
features for nilotinib are off-targets; four genes (i.e., ALK.MUT, MST1R,
PDGFRA, PIM1) of the 11 selected features for axitinib are off-targets.

Drug methotrexate is not related to all other drugs except axitinib in Figure
7. It is a generic chemotherapy agent and does not target a specific molecular
pathway. Therefore we do not expect to find strong effects of molecular
features. Although there are 13 features selected for methotrexate, they have
small effects ($<0.6$).

Figure 8: MPM estimator for coefficients. The marked red colour gene names are
target genes of the MEK inhibitors and the green ones are also identified for
the MEK inhibitors. The marked blue colour gene name is the target of the Bcr-
Abl inhibitors and the light blue colour gene names are also identified for
the Bcr-Abl inhibitors.

Looking at all estimated effects in Figure 8, including random effects for
cancer tissue types and regression coefficients for genes, we see that cancer
tissue types have stronger associations with drug-response than genes in
general. This is consistent with the results in [51]. Methotrexate has a
strong effect (-2.601) on blood cancer tissue relative to other tissue types
(Figure 9) and small negative effects on the thyroid (-0.446), lung (-0.083)
and bone (-0.294) cancer tissues. That means the cell lines of those cancer
tissues have decreasing log(IC50) when treating by methotrexate, since a
smaller IC50 value implies a faster response, i.e., that the cell line is
sensitive to the treatment with this drug. It is known to be an effective
chemotherapeutic agent for a number of cancers, including leukemia [38],
lymphoma [35], lung cancer [11] and osteosarcoma [23].

Figure 9: Posterior means of cancer tissue random effects for all drugs.
Random effects are centred around zero.

In Figure 9, Drugs RDEA119 and PD-0325901 both have negative effects on
digestive, thyroid and soft (e.g., rhabdomyosarcoma) cancer tissue types. [22,
32] showed RDEA119 as a highly selective MEK inhibitor for colon and thyroid
cancers respectively. Drug PD-0325901 also shows a weak negative effect
(-0.232) on blood tissue but a strong negative effect (-1.772) skin cancer
tissue type which was also revealed by [22]. It also has a negative effect
(-0.508) on the pancreas cancer tissue type. [21] showed that PD-0325901
decreased the pancreatic growth in mice. Also, a weak effect (-0.039) on
aerodigestive cancer tissue type. But [34] found its potential anti-tumor
effects in human head and neck squamous cell carcinoma. Drug CI-1040 has a
strong negative effect (-1.340) on skin cancer tissue but weak negative effect
(-0.321) on thyroid cancer tissue type. [40] assessed the selective CI-1040
(PD184352) for melanoma. However, drug AZD6244 only has weak negative effects
on the skin (-0.300) and soft (-0.147) cancer tissue types, although it was
discovered to treat non-small cell lung and thyroid cancers.

The blood tissue effect (-0.492) of drug nilotinib validates its known effect
on some blood cancers, in particular chronic myeloid leukemia [43]. Although
it has a negative thyroid tissue effect (-0.982), there is still no biological
or clinical evidence. Axitinib only has weak negative effects on urogenital
(-0.006), blood (-0.120) and bone (-0.185) cancer tissue types which were also
shown by [25, 18, 26].

##  5 Discussion

In this work, we have developed a structured multivariate Bayesian variable
selection model to predict multiple related response variables based on high-
dimensional structured predictors and identify important predictors (i.e.,
variable selection). With our model we want to (a) be able to borrow
information within and between responses and predictors, (b) learn the network
between responses and predictors, and (c) understand the residual covariance.
The proposed approach allows us to make use of partial network information
(i.e., on the relationships between responses and/or predictors) in the prior
for the selection indicator $\mathbf{\Gamma}$ and to further simultaneously
select sparse predictors and learn structures between response variables. In
addition, we have allowed for the accounting of sample heterogeneity through
random effects, which are excluded from the variable selection. Guidance for
specifying all hyperparameters in our model has also been provided.

Through the simulation studies, we have demonstrated that this approach can
recover the network structure between multiple response variables and
predictors and improve prediction accuracy even if given a little prior
knowledge but across most patterns (Figure 4(d)). Even without prior knowledge
of the network (i.e., $G=\bm{0}$), the SSUR model with MRF prior (with proper
hyper-parameters) is still able to recover the general patterns in
$\mathcal{G}$ and $\bm{\Gamma}$ (Figure 4(e)), but the estimated model might
be sparser than the truth. In this situation, the MRF prior degenerates to a
Bernoulli prior. Our approach is also robust to the noise information (i.e.,
false patterns) in the MRF prior (Table 2). More sensitivity analysis about
the hyper-inverse Wishart prior have been done in Supplementary S1, which
shows the SSUR model with a good MRF prior (i.e., good prior knowledge for the
relationship between responses and predictors) can reach good variable
selection and prediction performance even if the hyper-inverse Wishart prior
is misspecified. The SSUR model with hotspot prior is slightly sensitive to
the degree of freedom of the hyper-inverse Wishart prior than our approach.

The real data analysis has shown that our approach is able to identify true
targeted genes and certain off-targets for cancer drugs. The estimated
relationships between cancer drugs also validate the drug groups with common
targeted genes or pathways. The estimated cancer tissue effects can explain
the corresponding drugs’ effectiveness. However, drugs might also be effective
even if cancer tissues have non-negative effects, as long as the drugs’ target
genes have negative effects. Although we have addressed the heterogeneity of
different cancer tissue types through random effects, the estimated random
effect might be the total effect of multiple cancer subtypes within one
tissue. Alternatives can be to extend our model to multilevel random effects
or employ a flexible Dirichlet process for the random effect [31, 19].

In terms of computational efficiency, the junction tree sampler is efficient
to learn the relationships between response variables, although the assumption
of a decomposable graph might limit the searching graph space. [14] and [3]
show certain theoretical property and practical gains by the decomposable
graph. Our approach also inherits the ESS algorithm which might be time-
consuming if there are too many related response variables. An alternative is
variational inference that approximates posterior densities for Bayesian
models [7, 50]. For nonlinear models, a novel and efficient MCMC algorithm was
proposed by [1], which can explore the high-dimensional model space and
estimate the structure among multiple responses of diverse types (i.e., low-
intensity counts, binary, ordinal and continuous variables).

## Supplementary Materials

Software in the form of R package BayesSUR is available from the Comprehensive
R Archive Network (CRAN) at https://CRAN.R-project.org/package=BayesSUR and
all R codes for the simulated studies and real data analysis are available on
https://github.com/zhizuio/BayesSUR-RE. The attached supplementary materials
include the derivation of the expected upper bound of the hyperparameter $e$
in the MRF prior (Supplementary S1) and all estimated graph
$\hat{\mathcal{G}}$ for the sensitivity analysis in Section 3.3 (Supplementary
S2). More sensitivity analysis of the hyper-inverse Wishart prior in the SSUR
model with MRF prior based on simulation Algorithm 1 can also be found in the
Supplementary S3. The MCMC diagnostics are in Supplementary S4.

## Acknowledgements

#A. Lewin and M. Zucknick are joint last authors. The first author is
supported by the Faculty of Medicine at the University of Oslo, NORBIS and the
project BIG INSIGHT (237718) funded by the Research Council of Norway The
authors thank Dr. Leonardo Bottolo for discussions. Conflict of Interest: None
declared.

## References

  * [1] Angelos Alexopoulos and Leonardo Bottolo “Bayesian Variable Selection for Gaussian Copula Regression Models” In _Journal of Computational and Graphical Statistics_ 0.0 Taylor & Francis, 2020, pp. 1–16 DOI: 10.1080/10618600.2020.1840997
  * [2] Tomohiro Ando and Arnold Zellner “Hierarchical Bayesian analysis of the seemingly unrelated regression and simultaneous equations models using a combination of direct Monte Carlo and importance sampling techniques” In _Bayesian Analysis_ 5.1 International Society for Bayesian Analysis, 2010, pp. 65–95 DOI: 10.1214/10-BA503
  * [3] Marco Banterle et al. “Sparse variable and covariance selection for high-dimensional seemingly unrelated Bayesian regression” In _bioRxiv_ Cold Spring Harbor Laboratory, 2018, pp. 467019 DOI: 10.1101/467019
  * [4] Maria Maddalena Barbieri and James O. Berger “Optimal predictive model selection” In _The Annals of Statistics_ 32.3 The Institute of Mathematical Statistics, 2004, pp. 870–897 DOI: 10.1214/009053604000000238
  * [5] J. Barretina et al. “The Cancer Cell Line Encyclopedia enables predictive modelling of anticancer drug sensitivity” In _Nature_ 483.7391 Springer Nature Publishing AG, 2012, pp. 603–607 DOI: 10.1038/nature11003
  * [6] Anindya Bhadra and B.K. Mallick “Joint High-Dimensional Bayesian Variable and Covariance Selection with an Application to eQTL Analysis” In _Biometrics_ 69.2 The International Biometric Society, 2013, pp. 447–457 DOI: 10.1111/biom.12021
  * [7] David M. Blei, Alp Kucukelbir and Jon D. McAuliffe “Variational Inference: A Review for Statisticians” In _Journal of the American Statistical Association_ 112.518 Taylor & Francis, 2017, pp. 859–877 DOI: 10.1080/01621459.2017.1285773
  * [8] Leonardo Bottolo et al. “Bayesian detection of expression quantitative trait loci hot-spots” In _Genetics_ 189.4 Genetics, 2011, pp. 1449–1459 DOI: 10.1534/genetics.111.131425
  * [9] Leonardo Bottolo and Sylvia Richardson “Evolutionary Stochastic Search for Bayesian Model Exploration” In _Bayesian Analysis_ 5.3 International Society for Bayesian Analysis, 2010, pp. 583–618 DOI: 10.1214/10-BA523
  * [10] Carlos M. Carvalho, Helene Massam and Mike West “Simulation of Hyper-Inverse Wishart Distributions in Graphical Models” In _Biometrika_ 94.3 Oxford University Press, Biometrika Trust, 2007, pp. 647–659 DOI: 10.1093/biomet/asm056
  * [11] Yu-Hung Chen et al. “Methotrexate Conjugated to Gold Nanoparticles Inhibits Tumor Growth in a Syngeneic Lung Tumor Model” In _Molecular Pharmaceutics_ 4.5, 2007, pp. 713–722 DOI: 10.1021/mp060132k
  * [12] A Philip Dawid “Some matrix-variate distribution theory: notational considerations and a Bayesian application” In _Biometrika_ 68.1 Oxford University Press, Biometrika Trust, 1981, pp. 265–274 DOI: 10.1093/biomet/68.1.265
  * [13] Sameer K Deshpande, Veronika Ročková and Edward I George “Simultaneous variable and covariance selection with the multivariate spike-and-slab lasso” In _Journal of Computational and Graphical Statistics_ 28.4 Taylor & Francis, 2019, pp. 921–931
  * [14] A.M. Fitch, M.B. Jones and H. Massam “The performance of covariance selection methods that consider decomposable models only” In _Bayesian Analysis_ 9.3 International Society for Bayesian Analysis, 2014, pp. 659–684 DOI: 10.1214/14-BA874
  * [15] MJ. Garnett et al. “Systematic identification of genomic markers of drug sensitivity in cancer cells” In _Nature_ 483.7391 Springer Nature Publishing AG, 2012, pp. 570–575 DOI: 10.1038/nature11005
  * [16] Joe W. Gray and Gordon B. Mills “Large-Scale Drug Screens Support Precision Medicine” In _Cancer Discovery_ 5.11 American Association for Cancer Research, 2015, pp. 1130–1132 DOI: 10.1158/2159-8290.CD-15-0945
  * [17] Peter J Green and Alun Thomas “Sampling decomposable graphs using a Markov chain on junction trees” In _Biometrika_ 100.1 Oxford University Press, Biometrika Trust, 2013, pp. 91–110
  * [18] S Halbach et al. “Axitinib and sorafenib are potent in tyrosine kinase inhibitor resistant chronic myeloid leukemia cells” In _Cell Communication and Signaling_ 14 Springer Nature, 2016, pp. 6 DOI: 10.1186/s12964-016-0129-y
  * [19] Felix Heinzl, Ludwig Fahrmeir and Thomas Kneib “Additive mixed models with Dirichlet process mixture and P-spline priors” In _AStA Advances in Statistical Analysis_ 96, 2012, pp. 47–68 DOI: 10.1007/s10182-011-0161-6
  * [20] C.C. Holmes, D.G.T. Denison and B.K. Mallick “Accounting for Model Uncertainty in Seemingly Unrelated Regressionsa” In _Journal of Computational and Graphical Statistics_ 11.3 Taylor & Francis, Ltd., 2002, pp. 533–551 DOI: 10.1198/106186002475
  * [21] Bryan J. Holtz et al. “ERK activation is required for CCK-mediated pancreatic adaptive growth in mice” In _American Journal of Physiology-Gastrointestinal and Liver Physiology_ 307.7, 2014, pp. G700–G710 DOI: 10.1152/ajpgi.00163.2014
  * [22] Cory Iverson et al. “RDEA119/BAY 869766: A Potent, Selective, Allosteric Inhibitor of MEK1/2 for the Treatment of Cancer” In _Cancer Research_ 69.17 American Association for Cancer Research, 2009, pp. 6839–6847 DOI: 10.1158/0008-5472.CAN-09-0679
  * [23] Norman Jaffe “Recent advances in the chemotherapy of metastatic osteogenic sarcoma” In _Cancer_ 30.6, 1972, pp. 1627–1631
  * [24] Zhenyu Jia and Shizhong Xu “Mapping Quantitative Trait Loci for Expression Abundance” In _Genetics_ 176.1 Genetics, 2007, pp. 611–623 DOI: 10.1534/genetics.106.065599
  * [25] Jose A. Karama et al. “Phase 2 Trial of Neoadjuvant Axitinib in Patients with Locally Advanced Nonmetastatic Clear Cell Renal Cell Carcinoma” In _European Urology_ 66.5, 2014, pp. 874–880 DOI: 10.1016/j.eururo.2014.01.035
  * [26] Lauren T. Kerr, Jacqueline F. Donoghue, Alexander L. Wilding and Terrance G. Johns “Axitinib has antiangiogenic and antitumorigenic activity in myxoid liposarcoma” In _International Journal of Cancer_ 2016, 2016, pp. 3484673 DOI: 10.1155/2016/3484673
  * [27] Eiru Kim et al. “A network of human functional gene interactions from knockout fitness screens in cancer cells” In _Life Science Alliance_ 2.2 Life Science Alliance, 2019 DOI: 10.26508/lsa.201800278
  * [28] Kyu Ha Lee et al. “Multivariate Bayesian variable selection exploiting dependence structure among outcomes: Application to air pollution effects on DNA methylation” In _Biometrics_ 73.1 Genetics, 2017, pp. 232–241 DOI: 10.1111/biom.12557
  * [29] Alex Lewin et al. “MT-HESS: an efficient Bayesian approach for simultaneous association detection in OMICS datasets, with application to eQTL mapping in multiple tissues” In _Bioinformatics_ 32.4 Oxford University Press, 2016, pp. 523–532
  * [30] Fan Li and Nancy R. Zhang “Bayesian Variable Selection in Structured High-Dimensional Covariate Spaces With Applications in Genomics” In _Journal of the American Statistical Association_ 105.491 Taylor & Francis, 2010, pp. 1202–1214 DOI: 10.1198/jasa.2010.tm08177
  * [31] Yisheng Li, Xihong Lin and Peter Müller “Bayesian Inference in Semiparametric Mixed Models for Longitudinal Data” In _Biometrics_ 66.1, 2010, pp. 70–78 DOI: 10.1111/j.1541-0420.2009.01227.x
  * [32] Dingxie Liu, Joanna Xing, Barry Trink and Mingzhao Xing “BRAF mutation-selective inhibition of thyroid cancer cells by the novel MEK inhibitor RDEA119 and genetic-potentiated synergism with the mTOR inhibitor temsirolimus” In _International Journal of Cancer_ 127.12, 2010, pp. 2965–2973 DOI: 10.1002/ijc.25304
  * [33] Reza Mohammadi and Ernst Wit “BDgraph: An R package for Bayesian structure learning in graphical models” In _Journal of Statistical Software_ 89.3 Foundation for Open Access Statistics, 2019, pp. 1–30 DOI: 10.18637/jss.v089.i03
  * [34] Suresh Mohan et al. “MEK Inhibitor PD-0325901 Overcomes Resistance to PI3K/mTOR Inhibitor PF-5212384 and Potentiates Antitumor Effects in Human Head and Neck Squamous Cell Carcinoma” In _Clinical Cancer Research_ 21.17, 2015, pp. 3946–3956 DOI: 10.1158/1078-0432.CCR-14-3377
  * [35] Patrick G. Morris et al. “Rituximab, Methotrexate, Procarbazine, and Vincristine Followed by Consolidation Reduced-Dose Whole-Brain Radiotherapy and Cytarabine in Newly Diagnosed Primary CNS Lymphoma: Final Results and Long-Term Outcome” In _Journal of Clinical Oncology_ 31.31, 2013, pp. 3971–3979 DOI: 10.1200/JCO.2013.50.4910
  * [36] Tea Pemovska et al. “Individualized Systems Medicine Strategy to Tailor Treatments for Patients with Chemorefractory Acute Myeloid Leukemia” In _Cancer Discovery_ 3.12 American Association for Cancer Research, 2013, pp. 1416–1429 DOI: 10.1158/2159-8290.CD-13-0350
  * [37] Enrico Petretto et al. “New Insights into the Genetic Control of Gene Expression using a Bayesian Multi-tissue Approach” In _PLoS Computational Biology_ 6.4, 2010, pp. e1000737 DOI: 10.1371/journal.pcbi.1000737
  * [38] Bayard L. Powell et al. “Arsenic trioxide improves event-free and overall survival for adults with acute promyelocytic leukemia: North American Leukemia Intergroup Study C9710” In _Blood_ 116.19, 2010, pp. 3751–3757 DOI: 10.1182/blood-2010-02-269621
  * [39] Sylvia Richardson, Leonardo Bottolo and Jeffrey S. Rosenthal “Bayesian Models for Sparse Regression Analysis of High Dimensional Data” In _Bayesian Statistics 9_ Oxford Scholarship Online, 2011, pp. 539–568 DOI: 10.1093/acprof:oso/9780199694587.001.0001
  * [40] Victoria Roulstone et al. “BRAF- and MEK-Targeted Small Molecule Inhibitors Exert Enhanced Antimelanoma Effects in Combination With Oncolytic Reovirus Through ER Stress” In _Molecular Therapy_ 23.5 Elsevier Inc., 2014, pp. 931–942 DOI: 10.1038/mt.2015.15
  * [41] Hélène Ruffieux “Large-scale variational inference for Bayesian joint regression modelling of high-dimensional genetic data” In _PhD dissertation_ , 2019
  * [42] Daniel J. Russo et al. “A Tutorial on Thompson Sampling” In _Foundations and Trends in Machine Learning_ 11.1, 2018, pp. 1–96 DOI: 10.1561/2200000070
  * [43] Tomasz Sacha and Giuseppe Saglio “Nilotinib in the treatment of chronic myeloid leukemia” In _Future Oncology_ 15.9, 2019, pp. 953–965 DOI: 10.2217/fon-2018-0468
  * [44] Francesco C. Stingo, Yian A. Chen, Mahlet G. Tadesse and Marina Vannucci “Incorporating biological information into linear models: A Bayesian approach to the selection of pathways and genes” In _Ann. Appl. Stat._ 5.3 The Institute of Mathematical Statistics, 2011, pp. 1978–2002 DOI: 10.1214/11-AOAS463
  * [45] Jing Tang et al. “Drug target commons: a community effort to build a consensus knowledge base for drug-target interactions” In _Cell Chemical Biology_ 25.2, 2018, pp. 224–229 DOI: 10.1016/j.chembiol.2017.11.009
  * [46] Aki Vehtari, Andrew Gelman and Jonah Gabry “Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC” In _Statistics and Computing_ 27.5 Springer Nature Publishing AG, 2017, pp. 1413–1432 DOI: 10.1007/s11222-016-9696-4
  * [47] Hao Wang “Sparse Seemingly Unrelated Regression Modelling: Applications in Finance and Econometrics” In _Comput. Stat. Data Anal._ 54.11 Elsevier Science Publishers B. V., 2010, pp. 2866–2877 DOI: 10.1016/j.csda.2010.03.028
  * [48] W Yang et al. “Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells” In _Nucleic Acids Reserch_ 41.Database issue Oxford University Press, 2013, pp. D955–61 DOI: 10.1093/nar/gks1111
  * [49] Arnold Zellner and Tomohiro Ando “A direct Monte Carlo approach for Bayesian analysis of the seemingly unrelated regression model” In _Journal of Econometrics_ 159.1 Elsevier B.V., 2010, pp. 33–45 DOI: 10.1016/j.jeconom.2010.04.005
  * [50] C. Zhang, J. Butepage, H. Kjellstrom and S. Mandt “Advances in Variational Inference” In _IEEE Transactions on Pattern Analysis & Machine Intelligence_ 41.8 IEEE Computer Society, 2019, pp. 2008–2026 DOI: 10.1109/TPAMI.2018.2889774
  * [51] Zhi Zhao and Manuela Zucknick “Structured penalized regression for drug sensitivity prediction” In _Journal of the Royal Statistical Society, Series C (Applied Statistics)_ 69.3 Royal Statistical Society, Wiley, 2020, pp. 525–545

Supplementary Materials

### S1 Expected upper bound for $e$ in the MRF prior

Comparing the MRF prior (8) with a $mp$-dimensional multivariate normal
density

$(2\pi)^{-\frac{mp}{2}}|\bm{Q}|^{\frac{1}{2}}\exp\left\\{-\frac{1}{2}(\bm{\gamma}-\bm{\mu})^{\top}\bm{Q}(\bm{\gamma}-\bm{\mu})\right\\}\propto\exp\left\\{\bm{\mu}^{\top}\bm{Q}\bm{\gamma}-\frac{1}{2}\bm{\gamma}^{\top}\bm{Q}\bm{\gamma}\right\\},$

then let

$\displaystyle\bm{\mu}^{\top}\bm{Q}$ $\displaystyle=d\mathbbm{1}^{\top},$
$\displaystyle-\frac{1}{2}\bm{Q}$ $\displaystyle=eG.$

Thus we have

$-2e\bm{\mu}^{\top}G=d\mathbbm{1}^{\top}.$ (*)

If assuming the expected model sparsity $c_{2}$, we have

$\bm{\mu}=\mathbb{E}[\bm{\gamma}]=c_{2}\mathbbm{1}.$ (**)

Based on (* ‣ S1 Expected upper bound for $e$ in the MRF prior) and (** ‣ S1
Expected upper bound for $e$ in the MRF prior), we obtain

$e=-\frac{dmp}{2c_{2}}(\mathbbm{1}^{\top}G\mathbbm{1})^{-1}.$

### S2 Estimated $\mathcal{G}$ for the sensitivity analysis in Section 3.3

Figure S2.1: Sensitivity analysis case 1: delete edges uniformly

Figure S2.2: Sensitivity analysis case 2: delete blocks of edges Figure S2.3:
Sensitivity analysis case 3: add noise edges

### S3 More sensitivity analysis for the hyper-inverse Wishart prior

To evaluate the specification of the hyper-inverse Wishart prior for the SUR
model with a MRF prior, we show the results with different $\nu$, $\tau$ and
$\eta$ by simulating data according to Algorithm 1. Although a very large
$\nu=10^{4}$ enforces unrelated responses ($\hat{\mathcal{G}}$ in Figure S3.1,
the variable selection $\hat{\mathbf{\Gamma}}$ and prediction still perform
well ($\hat{\mathbf{\Gamma}}$ in FigureS3.1 and Table S3.1 if using the median
probability model. However, the posterior mean of $\mathbf{\Gamma}$ of the
SSUR model with hotspot prior ($\hat{\mathcal{G}}$ in Figure S3.5 is more
sensitive than the SSUR model with MRF prior. Figure S3.1-S3.4 show that other
hyperparameters of the hyper-inverse WIshart prior are not sensitive to the
structure recovery (i.e., $\hat{\mathbf{\Gamma}}$ and $\hat{\mathcal{G}}$) and
prediction performance (i.e., accuracy, MSE and MSPE) . So is the SSUR model
with hotspot prior (Figure S3.5-S3.8).

Figure S3.1: Posterior mean of the coefficients matrix $\hat{\bm{B}}$, latent indicator variable $\hat{\mathbf{\Gamma}}$ and learning structure $\hat{\mathcal{G}}$ by the SSUR model with MRF prior, $\mathcal{HIW}_{\mathcal{G}}({\bf 10^{4}},\tau\mathbbm{I})$, $\tau\sim\mathcal{G}amma(0.1,10)$ and $\eta\sim\mathcal{B}eta(0.1,1)$. Figure S3.2: Posterior mean of the coefficients matrix $\hat{\bm{B}}$, latent indicator variable $\hat{\mathbf{\Gamma}}$ and learning structure $\hat{\mathcal{G}}$ by the SSUR model with MRF prior, $\mathcal{HIW}_{\mathcal{G}}(22,\tau\mathbbm{I})$, $\tau\sim\mathcal{G}amma({\bf 10^{-6}},10)$ and $\eta\sim\mathcal{B}eta(0.1,1)$. Figure S3.3: Posterior mean of the coefficients matrix $\hat{\bm{B}}$, latent indicator variable $\hat{\mathbf{\Gamma}}$ and learning structure $\hat{\mathcal{G}}$ by the SSUR model with MRF prior, $\mathcal{HIW}_{\mathcal{G}}(22,\tau\mathbbm{I})$, $\tau\sim\mathcal{G}amma({\bf 10^{6}},10)$ and $\eta\sim\mathcal{B}eta(0.1,1)$. Figure S3.4: Posterior mean of the coefficients matrix $\hat{\bm{B}}$, latent indicator variable $\hat{\mathbf{\Gamma}}$ and learning structure $\hat{\mathcal{G}}$ by the SSUR model with MRF prior, $\mathcal{HIW}_{\mathcal{G}}(22,\tau\mathbbm{I})$, $\tau\sim\mathcal{G}amma(0.1,10)$ and $\eta\sim\mathcal{B}eta({\bf 10^{-6}},1)$. Figure S3.5: Posterior mean of the coefficients matrix $\hat{\bm{B}}$, latent indicator variable $\hat{\mathbf{\Gamma}}$ and learning structure $\hat{\mathcal{G}}$ by the SSUR model with hotspot prior, $\mathcal{HIW}_{\mathcal{G}}({\bf 10^{4}},\tau\mathbbm{I})$, $\tau\sim\mathcal{G}amma(0.1,10)$ and $\eta\sim\mathcal{B}eta(0.1,1)$. Figure S3.6: Posterior mean of the coefficients matrix $\hat{\bm{B}}$, latent indicator variable $\hat{\mathbf{\Gamma}}$ and learning structure $\hat{\mathcal{G}}$ by the SSUR model with hotspot prior, $\mathcal{HIW}_{\mathcal{G}}(22,\tau\mathbbm{I})$, $\tau\sim\mathcal{G}amma({\bf 10^{-6}},10)$ and $\eta\sim\mathcal{B}eta(0.1,1)$. Figure S3.7: Posterior mean of the coefficients matrix $\hat{\bm{B}}$, latent indicator variable $\hat{\mathbf{\Gamma}}$ and learning structure $\hat{\mathcal{G}}$ by the SSUR model with hotspot prior, $\mathcal{HIW}_{\mathcal{G}}(22,\tau\mathbbm{I})$, $\tau\sim\mathcal{G}amma({\bf 10^{6}},10)$ and $\eta\sim\mathcal{B}eta(0.1,1)$. Figure S3.8: Posterior mean of the coefficients matrix $\hat{\bm{B}}$, latent indicator variable $\hat{\mathbf{\Gamma}}$ and learning structure $\hat{\mathcal{G}}$ by the SSUR model with hotspot prior, $\mathcal{HIW}_{\mathcal{G}}(22,\tau\mathbbm{I})$, $\tau\sim\mathcal{G}amma(0.1,10)$ and $\eta\sim\mathcal{B}eta({\bf 10^{-6}},1)$. Table S3.1: Accuracy of variable selection and prediction performance by the SSUR models with MRF/hotspot prior and different hyper-inverse Wishart prior | accuracy | sensitivity | specificity | $\widehat{\text{elpd}_{\text{loo}}}$ | $\widehat{\text{elpd}_{\text{waic}}}$ | MSE | MSPE
---|---|---|---|---|---|---|---
SSUR model with MRF prior |  |  |  |  |  |  |
$\nu=22$, $\tau\sim\mathcal{G}amma(0.1,10)$, $\eta\sim\mathcal{B}eta(0.1,1)^{{\dagger}}$ | 0.989 | 0.996 | 0.988 | -6698 | -6752 | 0.102 | 0.271
$\nu=10^{4}$, $\tau\sim\mathcal{G}amma(0.1,10)$, $\eta\sim\mathcal{B}eta(0.1,1)$ | 0.978 | 0.988 | 0.976 | -19587 | -19587 | 0.064 | 0.271
$\nu=22$, $\tau\sim\mathcal{G}amma(10^{-6},10)$, $\eta\sim\mathcal{B}eta(0.1,1)$ | 0.989 | 0.997 | 0.988 | -6687 | -6739 | 0.098 | 0.265
$\nu=22$, $\tau\sim\mathcal{G}amma(10^{6},10)$, $\eta\sim\mathcal{B}eta(0.1,1)$ | 0.989 | 0.997 | 0.988 | -6681 | -6725 | 0.100 | 0.265
$\nu=22$, $\tau\sim\mathcal{G}amma(0.1,10)$, $\eta\sim\mathcal{B}eta(10^{-6},1)$ | 0.989 | 0.993 | 0.988 | -6705 | -6765 | 0.103 | 0.272
SSUR model with hotspot prior |  |  |  |  |  |  |
$\nu=22$, $\tau\sim\mathcal{G}amma(0.1,10)$, $\eta\sim\mathcal{B}eta(0.1,1)^{{\dagger}}$ | 0.985 | 0.930 | 0.995 | -3078 | -3078 | 0.101 | 0.309
$\nu=10^{4}$, $\tau\sim\mathcal{G}amma(0.1,10)$, $\eta\sim\mathcal{B}eta(0.1,1)$ | 0.970 | 0.926 | 0.979 | -19420 | -19420 | 0.063 | 0.282
$\nu=22$, $\tau\sim\mathcal{G}amma(10^{-6},10)$, $\eta\sim\mathcal{B}eta(0.1,1)$ | 0.984 | 0.933 | 0.994 | -3123 | -3123 | 0.103 | 0.312
$\nu=22$, $\tau\sim\mathcal{G}amma(10^{6},10)$, $\eta\sim\mathcal{B}eta(0.1,1)$ | 0.986 | 0.922 | 0.999 | -3076 | -3076 | 0.111 | 0.315
$\nu=22$, $\tau\sim\mathcal{G}amma(0.1,10)$, $\eta\sim\mathcal{B}eta(10^{-6},1)$ | 0.986 | 0.925 | 0.999 | -3057 | -3057 | 0.107 | 0.302
†These specifications are the same as the SSUR models in Table 1. |  |

### S4 Sensitivity analysis of the simulated data with different G-Wishart
distributions

Table S4: Accuracy of variable selection and prediction performance by the SSUR models with MRF prior and different G-Wishart distribution $M$ in G-Wishart distribution | accuracy | sensitivity | specificity | MSE | MSPE
---|---|---|---|---|---
off-diagnals 0.2 | 0.988 | 0.990 | 0.987 | 0.098 | 0.265
off-diagnals 0.5 | 0.988 | 0.989 | 0.987 | 0.097 | 0.265
Figure S4.1: Posterior mean of the coefficients matrix $\hat{\bm{B}}$, latent
indicator variable $\hat{\mathbf{\Gamma}}$ and learning structure
$\hat{\mathcal{G}}$ by the SSUR model with hotspot prior based on simulation
Algorithm 1 with off-diagonals 0.2 of $M$ in G-Wishart distribution. Figure
S4.2: Posterior mean of the coefficients matrix $\hat{\bm{B}}$, latent
indicator variable $\hat{\mathbf{\Gamma}}$ and learning structure
$\hat{\mathcal{G}}$ by the SSUR model with hotspot prior based on simulation
Algorithm 1 with off-diagonals 0.5 of $M$ in G-Wishart distribution.

### S5 Sensitivity analysis of the simulated data with error
$\bm{u_{i}}\sim\mathcal{N}(\bm{0},\Sigma_{u})$

To avoid tailoring to the approach with hyper-inverse Wishart prior, we
simulate correlated errors directly from two different multivariate normal
distributions $\bm{u_{i}}\sim\mathcal{N}(\bm{0},\Sigma_{u})$, where

$\Sigma_{u}=\begin{pmatrix}A_{u}&\bm{0}&\bm{0}\\\ \bm{0}&B_{u}&\bm{0}\\\
\bm{0}&\bm{0}&\mathbbm{I}_{8}\end{pmatrix},$

$A_{u}$ is 8-by-8 matrix with diagonals 1 and off-diagonals $u_{a}$, $B_{u}$
is 4-by-4 matrix with diagonals 1 and off-diagonals $u_{b}$. The first case
(i.e., case S6.1) has $(u_{a},u_{b})=(0.1,0.2)$. The second case (i.e., case
S6.2) has $(u_{a},u_{b})=(0.5,0.6)$. Table S6 shows that our approach has good
variable selection and prediction performance in both cases. However, our
approach cannot estimated weak correlations (i.e., $(u_{a},u_{b})=(0.1,0.2)$)
between response variables.

Table S5: Accuracy of variable selection and prediction performance by the SSUR models with MRF prior and different correlated $\bm{u_{i}}$ $(u_{a},u_{b})$ | accuracy | sensitivity | specificity | MSE | MSPE
---|---|---|---|---|---
(0.1,0.2) | 0.984 | 0.992 | 0.982 | 0.056 | 0.258
(0.5,0.6) | 0.983 | 0.991 | 0.981 | 0.051 | 0.262
Figure S5.1: Posterior mean of the coefficients matrix $\hat{\bm{B}}$, latent
indicator variable $\hat{\mathbf{\Gamma}}$ and learning structure
$\hat{\mathcal{G}}$ by the SSUR model with hotspot prior based on simulation
Algorithm 1 with $\bm{u_{i}}\sim\mathcal{N}(\bm{0},\Sigma_{u})$ in case S6.1.
Figure S5.2: Posterior mean of the coefficients matrix $\hat{\bm{B}}$, latent
indicator variable $\hat{\mathbf{\Gamma}}$ and learning structure
$\hat{\mathcal{G}}$ by the SSUR model with hotspot prior based on simulation
Algorithm 1 with $\bm{u_{i}}\sim\mathcal{N}(\bm{0},\Sigma_{u})$ in case S6.2.

### S6 MCMC diagnostics of the GDSC data analysis

Figure S6: Diagnostic plots of the MCMC sampler.