arXiv:1304.7284v2 [cs.LG] 16 Oct 2013

Supervised Heterogeneous Multiview Learning
for Joint Association Study and Disease
Diagnosis
Shandian Zhe
Purdue University
szhe@purdue.edu

Zenglin Xu
Purdue University
xu218@purdue.edu

Yuan Qi
Purdue University
alanqi@cs.purdue.edu
October 17, 2013
Abstract
Given genetic variations and various phenotypical traits, such as Magnetic Resonance Imaging (MRI) features, we consider two important and related tasks in
biomedical research: i)to select genetic and phenotypical markers for disease diagnosis and ii) to identify associations between genetic and phenotypical data.
These two tasks are tightly coupled because underlying associations between genetic variations and phenotypical features contain the biological basis for a disease. While a variety of sparse models have been applied for disease diagnosis and
canonical correlation analysis and its extensions have bee widely used in association studies (e.g., eQTL analysis), these two tasks have been treated separately.
To unify these two tasks, we present a new sparse Bayesian approach for joint association study and disease diagnosis. In this approach, common latent features
are extracted from different data sources based on sparse projection matrices and
used to predict multiple disease severity levels based on Gaussian process ordinal
regression; in return, the disease status is used to guide the discovery of relationships between the data sources. The sparse projection matrices not only reveal interactions between data sources but also select groups of biomarkers related to the
disease. To learn the model from data, we develop an efficient variational expectation maximization algorithm. Simulation results demonstrate that our approach
achieves higher accuracy in both predicting ordinal labels and discovering associations between data sources than alternative methods. We apply our approach to an
imaging genetics dataset for the study of Alzheimer’s Disease (AD). Our method
identifies biologically meaningful relationships between genetic variations, MRI
features, and AD status, and achieves significantly higher accuracy for predicting
ordinal AD stages than the competing methods.

1

1

Introduction

Recent advances in biomedical research have provided new opportunities to study diseases – for example, Alzheimer’s disease (AD), the most common neurodegenerative
disorder – from multiple data sources. For example, one data source contains genetic
variations, such as single nucleotide polymorphisms (SNPs), which can help us understand the genetic basis of diseases. Another data source can be molecular and clinical phenotypes, such as Magnetic Resonance Imaging (MRI) data, which can reveal
important phenotypic changes in patients.Finding associations between different data
sources can reveal unknown biological relationships and has a wide range of applications in computational biology [1], epidemiology [2], computational neural science [3],
and imaging genetics [4]. In addition to the genotypes and phenotypic traits, we have
valuable labeled information about disease stages from patient medical records. Thus
we face a new data analysis setting where the objective is two-fold: i) finding associations between different data sources and ii) selecting relevant (groups of) features from
all the sources to predict ordinal disease stages.
Many statistical approaches have been developed to discover associations or select
features (or variables) for prediction in a high dimensional problem. For association
studies, representative approaches are canonical correlation analysis (CCA) and its extensions [5, 6]. These approaches treat different data sources as separate linear projections from a common latent representation. These approaches have been widely used
in expression quantitative trait locus (eQTL) analysis. For example, Parkhomenko et
al. [7] applied sparse CCA (sCCA) to find relationships between genetic loci and gene
expression levels in Utah families; Witten and Tibshirani [8] used sCCA to reveal associations between gene expression and DNA copy variation; and Chen et al. [9] used
structured CCA for pathway selection. For disease diagnosis based on high dimensional biomarkers, popular approaches include lasso [10], elastic net [11], and group
lasso [12], and Bayesian automatic relevance determination [13, 14]. Here we treat
genotypes or phenotypes as predictors (i.e., biomarkers) and the disease status as the
response in a linear regression or classification setting. Non-zero estimated regression
or classification weights indicate relevant biomarkers for the disease [15, 16].
Despite their wide success in many applications, these approaches are limited by
the following factors:
• Most association studies neglect the supervision from the disease status. Because
many diseases, such as AD, are a direct result of genetic variations and often
highly correlated to clinical traits, the disease status provides useful yet currently
unutilized information for finding relationships between genetic variations and
clinical traits.
• For disease diagnosis, most sparse approaches use classification models and do
not consider the order of disease severity. For subjects in AD studies, there is
a natural severity order from being normal to mild cognitive impairment (MCI)
and then from MCI to AD. Classification models cannot capture the order in
AD’s severity levels. Furthermore, the classification approaches are often based
on conditional models (e.g., logistic regression) and ignore relationships between
multiple views.
2

• Most previous methods are not designed to handle heterogeneous data types. The
SNPs values are discrete (and ordinal based on an additive genetic model), while
the imaging features are continuous. Popular CCA or lasso-type methods simply
treat both of them as continuous data and overlook the heterogeneous nature of
the data.
To address these problems, we propose a new Bayesian approach that unifies multiview learning with sparse ordinal regression for joint association study and disease
diagnosis. In the new approach, genetic variations and phenotypical traits are generated from common latent features based on separate sparse projection matrices and
suitable link functions and the common latent features are used to predict the disease
status based on Gaussian process ordinal regression (See Section 2). To enforce sparsity in projection matrices, we assign spike and slab priors [17] over them; these priors
have been shown to be more effective than l1 penalty to learn sparse projection matrices
[18, 19]. The sparse projection matrices not only reveal critical interactions between
the different data sources but also identify groups of biomarkers in data relevant to disease status. Finding groups of biomarkers can avoid over-sparsification (i.e., selecting
one instead of multiple correlated features), thus boosting the accuracy for disease diagnosis. It can also help provide a better biological understanding because these groups
may form biologically units (i.e., pathways). Meanwhile, via its direct connection to
the latent features, the disease status influences the estimation of the projection matrices
so that it can guide the discovery of associations between heterogeneous data sources
relevant to the disease. Hence we name this new method Supervised Heterogeneous
Multiview Learning (SHML).
To learn the model from data, we develop a variational Bayesian expectation maximization (VB-EM) approach (See Section 3). It iteratively minimizes the Kullback
Leibler divergence between a tractable approximation and exact Bayesian posterior
distributions and provides an estimate to the model marginal likelihood. Maximizing
this estimate enables us to automatically choose a suitable dimension for the latent
features in a principled Bayesian framework.
In Section 4, we test our approach SHML on both synthetic and real datasets. On
synthetic data, SHML achieves both higher estimation accuracy in recovering true associations between different views than CCA and sparse CCA, and higher prediction
accuracy than multiple advanced alternative methods, such as the combination of CCA
and elastic net, and Gaussian process ordinal regression [20]. We then apply SHML
to an AD study. AD accounts for 60-80% of age-related dementia cases – one in eight
older Americans has AD – and there is no cure for AD till now. It is believed that its
underlying pathology precedes the onset of cognitive symptoms for many years [21].
Although AD studies have attracted a lot of attention from both academia and industry [22, 23], to our best knowledge, our paper presents the first (supervised) study to
uncover associations between genotypes and phenotypic traits relevant to AD. Our results on Alzheimer’s Disease Neuroimaging Initiative (ADNI) data show that SHML
achieves highest prediction accuracy among all the competing methods. Furthermore,
SHML finds biologically meaningful predictive relationships between SNPs, MRI features, and AD status.

3

y

Sh

H

U

G

Z

Sg

η

X

Figure 1: The graphical model of Supervised Heterogeneous Multiview Learning, where X is
the continuous view, Z is the ordinal view, and y are the labels.

2

Model

First, let us describe the data. We assume there are two heterogeneous data sources: one
contains continuous data – for example, MRI features – and one discrete ordinal data –
for instance, SNPs. Note that we can easily generalize our model below to handle more
views and other data types by adopting suitable link functions (e.g., a Possion model for
count data). Given data from n subjects, p continuous features and q discrete features,
we denote the continuous data by a p × n matrix X = [x1 , . . . , xn ], the discrete ordinal
data by a q × n matrix Z = [z1 , . . . , zn ], and the labels (i.e., the disease status) by a
n × 1 vector y = [y1 , . . . , yn ]> . For the AD study, we let yi = 0, 1, and 2 if the i-th
subject is in the normal, MCI or AD condition, respectively.
To link two data sources X and Z together, we introduce common latent features
U = [u1 , . . . , un ] and assume X and Z are generated from U by sparse projection.
The common latent feature assumption is sensible for association studies because both
SNPs and MRI features are biological measurements of the same subjects. Note that
ui is the latent feature for the i-th subject and its dimension k is estimated by evidence
maximization.
In a Bayesian framework, we give a Gaussian prior over U, p(U) =
Q
N
(u
|0,
I),
and specify the rest of the model (see Figure 1) as follows:
i
i
• Continuous data. Given U, X is generated from
p(X|U, G, η) =

n
Y

N (xi |Gui , η −1 I)

i=1

where G = [g1 , g2 , ...gp ]> is a p × k projection matrix, I is an identity matrix, and η −1 I is the precision matrix of the Gaussian distribution. We assign a
Gamma prior over η, p(η|r1 , r2 ) = Gamma(η|r1 , r2 ) where r1 and r2 are the
hyperparameters and set to be 10−3 in our experiments.
• Ordinal data. For an ordinal observation z ∈ {0, 1, . . . ,
R − 1} its value is decided by which region an auxiliary variable c falls in
−∞ = b0 < b1 < . . . < bR = ∞.
If c falls in [br , br+1 ), z is set to be r. For the AD study, the SNPs Z takes
values in {0, 1, 2} and therefore R = 3. Given a q × k projection matrix H =
[h1 , h2 , ...hq ]> , the auxiliary variables C = {cij } and the ordinal data Z are
4

generated from
p(Z, C|U, H) =

q Y
n
Y

p(cij |hi , uj )p(zij |cij )

i=1 j=1

where
p(cij |hi , uj ) = N (cij |h>
i uj , 1)
p(zij |cij ) =

2
X

δ(zij = r)δ(br ≤ cij < br+1 ),

r=0

where δ(a) = 1 if a is true and δ(a) = 0 otherwise.
• Labels. The disease statuses y are ordinal variables too. To generate y, we use
a Gaussian process ordinal regression model [20] based the latent representation
U,
p(y|U) = p(y|f )p(f |U),
where
p(f |U) = N (f |0, K),
p(y|f ) =

2
X

δ(yi = r)δ(br ≤ fi < br+1 ),

r=0

where Kij = k(ui , uj ) is the cross-covariance between ui and uj . We can
choose k from a rich family of kernel functions such as linear, polynomial, and
Gaussian kernels to model relationships between the labels y and the latent features U.
Note that the labels y are linked to the data X and Z via the latent features U and
the projection matrices H and G. Due to the sparsity in H and G, essentially
only a few groups of variables in X and Z are selected to predict y. Note that
each of group is linked to a feature in U.
• Sparse Priors. Because we want to identify a few critical interactions between
different data sources, we use spike and slab prior distributions [17] to sparsify
the projection matrices G and H. Specifically, we use a p × k matrix Sg to
represent the selection of elements in G: if sij = 1, gij is selected and follows
a Gaussian prior distribution with variance σ12 ; if sij = 0, gij is not selected and
forced to almost zero (i.e., sampled from a Gaussian with a very small variance
σ22 ). Specifically, we have the following prior over G:
p(G|Sg , Πg ) =

p Y
k
Y

ij ij
p(gij |sij
g )p(sg |πg )

i=1 j=1

where
ij
2
ij
2
p(gij |sij
g ) = sg N (gij |0, σ1 ) + (1 − sg )N (gij |0, σ2 ),
ij
ij
p(sij
g |πg ) = πg

sij
g

ij

(1 − πgij )1−sg ,
5

2
2
where πgij in Πg is the probability of sij
g = 1, and σ1  σ2 (in our experiment,
2
2
−6
we set σ1 = 1 and σ2 = 1o ). To reflect our uncertainty about Πg , we assign
a Beta hyperprior distribution:

p(Πg |l1 , l2 ) =

p Y
k
Y

Beta(πgij |l1 , l2 ),

i=1 j=1

where l1 and l2 are hyperparameters. We set l1 = l2 = 1 in our experiments.
Similarly, H is sampled from
p(H|Sh , Πh ) =

q Y
k
Y

ij ij
p(hij |sij
h )p(sh |πh ),

i=1 j=1

where
ij
ij
2
2
p(hij |sij
h ) = sh N (hij |0, σ1 ) + (1 − sh )N (hij |0, σ2 ),
ij
ij
p(sij
h |πh ) = πh

sij
h

ij

(1 − πhij )1−sh ,

where Sh are binary selection variables and πhij in Πh is the probability of sij
h =
1. Again, we use a Beta hyperprior distribution:
p(Πh |d1 , d2 ) =

q Y
k
Y

Beta(πhij |d1 , d2 ),

i=1 j=1

where d1 and d2 are hyperparameters. We set d1 = d2 = 1 in our experiments.
Based on all these specifications, the joint distribution of our model, SHML, is
p(X, Z, y, U, G, Sg , Πg , η, C, H, Sh , Πh , f , )
= p(X|U, G, η)p(G|Sg )p(Sg |Πg )p(Πg |l1 , l2 )p(η|r1 , r2 )
· p(Z, C|U, H)p(H|Sh )p(Sh |Πh )p(Πh |d1 , d2 )
· p(y|f )p(f |U)p(U).

3

(1)

Estimation

Given the model specified in the previous section, now we present an efficient, principled method to estimate the latent features U, the projection matrices H and G, the
selection indicators Sg and Sh , the selection probabilities Πg and Πh , the variance η,
the auxiliary variables C for generating ordinal data Z, and the auxiliary variables f
for generating the labels y. In a Bayesian framework, this estimation task amounts to
computing their posterior distributions.
However, computing the exact posteriors turns out to be infeasible since we cannot calculate the normalization constant of the posteriors based on Equation (1). Thus,
we resort to a variational Bayesian Expectation Maximization (VB-EM) approach [24].
6

More specifically, in the E step, we approximate the posterior distributions of H, G, Sg , Sh , Πg ,
Πh , η, C and f by a factorized distribution
Q(H)Q(G)Q(Sg )Q(Sh )Q(Πg )Q(Πh )Q(η)Q(C)Q(f )
and then use the approximate distributions to compute expectations in the M step to
optimize the latent features U.
To obtain the variational approximation, we minimize the Kullback-Leibler (KL)
divergence between the approximate and the exact posteriors, KL(Q||P ) where P represents the exact joint posterior distributions. To this end, we use coordinate descent;
we update an approximate distribution, say, Q(H), while fixing the other approximate
distributions, and iteratively refine all the approximate distributions. The detailed updates are given in the following paragraphs.

3.1

Updating variational distributions for continuous data

For the continuous data X, the approximate distributions of the projection matrix G,
the noise variance η, the selection indicators Sg and the selection probabilities Πg are
Q(G) =

Q(Sg ) =

p
Y

N (gi ; λi , Ωi ),

i=1
p Y
k
Y

(2)

sij

ij

βijg (1 − βij )1−sg ,

(3)

Beta(πgij |˜l1ij , ˜l2ij ),

(4)

i=1 j=1

Q(Πg ) =

p Y
k
Y
i=1 j=1

Q(η) = Gamma(η|r̃1 , r̃2 ).

(5)

The mean and covariance of gi are calculated as follows:
Ωi = hηiUU> +

−1
1
1
diag(hsig i) + 2 diag(1 − hsig i)
,
2
σ1
σ2

λi = Ωi (hηiUx̃i ),
where h·i means expectation over a distribution, x̃i and sig are the transpose of the i-th
2
rows of X and Sg , hsig i = [βi1 , . . . , βik ]> , and hgij
i is the j-th diagonal element in
Ωi .
ij
The parameter βij in Q(sij
g ) is calculated as βij = 1/ 1 + exp(hlog(1 − πg )i −
2

σ
2
hlog(πgij )i + 21 log( σ12 ) + 12 hgij
i( σ12 − σ12 )) . The parameters of the Beta distribution
2
1
2
Q(π ij ) is given by ˜lij = βij + l1 and ˜lij = 1 − βij + l2 . The parameters of the
g

1

2

Gamma distribution Q(η) are updated as r̃1 = r1 +
tr(hGiUX> ) + 12 tr(UU> hG> Gi).

7

np
2

and r̃2 = r2 + 21 tr(XX> ) −

The moments required in the above distributions are calculated as hηi =

r̃1
r̃2

and

hlog(πgij )i = ψ(˜l1ij ) − ψ(˜l1ij + ˜l2ij ),
hlog(1 − πgij )i = ψ(˜l2ij ) − ψ(˜l1ij + ˜l2ij ),
hG> Gi =

p
X

Ωi + λ i λ >
i ,

i=1

hGi = [λ1 , . . . , λp ]> ,
where ψ(x) =

3.2

d
dx

(6)

ln Γ(x).

Updating variational distributions for ordinal data

For the ordinal data Z, we update the approximate distributions of the projection matrix H, the auxiliary variables C, the sparse selection indicators Sh and the selection
probabilities Πh . Specifically, the variational distributions of C and H are
Q(C) =

q Y
k
Y

Q(cij ),

(7)

i=1 j=1

Q(cij ) ∝ δ(bzij ≤ cij < bzij +1 )N (cij |c̄ij , 1),
Q(H) =

q
Y

N (hi ; γ i , Λi ),

(8)
(9)

i=1

where c̄ij = γ >
i uj and
Λi = UU> +

−1
1
1
diag(hsih i) + 2 diag(h1 − sih i)
,
2
σ1
σ2

γ i = Λi (Uhc̃i i),
where c̃i is the transpose of the i-th row of C.
The variational distributions of Sh and Πh are
Q(Sh ) =

q Y
k
Y

sij

ij

αijh (1 − αij )1−sh ,

(10)

Beta(πhij |d̃ij1 , d̃ij2 ),

(11)

i=1 j=1

Q(Πh ) =

q Y
k
Y
i=1 j=1


σ2
where αij = 1/ 1+exp(hlog(1−πhij )i−hlog(πhij )i+ 21 log( σ12 )+ 12 hh2ij i( σ12 − σ12 )) ,
2
1
2
i
>
2
˜ij
d˜ij
1 = αij + d1 , d2 = 1 − αij + d2 , hsh i = [αi1 , . . . , αik ] , and hhij i is the j-th
diagonal element in Λi .

8

The required moments for updating the above distributions can be calculated as
follows:
hc̃i i = [hci1 i, . . . , hcin i]> ,
hlog(π ij )i = ψ(d˜ij ) − ψ(d˜ij + d˜ij ),
hlog(1 −

h
πhij )i

1

1

2

˜ij ˜ij
= ψ(d˜ij
2 ) − ψ(d1 + d2 ),

hcij i = c̄ij −

N (bzij +1 |c̄ij , 1) − N (bzij |c̄ij , 1)
,
Φ(bzij +1 − c̄ij ) − Φ(bzij − c̄ij )

where Φ(·) is the cumulative distribution function of a standard Gaussian distribution.
Note that in Equation (8), Q(cij ) is a truncated Gaussian and the truncation is controlled by the observed ordinal data zij .

3.3

Updating variational distributions for labels

We update the variational distribution of the auxiliary variables f as follows:
Q(f ) =

n
Y

Q(fi ),

(12)

i=1

Q(fi ) ∝ δ(byi ≤ fi < byi +1 )N (fi |f¯i , σf2i ),

(13)

where
f¯i = Ki,¬i K−1
¬i,¬i hf¬i i,

(14)

σf2i = Ki,i − Ki,¬i K−1
¬i,¬i K¬i,i ,

(15)

where Ki,¬i is the covariance between ui and U¬i , K¬i,¬i is the covariance on U¬i
(U¬i = [u1 , · · · ui−1 , ui+1 , · · · un ]), hf¬i i = [hf1 i, · · · , hfi−1 i, hfi+1 i, · · · , hfn i]> ,
and each hfi i is
N (byi +1 |f¯i , σf2i ) − N (byi |f¯i , σf2i )
.
(16)
hfi i = f¯i − σf2i ·
b
−f¯i
byi −f¯i
Φ( yi +1
)
−
Φ(
)
σf
σf
i

i

Note that Q(fi ) is also a truncated Gaussian and the truncated region is decided by the
ordinal label yi . In this way, the supervised information from y is incorporated into
estimation of f and then estimation of the other quantities by the recursive updates.

3.4

Optimizing the latent representation U

After the expectations of the other variables are calculated, we optimize U by maximizing the following variational lower bound
1
F (U) = − tr(UU> ) + hηitr(X> hGiU)
2
1
1
1
− tr(hH> HiUU> ) − log|K| − tr(hff > iK−1 )
2
2
2
hηi
−
tr(hG> GiUU> ) + tr(hCi> hHiU) + constant,
(17)
2
9

where
hH> Hi =

p
X

Λi + γ i γ >
i ,

hHi = [h1 , . . . , hq ]> ,

(18)

i=1

hff > i = hf ihf i> − diag(hf i2 ) + diag(hf 2 i),
hfi2 i

2

= hfi i +
+ σf2i ·
− σf2i ·

(19)

σf2i

(byi − hfi i)N (byi |hfi i, σf2i )
Φ(

byi +1 −hfi i
)
σfi

− Φ(

byi −hfi i
)
σfi

(byi +1 − hfi i)N (byi +1 |hfi i, σf2i )
Φ(

byi +1 −hfi i
)
σfi

− Φ(

byi −hfi i
)
σfi

,

(20)

and the constant means a value independent of U so that it is irrelevant for optimizing
U. Note that we can optimize the dimension k by maximizing the full variational lower
bound of our model, which involves other quantities as well, such as hHi and hGi. To
save space, we do not present the long equation for the full lower bound (which can be
easily derived based on what we have presented).
The other required moments are given in Equations (6) and (16). We use the LBFGS algorithm to maximize the cost function F over U. The gradient of U is given
by
∂F
= hηihGi> X + hHi> hCi − I + hηihG> Gi
∂U

 ∂K
1
1
+ hH> Hi U − K−1 − K−1 hff > iK−1
.
2
2
∂U
Note that

∂K
∂U

(21)

depends on the form of the kernel function k(ui , uj ).

Algorithm 1 VB-EM for model estimation
1. Initialize U, the hyperparamters, and the
moments of all the approximate distributions.
2. Loop until convergence:
E-step Update all the approximate distributions according to (2-5, 7-11, 12-13).
M-step Use L-BFGS to optimize U.
3. Output U and all the approximate posterior
distributions.

3.5

Prediction

Let us denote the training data as Dtrain = {Xtrain , Ztrain ,
ytrain } and the test data as Dtest = {Xtest , Ztest }. The prediction task needs the latent
representation Utest for Dtest . There are two candidate strategies for obtaining Utest .
10

The first one is separate learning: we first learn the projection matrices from Dtrain ,
i.e., Q(H) and Q(G), and then fix them in the variational EM procedure on Dtest to
learn Utest . Note that there are no updates for ordinal label part on Dtest and the terms
regarding ordinal labels should also be removed from Equation (17) and (21). The
second strategy is joint learning, where we carry out variational EM simultaneously on
Dtrain and Dtest . A drawback of the first strategy is that the (distributions of) loading
matrices are fixed when learning latent representation Utest . Therefore, we adopt the
the second strategy; in other words, the variational EM algorithm uses all the data to
update the variation distributions, except Q(f ) where only labels in the training set are
used. After both Utest and Utrain are obtained from the M-step, we predict the labels
for test data as follows:


ftest = K Utest , Utrain K−1 Utrain , Utrain hftrain i,
(22)
i
ytest

=

R−1
X

i
r · δ(br ≤ ftest
< br+1 ),

(23)

r=0
i
where ytest
is the prediction for i-th test sample.

4

Related Work

The proposed SHML model is related to a broad family of probabilistic latent variable models, including probabilistic principle component analysis [25], probabilistic
canonical correlation analysis [26] and their extensions [27, 28, 29, 30].They all learn
a latent representation whose projection leads to the observed data. Recent studies on
probabilistic factor analysis methods put more focus on the sparsity-inducing priors to
the projection matrix. Among them, Guan et al. [27] used the Laplace prior, the Jeffrey’s prior, and the inverse-Gaussian prior; Archambeau & Bach [29] employed the
inverse-Gamma prior; and Virtanen et al. [30] used the Automatic Relevance Determination(ARD) prior. Despite their success, these sparsity-inducing priors have their own
disadvantages – they confound the degree of sparsity with the degree of regularization
on both relevant and irrelevant variables, while in practical settings there is little reason
that these two types of complexity control should be so tightly bounded together. Although the inverse-Gaussian prior and the inverse-Gamma prior provide more flexibility of controlling the sparsity, they suffer from being highly sensitive to the controlling
parameters and thus lead to unstable solutions. In contrast, our model adopts the spike
and slab prior, which has been recently used in multi-task multiple kernel learning [31],
sparse coding [18], and latent factor analysis [32]. Note that while our Beta priors over
the selection indicators lead to simple yet effective variational updates, the hierarchical
prior in [32] can better handle the selection uncertainty. Regardless what priors are
assigned to the spike and slab models, they generally avoid the confounding issue by
separately controlling the projection sparsity and the regularization effect over selected
elements.
SHML is also connected with many methods on learning from multiple sources
or views [33]. Multiview learning methods are often used to learn a better classifier
for multi-label classification – usually in text mining and image classification domains
11

– based on correlation structures among the training data and the labels [28, 30, 34].
However, in medical analysis and diagnosis, we meet two separate tasks – the association discovery between genetic variations and clinical traits, and the diagnosis on
patients. Our proposed SHML conducts these two tasks simultaneously: it employs the
diagnosis labels to guide association discovery, while leveraging the association structures to improve the diagnosis. In particular, the diagnosis procedure in SHML leads to
an ordinal regression model based on latent Gaussian process models. The latent Gaussian process treatment differentiates ours from multiview CCA models [35]. Moreover,
most multiview learning methods do not model the heterogeneous data types from different views, and simply treat them as continuous data. This simplification can degrate
the predictive performance. Instead, based on a probabilistic framework , SHML uses
suitable link functions to fit different types of data.

5

Experiments

In this section, we demonstrate the effectiveness of SHML on both synthetic and real
data for AD study.

5.1

Simulation Study

We first design a simulation study to examine the performance of SHML in terms of
(i) estimation accuracy in finding associations between two views and (ii) prediction
accuracy on ordinal labels.
Simulation data. To generate the ground truth, we set n = 200 (200 instances),
p = q = 40, and k = 5. We designed G, the 40×5 projection matrix for the continuous
data X, to be a block diagonal matrix; each column of G had 8 elements being ones
and the rest of them were zeros, ensuring each row with only one nonzero element. We
designed H, the 40 × 5 projection matrix for the ordinal data Z, to be a block diagonal
matrix; each of the first four columns of H had 10 elements being ones and the rest
of them were zeros, and the fifth column contains only zeros. We randomly generated
the latent representations U ∈ Rk×n with each column ui ∼ N (0, I). To generate
Z, we first sampled the auxiliary variables C with each column ci ∼ N (Hui , 1), and
then decided the value of each element zij in Z by the region cij falls in – in other
P2
words, zij = r=0 rδ(br < cij ≤ br+1 ) where b = {− inf, −1, 1, inf}. Similarly, to
generate y, we sampled the auxiliary variables f from N (0, U> U + I) and then each
yi was generated by p(yi |fi ) = δ(yi = 0)δ(fi ≤ 0) + δ(yi = 1)δ(fi > 0).
Comparative methods. We compared SHML with several state-of-the-art methods including (1) CCA [6], which finds the projection directions that maximize the
correlation between two views, (2) sparse CCA [36, 8], where sparse priors are put
on the CCA directions, and (3) Multiple Regression with lasso (MRLasso) [37] where
each column of the second view (Z) is regarded as the output of the first view (X). We
did not include results from the sparse probabilistic projection approach [29] because
it performed unstably in our experiments. Regarding the software implementation, we
used the built-in Matlab Matlab routine for CCA and the code by [36] for sparse CCA.

12

We implemented MRLasso based on the Glmnet package (cran.r-project.org/
web/packages/glmnet/index.html).
To compare accuracy on predicting labels y, we compared our method with the following ordinal or multinomial regression methods: (1) lasso for multinomial regression
[10], (2) elastic net for multinomial regression [11], (3) sparse ordinal regression with
the splike and slab prior, (4) CCA + lasso, for which we first ran CCA to obtain the
latent features H and then applied lasso to predict y, (5) CCA + elastic net, for which
we first ran CCA to obtain the projection matrices and then applied elastic net on the
projected data, (6) Gaussian Process Ordinal Regression (GPOR) [20], which employs
Gaussian processes to learn the latent function for ordinal regression, and (7) Laplacian Support Vector Machine (LapSVM) [38], a semi-supervised SVM classification
method. We used the Glmnet package for lasso and elastic net, the GPOR package by
[20], and the LapSVM package by [38]. For all the methods, we used 10-fold cross
validation to tune free parameters for each run; for example, we used extensive crossvalidation to choose the kernel form (Gaussian or Polynomials) and its parameters (the
kernel width or polynomial orders) for SHML, GPOR, and LapSVM. Note that all
these methods, except SHML, stack X and Z together into one data matrix and ignore
their heterogeneous nature.
Because alternative methods cannot learn the dimension automatically from the
data, for fair comparison, we provided the dimension of the latent representation to
all the methods we tested in our simulations. For each run in our experiment, we
partitioned the data into 10 subsets and used 9 of them for training and 1 subset for
testing. We repeated the procedure 10 times to generate the averaged results.
Results. To estimate linkage (i.e., interactions) between X and Z, we calculated the
cross covariance matrix GH> . We then computed the precision and the recall based
on the ground truth. The the precision-recall curves are shown in Figure 2. Clearly,
our method successfully recovered almost all the links and significantly outperformed
all the competing methods. This improvement may come from i) the use of the spike
and slab priors, which not only remove irrelevant elements in the projection matrices
but also avoid over-penalize the active association structures (the Laplace prior used
in sparse CCA does over penalize the relevant ones) and ii) more importantly, the
supervision from the labels y, which is probably the biggest difference between ours
and the other methods for the association study.
The prediction accuracies on unknown y and their standard errors are shown in
Figure 3. Our proposed SHML model achieves significant improvement over all the
other methods. In particular, it reduces the prediction error of elastic net (which ranks
the second best) by 25%, and reduces the error of LapSVM (which ranks the last), by
48%. Note that although utilizing the information from the unlabeled data, LapSMV
lacks the capability to utilize underlying interaction structures and sparsify the model
parameters, which may contribute its poor performance in the experiments.
In summary, the simulation results confirm the power of SHML in both discovering
true associations between heterogeneous data sources and predicting unknown labels.

13

1

SHML
0.9
0.8

Precision

0.7

Sparse CCA

0.6

MRLasso

0.5
0.4
0.3

CCA

0.2
0.1
0

0

0.2

0.4

0.6

0.8

1

Recall

Figure 2: The precision-recall curves for association discovery.

Precision

0.9
0.85
0.8

ne
t
ar
se
O
R
G
P
C
O
C
A+ R
C
La
C
A+
ss
El
o
as
tic
N
O
et
ur
M
od
el

o
ss

Sp

La

tic
as
El

La

pS

VM

0.75

Figure 3: Prediction accuracies on the simulation data. The results are averaged over 10 runs.

5.2

Study of Alzheimer’s Disease

We conducted association analysis and diagnosis of AD based on a dataset from Alzheimer’s
Disease Neuroimaging Initiative(ADNI). The ADNI study is a longitudinal multisite
observational study of elderly individuals with normal cognition, mild cognitive impairment, or AD. AD is the most common form of dementia with about 30 million
patients worldwide and payments for care are estimated to be $200 billion in 2012.1 .
In this analysis, we used SHML to study the associations of genotypes and brain atrophy measured by MRI and to predict the subject status (normal vs MCI vs AD). Note
that the labels are ordinal since the three states represent increasing severity levels of
the dementia.
The dataset was downloaded from http://adni.loni.ucla.edu/. After
removing missing data, it consists 618 subjects (183 normal, 308 MCI and 134 AD),
and for each patient, there are 924 SNPs (selected as the top SNPs to separate normal subjects from AD in ADNI) and 328 MRI features measuring the brain atrophies
in different brain regions based on cortical thickness, surface area or volume using
1 www.alz.org/downloads/facts_figures_2012.pdf

14

FreeSurfer software.
We compared SHML with the alternative methods on accuracy of predicting whether
a subject is in the normal or MCI or AD condition. We randomly split the dataset into
556 training and 62 test samples 10 times and ran all the competing methods on each
partition. As for the simulation study, we used the 10-fold cross validation for each
run to tune free parameters. In SHML, in order to determine dimension k for the latent
representation U, we computed the variational lower bounds as an approximation to
the model marginal likelihood (i.e., evidence), with various k values {10, 20, 40, 60}.
We chose the value with the largest approximate evidence, which led to k = 20 (see
Figure 4). Our experiments confirmed that with k = 20, SHML achieved highest
5

Lower Bound of Evidence

x 10
−7.2

−7.22

−7.24

−7.26

10

20
40
60
Dimension of latent features

Figure 4: The variational lower bound for the model marginal likelihood.
prediction accuracy, demonstrating the benefit of evidence maximization.
The accuracies for predicting unknown labels y and their standard errors are shown
in Figure 5. Our method achieved the highest prediction accuracy, higher than that of
the second best method, GP ordinal Regression, by 10% and than that of the worst
method, CCA+lasso, by 22%.

Precision

0.65
0.6
0.55
0.5

R
La
ss
El
o
as
tic
N
O
et
ur
M
od
el

R

+

A

+

C

O

PO
G

C

C

A

C

Sp

ar

se

N

et

so

tic

La
s

as

El

La

pS
VM

0.45

Figure 5: The prediction accuracy with standard errors on the real data.
We also examined the strongest associations discovered by SHML based on this
dataset. First of all, the ranking of MRI features in terms of their prediction power
of three different disease populations (normal, MCI and AD) demonstrate that most

15

(a)
1.5

0.5

0

−0.5

−1

S5
A
PD

2
2A
C
SL

T3

AP
T
M

N
W

N

SF

−1.5

A

2A
2

S5
PD

SL
C

T3

M
AP
T

(b)

CT std of L. PosteriorCingulate
CT std of R. InferiorTemporal
CT std of R. InferiorParietal
CT std of R. Fusiform
CT std of R. FrontalPole
CT std of R. IsthmusCingulate
CT std of R. CaudalMiddleFrontal
CT std of R. CaudalAnteriorCingulate
CT std of R. Bankssts
CT std of L. TransverseTemporal
CT std of L. Supramarginal
CT std of L. SuperiorTemporal
CT std of L. SuperiorParietal
CT std of L. SuperiorFrontal
CT std of L. RostralMiddleFrontal
CT std of R. LateralOccipital
CT std of L. Precuneus
CT std of R. LateralOrbitofrontal
CT std of R. MiddleTemporal
CT std of L. Precentral
CT std of L. RostralAnteriorCingulate
CT std of L. TemporalPole
CT Avg of L. RostralAnteriorCingulate
CT std of R. MedialOrbitofrontal
Vol (WMP) of R. CerebellumCortex
Vol (WMP) of CorpusCallosumPosterior

1

−1

N

2

S5
A

2A

PD

SL
C

T3

M

AP
T

N

N

W

SF

−1.5

−0.5

SF

−1

0

N

0
−0.5

0.5

(c)

Vol (CP) of R. CaudalAnteriorCingulate

Surf Area of R. CaudalAnteriorCingulate

MA
MAP3K
P 1
M 3K1
MA IER
P 3
BC 3K1
BCAR3
BCAR3
BCAR3
BCAR3
NC AR3
NCOA2
NCOA2
NCOA2
L OA
LR PPR2
1
LRRTM
XP RTM4
NP 4
EP
ST 3
13

1
0.5

Vol (WMP) of CorpusCallosumAnterior
Vol (WMP) of CorpusCallosumMidAnterior
Vol (WMP) of CorpusCallosumCentral
Vol (CP) of L. PosteriorCingulate
CT Avg of L. PosteriorCingulate
Vol (WMP) of L. CerebellumCortex
CT std of L. CaudalAnteriorCingulate
CT std of L. ParsOpercularis
CT std of L. Parahippocampal
CT std of L. Paracentral
CT std of L. MiddleTemporal
CT std of L. Lingual
CT std of L. LateralOrbitofrontal
CT std of L. LateralOccipital
CT std of L. IsthmusCingulate
CT std of L. InferiorTemporal
CT std of L. InferiorParietal
CT std of L. Fusiform
CT std of L. FrontalPole
CT std of L. Entorhinal
CT std of L. Cuneus
CT std of L. ParsOrbitalis
CT std of L. Postcentral
CT std of L. ParsTriangularis
CT std of L. CaudalMiddleFrontal

1

W

CT Avg of R. Pericalcarine
CT std of R. TemporalPole
CT Avg of L. CaudalAnteriorCingulate
Vol (CP) of L. CaudalAnteriorCingulate
CT Avg of R. RostralAnteriorCingulate
CT std of R. RostralAnteriorCingulate
CT std of R. Paracentral
CT std of L. Bankssts
CT std of R. Insula
CT std of R. Parahippocampal
CT std of L. Insula
CT std of R. ParsOrbitalis
CT std of R. ParsTriangularis
CT std of R. Postcentral
CT std of R. PosteriorCingulate
CT std of R. Precentral
CT std of R. Precuneus
CT std of R. TransverseTemporal
CT std of R. SuperiorTemporal
CT std of R. RostralMiddleFrontal
CT std of R. SuperiorFrontal
CT std of R. SuperiorParietal
CT std of R. ParsOpercularis
CT std of R. Supramarginal
Surf Area of R. RostralAnteriorCingulate

1.5

(d)

Figure 6: The estimated associations between MRI features and SNPs. In each subfigure, the MRI features are listed on the right and the SNP names are given at the
bottom.
of top ranked features are based on the cortical thickness measurement. On the other
hand, the features based on volume and surface area estimation of the same brain structures are less predictive. Particularly, thickness measurements of middle temporal lobe,
precuneus, and fusiform were found to be most predictive compared with other brain
regions. These findings are consistent with the memory-related function in these regions and findings in the literature for their prediction power of AD. We also found
that measurements of the same structure on the left and right side have similar weights,
indicating that the algorithm can automatically select correlated features in groups,
since no asymmetrical relationship has been found for the brain regions involved in
AD.
Secondly, the analysis of associating genotype to AD disease prediction also generated interesting results. Similar to the MRI features, SNPs that are in the vicinity
of each other often listed together, indicating the group selection characteristics of the
algorithm. For example, the top ranks SNPs are associated with a few genes including
PSMC1P12 (proteasome 26S subunit, ATPase), NCOA2 (The nuclear receptor coactivator 2), and WDR52(WD repeat domain 52), which have been studied intensively in
cancer research.
At last, biclustering of the gene-MRI association, as shown in Figure 6 reveal interesting pattern in terms of the relationship between genetic variations and brain at16

rophy measured by structural MRI. For example, the top ranks SNPs are associated
with a few genes including BCAR3 (Breast cancer anti-estrogen resistance protein 3)
and NCOA2, which have been studied more carefully in cancer research. One of the
genes associated with this set of SNPs is MATP (microtubule-associated protein tau),
which codes the tau gene that are associated closely with the AD. These findings reveal strong association between MATP gene and atrophy in the memory-related brain
regions. Moreover, the same set of SNPs are also highly associated with cingulate, but
in an opposite direction. These results indicate an opposite effect of genotype to the
cingulate region, which is part of the limbic system and involve in emotion formation
and processing, compared with other structures such as temporal lobe, which plays a
more important role in the formation of long-term memory.
In summary, SHML discovered synergistic predictive relationships between brain
atrophy, genetic variations and the disease status.

6

Conclusions

We have presented, SHML, a new Bayesian multiview learning framework. SHML
simultaneously finds key associations between data sources (i.e., genetic variations and
phenotypic traits) and to predict unknown ordinal labels. Experimental results on the
ADNI data indicate that SHML found biologically meaningful associations between
SNPs and MRI features and led to significant improvement on predicting the ordinal
AD stages over the alternative classification and ordinal regression methods. Although
we have focused on the AD study, we expect that SHML, as a powerful extension
of CCA, can be applied to a wide range of applications in biomedical research – for
example, eQTL analysis supervised by additional labeling information.

7

Acknowledgments

This work was supported by NSF IIS-0916443, NSF CAREER award IIS-1054903,
and the Center for Science of Information (CSoI), an NSF Science and Technology
Center, under grant agreement CCF-0939370.

References
[1] L. Consoli et al. QTL analysis of proteome and transcriptome variations for dissecting the
genetic architecture of complex traits in maize. Plant Mol Biol., 48(5):575–581, 2002.
[2] D. Hunter. Lessons from genome-wide association studies for epidemiology. Epidemiology, (3):363–367, 2012.
[3] S. Gandhi and N. Wood. Genome-wide association studies: the key to unlocking neurodegeneration? Nature Neuroscience, 13:789–794, 2010.
[4] J. Liu et al. Combining fMRI and SNP Data to Investigate Connections Between Brain
Function and Genetics Using Parallel ICA. Hum Brain Mapp., (1):1–30, 2009.
[5] H. Harold. Relations between two sets of variates. Biometrika, 28:321–377, 1936.

17

[6] F. Bach and M. Jordan. A probabilistic interpretation of canonical correlation analysis.
Technical report, UC Berkeley, 2005.
[7] E. Parkhomenko, D. Tritchler, and J. Beyene. Genome-wide sparse canonical correlation
of gene expression with genotypes. BMC Proc., (Suppl 1), 2007.
[8] M. Daniela and R. Tibshirani. Extensions of sparse canonical correlation analysis, with
applications to genomic data. Statistical Applications in Genetics and Molecular Biology,
383(1), 2009.
[9] X. Chen, H. Liu, and J. Carbonell. Structured sparse canonical correlation analysis. In
AISTATS’12, volume 22, pages 199–207, 2012.
[10] R. Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal
Statistical Society, Series B, 58:267–288, 1994.
[11] H. Zou and T. Hastie. Regularization and variable selection via the elastic net. Journal of
the Royal Statistical Society, Series B, 67:301–320, 2005.
[12] M. Yuan and Y. Lin. Model selection and estimation in regression with grouped variables.
Journal of the Royal Statistical Society, Series B, 68(1):49–67, 2007.
[13] D. MacKay. Bayesian interpolation. Neural Computation, 4:415–447, 1991.
[14] R. Neal. Bayesian Learning for Neural Networks. 1996.
[15] P. Yu et al. Enriching amnestic mild cognitive impairment populations for clinical trials: optimal combination of biomarkers to predict conversion to dementia. Journal of
Alzheimer’s Disease, 32(2):373–385, 2012.
[16] L. Shen, Y. Qi, et al. Sparse bayesian learning for identifying imaging biomarkers in AD
prediction. In Med Image Comput Comput Assist Interv., volume 13, pages 611–618, 2010.
[17] E. George and R. McCulloch. Approaches for bayesian variable selection. Statistica Sinica,
7(2):339–373, 1997.
[18] I. Goodfellow, A. Couville, and Y. Bengio. Large-scale feature learning with spike-andslab sparse coding. In ICML’12. 2012.
[19] S. Mohamed, K. Heller, and Z. Ghahramani. Bayesian and L1 approaches for sparse unsupervised learning. In ICML’02, 2012.
[20] W. Chu and Z. Ghahramani. Gaussian processes for ordinal regression. Journal of Machine
Learning Research, 6:1019–1041, 2005.
[21] M. Braskie et al. Plaque and tangle imaging and cognition in normal aging and Alzheimer’s
diseasee. Neurobiology of AgingArchives of Neurology, 31(10):1669–1678, 2010.
[22] J. Zhou, J. Liu, V. Narayan, and J. Ye. Modeling disease progression via fused sparse group
lasso. In KDD’12, pages 1095–1103, 2012.
[23] L. Yuan, Y. Wang, P. Thompson, V. Narayan, and J. Ye. Multi-source learning for joint
analysis of incomplete multi-modality neuroimaging data. In KDD’12, pages 1149–1157,
2012.
[24] M. Beal. Variational Algorithms for Approximate Bayesian Inference. Ph.d. thesis, Gatsby
Computational Neuroscience Unit, University College London, 2003.
[25] M. Tipping and C. Bishop. Probabilistic principal component analysis. Journal of The
Royal Statistical Society Series B-statistical Methodology, 61:611–622, 1999.
[26] F. Bach and M. Jordan. A probabilistic interpretation of canonical correlation analysis.
Technical report, UC Berkley, 2005.

18

[27] Y. Guan and J. Dy. Sparse probabilistic principal component analysis. Journal of Machine
Learning Research - Proceedings Track, 5:185–192, 2009.
[28] S. Yu, K. Yu, V. Tresp, H. Kriegel, and M. Wu. Supervised probabilistic principal component analysis. In KDD’06, pages 464–473, 2006.
[29] C. Archambeau and F. Bach. Sparse probabilistic projections. In Advances in Neural
Information Processing Systems 21, pages 73–80. 2009.
[30] S. Virtanen, A. Klami, and S. Kaski. Bayesian CCA via group sparsity. In ICML’11, pages
457–464, June 2011.
[31] M. Titsias and M. Lázaro-Gredilla. Spike and slab variational inference for multi-task and
multiple kernel learning. In NIPS’11, pages 2339–2347, 2011.
[32] C. Carvalho, J. Chang, J. Lucas, JR Nevins, Q. Wang, and M. West. High-dimensional
sparse factor modeling: applications in gene expression genomics. Journal of the American
Statistical Association, 103(484):1438–1456, 2008.
[33] D. Hardoon, G. Leen, S. Kaski, and J. Shawe-Taylor, editors. NIPS Workshop on Learning
from Multiple Sources, 2008.
[34] I. Rish, G. Grabarnik, G. Cecchi, F. Pereira, and G. Gordon. Closed-form supervised
dimensionality reduction with generalized linear models. In ICML’08, pages 832–839,
2008.
[35] J. Rupnik and J. Shawe-Taylor. Multi-view canonical correlation analysis. In KDD’10,
2010.
[36] L. Sun, S. Ji, and J. Ye. Canonical correlation analysis for multi-label classification: A
least squares formulation, extensions and analysis. IEEE Transactions on Pattern Analysis
and Machine Intelligence, 33(1):194–200, 2011.
[37] S. Kim, K. Sohn, and E. Xing. A multivariate regression approach to association analysis
of a quantitative trait network. Bioinformaics, 25(12):204–212, 2009.
[38] S. Melacci and B. Mikhail. Laplacian support vector machines trained in the primal. Journal of Machine Learning Research, 12:1149–1184, March 2011.

19

