# Kernel learning approaches for summarising and combining posterior
similarity matrices

Alessandra Cabassi, Sylvia Richardson, Paul D. W. Kirk

###### Abstract

Summary: When using Markov chain Monte Carlo (MCMC) algorithms to perform
inference for Bayesian clustering models, such as mixture models, the output
is typically a sample of clusterings (partitions) drawn from the posterior
distribution. In practice, a key challenge is how to summarise this output.
Here we build upon the notion of the posterior similarity matrix (PSM) in
order to suggest new approaches for summarising the output of MCMC algorithms
for Bayesian clustering models. A key contribution of our work is the
observation that PSMs are positive semi-definite, and hence can be used to
define probabilistically-motivated kernel matrices that capture the clustering
structure present in the data. This observation enables us to employ a range
of kernel methods to obtain summary clusterings, and otherwise exploit the
information summarised by PSMs. For example, if we have multiple PSMs, each
corresponding to a different dataset on a common set of statistical units, we
may use standard methods for combining kernels in order to perform integrative
clustering. We may moreover embed PSMs within predictive kernel models in
order to perform outcome-guided data integration. We demonstrate the
performances of the proposed methods through a range of simulation studies as
well as two real data applications.
Availability: R code is available at https://github.com/acabassi/combine-psms.
Contact<EMAIL_ADDRESS><EMAIL_ADDRESS>

Kernel learning approaches for summarising and combining posterior similarity
matrices

Alessandra Cabassi1, Sylvia Richardson1, and Paul D. W. Kirk1,2
_1 MRC Biostatistics Unit_
_2 Cambridge Institute of Therapeutic Immunology & Infectious Disease
University of Cambridge, U.K._

Preprint,

##  1 Introduction

Clustering techniques aim to partition a set of statistical units in such a
way that items in the same group are similar and items in different groups are
dissimilar. The definition of similarity between observations varies according
to the application at hand. In biomedical applications, for example,
clustering can be used to define groups of patients who have similar genotypes
or phenotypes and therefore are more likely to respond similarly to a certain
treatment and/or have similar prognoses. Statistical methods for clustering
can be divided into two main categories: heuristic approaches as, for
instance, $k$-means (Hartigan and Wong,, 1979) and hierarchical clustering
(Kaufman and Rousseeuw,, 1990), and model-based techniques such as mixture
models (McLachlan and Peel,, 2004). Here we are interested in particular in
Bayesian mixture models. Inference on the parameters of these models can be
performed either via deterministic approximate inference methods based on
variational inference (Blei et al.,, 2006, Bishop,, 2006) or using Markov
chain Monte Carlo (MCMC) schemes to sample from the posterior distribution
(Gelfand and Smith,, 1990). One of the advantages of Bayesian model-based
clustering is that, if infinite mixtures are used, one does not need to
specify the number of clusters a priori (Rasmussen,, 2000). On the other hand,
summarising the output of MCMC algorithms can be challenging, as it includes a
large number of partitions of the data sampled from the posterior
distribution. The labels of each mixture component can change at each
iteration of the MCMC, because the model likelihood is invariant under
permutations of the indices. This phenomenon is known as _label switching_
(Celeux et al.,, 2000). One way of summarising the cluster allocations that
circumvents this problem is to compute a posterior similarity matrix (PSM)
containing the probabilities of each pair of observations of belonging to the
same cluster (more details about Bayesian mixture models and PSMs can be found
in Section 2). In addition to that, however, one is often interested in
finding one clustering of the data that best matches the information contained
in the PSM (Fritsch and Ickstadt,, 2009).

We propose a new method to summarise the PSMs derived from the MCMC output of
Bayesian model-based clustering, showing that they are valid kernel matrices.
Consequently, we are able to use kernel methods such as kernel $k$-means to
find a summary clustering, and can make use of machine learning algorithms
developed for pattern analysis (see e.g. Shawe-Taylor and Cristianini,, 2004)
to combine the PSMs obtained from different datasets. We assume that we have a
similarity matrix for each dataset, from initial clustering analyses performed
for each dataset independently. We show how these multiple kernel learning
(MKL) techniques allow us to find a global clustering structure and assign
different weights to each PSM, depending on how much information each provides
about the clustering structure. We also show how we may include a response
variable, if it is available, in order to perform outcome-guided integrative
clustering. In particular, we show that, in the unsupervised framework, we can
use the localised multiple kernel $k$-means approach of Gönen and Margolin,
(2014). We further demonstrate that, if a response variable is available for
each data point, it is possible to incorporate this information using the
_simpleMKL_ algorithm of Rakotomamonjy et al., (2008) in order to perform
outcome-guided integrative clustering. Both algorithms assign a weight to each
PSM, that is output together with the global cluster assignment.

This work therefore contributes to the many ways of summarising the
clusterings sampled from the posterior distribution that have already been
proposed (Fritsch and Ickstadt,, 2009, Wade and Ghahramani,, 2018). Our
approach performs equally well in the case of one dataset and has the
advantage of being easily extended to the case of multiple data sources.

From a different perspective, this work also suggests a new rational way to
define kernels that are appropriate when the data are believed to possess a
clustering structure. Lanckriet et al., (2004) applied MKL methods to the
problem of genomic data fusion, trying different kernels for each data source.
However, how to define a kernel in general remained an open problem. Cabassi
and Kirk, (2020) suggested the use of methods based on multiple kernel
learning to summarise the output of multiple similarity matrices, each
resulting from applying consensus clustering (Monti et al.,, 2003) to a
different dataset. Similarly, using PSMs as kernels ensures that the
similarities between data points that we consider correctly reflect the
clustering structure present in the data.

The manuscript is organised as follows. In Section 2 we introduce the problem
of summarising PSMs. In Section 3 we recall the theory of kernel methods,
prove that PSMs are valid kernel matrices and explain how this allows us to
apply kernel methods to the problems mentioned above. We also present the
method that we use to choose the best number of clusters in the final summary.
In Section 4 we present some simulation studies demonstrating the performances
of kernel methods to summarise PSMs and to integrate multiple datasets in an
unsupervised and outcome-guided fashion. Finally, in Section 5 we introduce
our motivating examples and illustrate how our methodology can be applied to
it.

##  2 Background

In this section we briefly recall the concept of Bayesian mixture models
(Section 2.1) and we expose the problem of summarising the posterior
distribution on the cluster allocations (Section 2.2).

###  2.1 Bayesian mixture models

Mixture models assume that the data are drawn from a mixture of distributions:

$p(\boldsymbol{x})=\sum_{k=1}^{K}\pi_{k}f_{\boldsymbol{x}}(\boldsymbol{x}|\phi_{k}).$
(1)

where $f_{\boldsymbol{x}}$ is a parametric density that depends on the
parameter(s) $\phi_{k}$, and $\pi_{k}$ are the mixture weights, which must
satisfy $0\leq\pi_{k}\leq 1$ and $\sum_{k=1}^{K}\pi_{k}=1$. In the Bayesian
framework, we assign a prior distribution to the set of all parameters,
$\Pi=[\pi_{1},\dots,\pi_{k}]$ and $\Phi=[\phi_{1},\dots,\phi_{k}]$. In the
finite mixtures case, methods to estimate the posterior distributions when the
true number of mixture components is unknown include the MCMC-based algorithms
proposed by Ferguson, (1973) and Richardson and Green, (1997). MCMC approaches
also exist for so-called infinite mixture models (Rasmussen,, 2000), such as
Dirichlet process mixtures and their generalisations (Antoniak,, 1974).

###  2.2 Posterior similarity matrices

When using MCMC methods in order to perform Bayesian clustering on a dataset
$X=[\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{N}]$, one obtains a vector of
cluster assignments $\boldsymbol{c}^{(b)}=[c_{1}^{(b)},\dots,c_{N}^{(b)}]$
from the posterior distribution for each iteration $b=1,\dots,B$ of the
algorithm (see, for example, Neal,, 2000). From this, it is possible to obtain
a Monte Carlo estimate of the probability that observations $i$ and $j$ belong
to the same cluster as follows:

$P(c_{i}=c_{j}|X)\approx\frac{1}{B}\sum_{b=1}^{B}\mathbb{I}_{\\{c_{i}^{(b)}=c_{j}^{(b)}\\}}=\vcentcolon\Delta_{ij}.$
(2)

We denote by $\Delta$ the posterior similarity matrix that is the matrix that
has $ij$th entry $\Delta_{ij}$ equal to the right hand side of Equation (2).

Many ways to find a final clustering using the PSM $\Delta$ have been proposed
(Binder,, 1978, Wade and Ghahramani,, 2018, Dahl,, 2006, Fritsch and
Ickstadt,, 2009, Medvedovic and Sivaganesan,, 2002). A simple solution is to
choose, among the $\boldsymbol{c}^{(b)}$, the one that maximises the posterior
density. The problem with this approach is that many clusterings are
associated with very similar posterior densities (Fritsch and Ickstadt,,
2009). A more principled approach is to define a loss function
$L(\boldsymbol{c},\hat{\boldsymbol{c}})$ measuring the loss of information
that occurs when estimating the true clustering $\boldsymbol{c}$ with
$\hat{\boldsymbol{c}}$ (Binder,, 1978). The optimal clustering
$\boldsymbol{c}^{*}$ is then defined as the one minimising the posterior
expected loss:

$\boldsymbol{c}^{*}=\operatorname*{arg\,min}_{\hat{\boldsymbol{c}}}\mathbb{E}\left[L(\boldsymbol{c},\hat{\boldsymbol{c}})|X\right]=\operatorname*{arg\,min}_{\hat{\boldsymbol{c}}}\sum_{\boldsymbol{c}}L(\boldsymbol{c},\hat{\boldsymbol{c}})p\left(\boldsymbol{c}|X\right).$
(3)

Binder, (1978), for instance, suggested choosing the clustering
$\hat{\boldsymbol{c}}$ that minimises the loss function

$L_{\text{Binder}}(\boldsymbol{c},\hat{\boldsymbol{c}})=\sum_{i<j}\left[l_{1}\mathbb{I}(c_{i}=c_{j})\mathbb{I}(\hat{c}_{i}\not=\hat{c}_{j})+l_{2}\mathbb{I}(c_{i}\not=c_{j})\mathbb{I}(\hat{c}_{i}=\hat{c}_{j})\right],$
(4)

where $l_{1}$ and $l_{2}$ are positive constants determining whether assigning
observations that belong to the same clusters to different clusters is
penalised more highly than assigning observations that belong to different
clusters to the same cluster ($l_{1}/l_{2}>1$) or vice versa ($l_{1}/l2<1$)
and $\mathbb{I}$ is the indicator function. If $l_{1}=l_{2}$, then

$\boldsymbol{c}^{*}_{\text{Binder}}=\operatorname*{arg\,min}_{\hat{\boldsymbol{c}}}\sum_{i<j}\lvert\mathbb{I}_{\\{\hat{c}_{i}=\hat{c}_{j}\\}}-\Delta_{ij}\rvert.$
(5)

More recently, Wade and Ghahramani, (2018) proposed an alternative to Binder’s
loss function based on the variation of information of Meilă, (2007):

$L_{\text{VI}}(\boldsymbol{c},\hat{\boldsymbol{c}})=H(\boldsymbol{c})+H(\hat{\boldsymbol{c}})-2I(\boldsymbol{c},\hat{\boldsymbol{c}}),$
(6)

where $H(\boldsymbol{c})$ and $H(\hat{\boldsymbol{c}})$ represent the entropy
of clusterings $\boldsymbol{c}$ and $\hat{\boldsymbol{c}}$ respectively,
$I(\boldsymbol{c},\hat{\boldsymbol{c}})$ is the mutual information between
clusterings $\boldsymbol{c}$ and $\hat{\boldsymbol{c}}$. Defining by $C_{i}$
the set of elements that belong to cluster $i$ in clustering $\boldsymbol{c}$,
$\hat{C}_{j}$ the set of elements that belong to cluster $j$ in clustering
$\hat{\boldsymbol{c}}$, $n_{ij}=|C_{i}\cap\hat{C}_{j}|$,
$n_{i+}=\sum_{j}n_{ij}$, $n_{+j}=\sum_{i}n_{ij}$, this is

$L_{\text{VI}}=-\sum_{i=1}^{K_{N}}\frac{n_{i+}}{N}\log_{2}\left(\frac{n_{i+}}{N}\right)-\sum_{j=1}^{\hat{K}_{N}}\frac{n_{+j}}{N}\log_{2}\left(\frac{n_{+j}}{N}\right)-2\sum_{i=1}^{K_{N}}\sum_{j=1}^{\hat{K}_{N}}\frac{n_{ij}}{N}\log_{2}\left(\frac{n_{ij}N}{n_{i+}n_{+j}}\right),$
(7)

where $K_{N}$ and $\hat{K}_{N}$ represent the number of clusters in
$\boldsymbol{c}$ and $\hat{\boldsymbol{c}}$ respectively.

Dahl, (2006) advanced the idea to choose, among all the clustering vectors
$\boldsymbol{c}^{(b)}$, the one that minimises the least-squared distance to
the PSM:

$\boldsymbol{c}^{*}_{\text{Dahl}}=\operatorname*{arg\,min}_{\hat{c}}\sum_{i<j}\left(\mathbb{I}_{\\{\hat{c}_{i}=\hat{c}_{j}\\}}-\Delta_{ij}\right)^{2}.$
(8)

This turned out to be equivalent to minimising Binder’s loss function
(Equation 5). Fritsch and Ickstadt, (2009) improved on the methods of Binder
and Dahl by maximising the posterior expected adjusted Rand index of Hubert
and Arabie, (1985).

Moreover, Medvedovic and Sivaganesan, (2002) applied the complete linkage
approach of Everitt, (1993) to the matrix of pseudo-distances $1-\Delta$,
while Molitor et al., (2010) used the partitioning around medoids algorithm of
Kaufman and Rousseeuw, (1990).

All of these methods are only applicable with one PSM. In what follows we
describe a new way to find a clustering using posterior similarity matrices
that also allows us to summarise multiple similarity matrices $\Delta_{m}$ and
find a global clustering.

##  3 Methods

In this section we explain how the output of the MCMC algorithms for Bayesian
mixture models can be modified in order to obtain valid kernel matrices
(Section 3.1). After that, in Section 3.2, we show that this allows us to (i)
use the kernel $k$-means algorithm to summarise a PSM and find a summary
clustering of the data; (ii) combine multiple PSMs to perform integrative
clustering of multiple datasets using an extension of kernel $k$-means; (iii)
use an external response variable to determine the weight of each dataset in
the integrative setting, by using predictive kernel methods such as Support
Vector Machines (SVMs). In the last part of the section we explain how we
choose the number of clusters in the final clustering.

###  3.1 Kernels

Kernel methods allow non-linear relationships between data points with low
computational complexity to be modelled, thanks to the so-called _kernel
trick_ (Shawe-Taylor and Cristianini,, 2004). For this reason, kernel methods
have been widely used to extend many traditional algorithms to the non-linear
framework, such as the PCA (Schölkopf et al.,, 1998), linear discriminant
analysis (Mika et al.,, 1999, Roth and Steinhage,, 2000, Baudat and Anouar,,
2000) and ridge regression (Friedman et al.,, 2001, Shawe-Taylor and
Cristianini,, 2004).

Kernel methods proceed by embedding the observations into a higher-dimensional
feature space $\mathcal{X}$ endowed with an inner product
$\langle\cdot,\cdot\rangle_{\mathcal{X}}$ and induced norm
$\lVert\cdot\rVert_{\mathcal{X}}$, making use of a map
$\phi:\mathbb{R}^{P}\to\mathcal{X}$.

######  Definition 3.1.

A _positive definite kernel_ or, more simply, a _kernel_ $\delta$ is a
symmetric map $\delta:\mathbb{R}^{P}\times\mathbb{R}^{P}\to\mathbb{R}$ for
which for all
$\mathbf{x}_{1},\mathbf{x}_{2},\dots,\mathbf{x}_{N}\in\mathbb{R}^{P}$, the
matrix $\Delta$ with entries
$\Delta_{ij}=\delta(\mathbf{x}_{i},\mathbf{x}_{j})$ is positive semi-definite.
The matrix $\Delta$ is called the _kernel matrix_ or _Gram matrix_.

(The domain of the map $\phi$ is set to be $\mathbb{R}^{P}$ for simplicity,
but the definition also holds for more general sets of departure for $\phi$.)
Using this definition, it is possible to prove that any inner product of
feature maps gives rise to a kernel:

######  Proposition 3.1.

The map $\delta:\mathbb{R}^{P}\times\mathbb{R}^{P}\to\mathbb{R}$ defined by
$\delta(\mathbf{x},\mathbf{x}^{\prime})=\langle\phi(\mathbf{x}),\phi(\mathbf{x}^{\prime})\rangle_{\mathcal{X}}$
is a kernel.

Moreover, using Mercer’s theorem, it can be shown that for any positive semi-
definite kernel function, $\delta$, there exists a corresponding feature map,
$\phi$ (see e.g. Vapnik,, 1998). That is,

######  Theorem 3.1.

For each kernel $\delta$, there exists a feature map $\phi$ taking value in
some inner product space $\mathcal{X}$ such that
$\delta(\mathbf{x},\mathbf{x}^{\prime})=\langle\phi(\mathbf{x}),\phi(\mathbf{x}^{\prime})\rangle_{\mathcal{X}}$.

In practice, it is therefore often sufficient to specify a positive semi-
definite kernel matrix, $\Delta$, in order to allow us to apply kernel methods
such as those presented in the following sections. For a more detailed
discussion of kernel methods, see e.g. Shawe-Taylor and Cristianini, (2004).

###  3.2 Identifying PSMs as kernel matrices

It has been shown elsewhere that co-clustering matrices are valid kernel
matrices (Cabassi and Kirk,, 2020). We show here that this result also holds
for PSMs, and hence they can be used as input for any kernel-based model. PSMs
are convex combinations of co-clustering matrices $C^{(b)}$, where each matrix
$C^{(b)}$ is defined as follows:

$C_{ij}^{(b)}=\left\\{\begin{array}[]{ll}1&\text{if
}c^{(b)}_{i}=c^{(b)}_{j},\\\ 0&\textit{otherwise},\end{array}\right.$

where $C_{ij}^{(b)}$ indicates whether the statistical units $i$ and $j$ are
assigned to the same cluster at iteration $b$ of the MCMC chain. Indicating by
$K$ the total number of clusters and reordering the rows and column, each
$C_{ij}^{(b)}$ can be written as a block-diagonal matrix where every block is
a matrix of ones:

$C^{(b)}=\left[\begin{array}[]{c c c c c}J_{1}&0&0&\dots&0\\\
0&J_{2}&0&\dots&0\\\ \vdots&&&\ddots&\vdots\\\
0&0&0&\dots&J_{K}\end{array}\right]$ (9)

where $J_{k}$ is an $n_{k}\times n_{k}$ matrix of ones, with $n_{k}$ being the
number of items in cluster $k$. The eigenvalues of a block diagonal matrix are
simply the eigenvalues of its blocks, which, in this case, are nonnegative.
Therefore all $C^{(b)}$, with $b=1,\dots,B$, are positive semidefinite. Now,
if $\lambda$ is a nonnegative scalar, and $C$ is positive semidefinite, then
$\lambda C$ is also positive semidefinite. Moreover, the sum of positive
semidefinite matrices is a positive semidefinite matrix. Therefore, given any
set of nonnegative $\lambda_{b}$, $b=1,\dots,B$,
$\sum_{b=1}^{B}\lambda_{b}C^{(b)}$ is positive semidefinite. We can conclude
that any PSM is positive definite.

###  3.3 Summarising and combining PSMs using kernel methods

Here we show that the fact that all posterior similarity matrices are valid
kernels allows us to use kernel methods to find a clustering of the data that
summarises a sample of clusterings
$\boldsymbol{c}^{(1)},\dots,\boldsymbol{c}^{(B)}$ from the posterior
distribution of an MCMC algorithm for Bayesian clustering. To do this we
suggest to use an extension of the well-known $k$-means algorithm that only
needs as input a kernel matrix.

Moreover, this method can be easily extended to allow us to combine multiple
PSMs. This can be a useful feature under many circumstances. For instance, as
it is the case in our motivating examples, one could have different types of
information relative to the same statistical observations. In this situation,
it may be appropriate to define and fit different mixture models on each data
type, and then summarise the two posterior samples of clusterings at a later
stage. This can be achieved by using multiple kernel $k$-means algorithms,
that allow us to combine multiple kernels to find a global clustering. On top
of that, these techniques also assign different weights to each kernel. These
can be used to assess how much each dataset contributed to the final
clustering and therefore to get an idea of how much information is present in
each data type about the clustering structure.

The problem with combining multiple kernels is, however, that it is not always
clear whether they all have the same clustering structure. To overcome this
issue, we also propose an outcome-guided algorithm to summarise multiple PSMs.
The idea is that, instead of choosing the weight of each kernel in an
unsupervised way, if we have a variable available which is closely related to
the outcome of interest, we should weight more highly the kernels in which
statistical units that have similar outcomes are closer to each other. In
mathematical terms, this corresponds to using support vector machines to find
the kernel weights, where the response variable is our proxy for the outcome.

####  3.3.1 Summarising PSMs

In order to illustrate our method for summarising PSMs, first we recall the
main ideas behind $k$-means clustering and then we present its extension to
the kernel framework.

#####  $k$-means clustering

$k$-means clustering is a widely used clustering algorithm, first introduced
by Steinhaus, (1956). Let $\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{N}$
indicate the observed dataset, with $\boldsymbol{x}_{n}\in\mathbb{R}^{p}$ and
$z_{nk}$ be the corresponding cluster labels, where $\sum_{k}z_{nk}=1$ and

$z_{nk}=\left\\{\begin{array}[]{l l}1,&\text{if }\boldsymbol{x}_{n}\text{
belongs to cluster }k,\\\ 0,&\text{otherwise}.\\\ \end{array}\right.$ (10)

We denote by Z the $N\times K$ matrix with $ij$th element equal to $z_{ij}$.
The goal of the $k$-means algorithm is to minimise the sum of all squared
distances between the data points $\boldsymbol{x}_{n}$ and the corresponding
cluster centroid $\boldsymbol{m}_{k}$. The optimisation problem is

(11a) (11b) (11c)

where $\lVert\cdot\rVert_{2}$ indicates the Euclidean norm.

##### Kernel $k$-means clustering

Now we can show how the kernel trick works in the case of the $k$-means
clustering algorithm (Girolami,, 2002). Redefining the objective function of
Equation (11) based on the distances between observations and cluster centres
in the feature space $\mathcal{X}$, the optimisation problem becomes:

(11la) (11lb) (11lc)

where we indicated by $\tilde{\boldsymbol{m}}_{k}$ the cluster centroids in
the feature space $\mathcal{X}$. Using this kernel, each term of the sum in
Equation (11l) can be written as

$\displaystyle\lVert\phi(\boldsymbol{x}_{n})-\tilde{\boldsymbol{m}}_{k}\rVert^{2}_{\mathcal{X}}$
$\displaystyle=\langle\phi(\boldsymbol{x}_{n})-\tilde{\boldsymbol{m}}_{k},\phi(\boldsymbol{x}_{n})-\tilde{\boldsymbol{m}}_{k}\rangle_{\mathcal{X}}$
(11lm)
$\displaystyle=\langle\phi(\boldsymbol{x}_{n}),\phi(\boldsymbol{x}_{n})\rangle_{\mathcal{X}}-\frac{2}{N_{k}}\sum\limits_{i=1}^{N}z_{ik}\langle\phi(\boldsymbol{x}_{n}),\phi(\boldsymbol{x}_{i})\rangle_{\mathcal{X}}$
(11ln)
$\displaystyle\quad+\frac{1}{N_{k}^{2}}\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{N}z_{ik}z_{jk}\langle\phi(\boldsymbol{x}_{i}),\phi(\boldsymbol{x}_{j})\rangle_{\mathcal{X}}$
(11lo)
$\displaystyle=\delta(\boldsymbol{x}_{n},\boldsymbol{x}_{n})-\frac{2}{N_{k}}\sum\limits_{i=1}^{N}z_{ik}\delta(\boldsymbol{x}_{n},\boldsymbol{x}_{i})+\frac{1}{N_{k}^{2}}\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{N}z_{ik}z_{jk}\delta(\boldsymbol{x}_{i},\boldsymbol{x}_{j}).$
(11lp)

Therefore, we do not need to evaluate the map $\phi$ at every point
$\boldsymbol{x}_{i}$ to compute the objective function of Equation (11l).
Instead, we just need to know the values of the kernel evaluated at each pair
of data points $\delta(\boldsymbol{x}_{i},\boldsymbol{x}_{j})$,
$i,j=1,\dots,N$. This is what is commonly referred to as the kernel trick.

We have seen how kernels can be used to perform $k$-means clustering. Now, if
we have a sample of clusterings from the posterior, we can easily exploit this
technique to find a summary clustering. Once we compute our PSM $\Delta$, this
will be our kernel matrix. The clustering of interest will then be the one
given by kernel $k$-means in the form of $z_{nk}$, $n=1,\dots,N$,
$k=1,\dots,K$. The number of clusters $K$ is chosen as explained in Section
3.4.

####  3.3.2 Combining PSMs to perform integrative clustering

To combine multiple PSMs relative to the same statistical units, all we need
to do is to use the extension of kernel $k$-means to the case of multiple
kernels.

##### Multiple kernel $k$-means clustering

Gönen and Margolin, (2014) extended the kernel $k$-means approach to the case
of multiple kernels. We consider multiple datasets $X_{1},\dots,X_{M}$ each
with a different mapping function $\phi_{m}:\mathbb{R}^{P}\to\mathcal{X}_{m}$
and corresponding kernel
$\delta_{m}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})=\langle\phi_{m}(\boldsymbol{x}_{i}),\phi_{m}(\boldsymbol{x}_{j})\rangle_{\mathcal{X}_{m}}$
and kernel matrix $\Delta_{m}$. Then, if we define

$\phi_{\Theta}(\boldsymbol{x}_{i})=[\theta_{i1}\phi_{1}(\boldsymbol{x}_{i})^{\prime},\theta_{i2}\phi_{2}(\boldsymbol{x}_{i})^{\prime},\dots,\theta_{iM}\phi_{M}(\boldsymbol{x}_{i})^{\prime}]^{\prime},$
(11lq)

where $\Theta\in\mathbb{R}^{N\times M}_{+}$ is a vector of kernel weights such
that $\theta_{im}$ is the weight of observation $\boldsymbol{x}_{i}$ in
dataset $m$ and $\sum_{m}\theta_{im}=1$ for all $i=1,\dots,N$ and
$\theta_{im}\geq 0$, the kernel function of this multiple feature problem is a
convex sum of the single kernels:

$\displaystyle\delta_{\Theta}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})$
$\displaystyle=\langle\phi_{\Theta}(\boldsymbol{x}_{i}),\phi_{\Theta}(\boldsymbol{x}_{j})\rangle_{\mathcal{X}_{m}}$
(11lr)
$\displaystyle=\sum_{m=1}^{M}\theta_{im}\theta_{jm}\langle\phi_{m}(\boldsymbol{x}_{i}),\phi_{m}(\boldsymbol{x}_{j})\rangle_{\mathcal{X}_{m}}$
(11ls)
$\displaystyle=\sum_{m=1}^{M}\theta_{im}\theta_{jm}\delta_{m}(\boldsymbol{x}_{i},\boldsymbol{x}_{j}).$
(11lt)

We denote the corresponding kernel matrix by $\Delta_{\Theta}$. The
optimisation strategy proposed by Gönen and Margolin, (2014) is based on the
idea that, for some fixed vector of weights $\boldsymbol{\theta}$, the problem
is equivalent to the one of Equation (11l), where we had only one kernel.
Therefore, they develop a two-step optimisation strategy: (1) given a fixed
matrix of weights $\Theta$, solve the optimisation problem as in the case of
one kernel, with kernel matrix $\delta_{\Theta}$ and then (2) minimise the
objective function with respect to the kernel weights, keeping the assignment
variables fixed. This is a convex quadratic programming (QP) problem that can
be solved with any standard QP solver up to a moderate number of kernels $M$.

Similarly to before, once we have defined the kernels $\Delta_{m}$ to be equal
to each of our PSMs, the labels found through multiple kernel $k$-means
constitute the clustering that we are looking for. Moreover, the kernel
weights give us an indication of how each PSM contributed to the final
clustering.

DatasetsStep 1MCMC, CC, …KernelsStep 2Multiple kernel $k$-meansCluster labels
Figure 1: Schematic representation of the MKL-based integrative clustering
approach. Each colour indicates a different dataset/kernel. First, a mixture
model is fit using MCMC on each dataset separately. The resulting PSMs are
valid kernels that can be used as input to kernel $k$-means to find a global
clustering of the data. Note that, as shown elsewhere (Cabassi and Kirk,,
2020), other similarity matrices, such as the similarity matrices given by
consensus clustering (CC), define valid kernels. For this reason, they can be
combined with PSMs via MKL.

####  3.3.3 Embedding PSMs within predictive kernel machines to perform
outcome-guided integration

Suppose now that, in addition to the posterior similarity matrices
$\Delta_{1},\dots,\Delta_{M}$, we also have a categorical response variable
$y_{n}$ associated with each observation $\boldsymbol{x}_{n}$. As we explained
above, we would like to use this information to guide our clustering
algorithm. We can use the _simpleMKL_ algorithm described in the remainder of
this section to find the kernel weights $\theta_{1},\dots,\theta_{M}$ and then
use kernel $k$-means on the weighted kernel
$\Delta=\sum_{m}\theta_{m}\Delta_{m}$ to find the final clustering (Figure 2).

##### Support vector machines

We briefly recall here the concept of the support vector machine (Boser et
al.,, 1992) that is widely used for solving problems in classification and
regression (Schölkopf and Smola,, 2001, Bishop,, 2006).

In its simplest form, this method is applied to a binary classification
problem, in which the data points
$\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{N}\in\mathbb{R}^{P}$ in the training
set are assigned to two classes indicated by the target values
$y_{n}\in\\{-1,1\\},n=1,\dots,N$. We consider a feature map
$\phi:\mathbb{R}^{P}\to\mathcal{X}$ and the associated kernel
$\delta(\cdot,\cdot):\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ such that
$\delta(\boldsymbol{x}_{i},\boldsymbol{x}_{j})=\langle\phi(\boldsymbol{x}_{i}),\phi(\boldsymbol{x}_{j})\rangle_{\mathcal{X}}$.
Suppose that there exist some values of $\alpha_{n}$ and $b$ such that

$f(\boldsymbol{x})=\sum_{n=1}^{N}\alpha_{n}\delta(\boldsymbol{x},\boldsymbol{x}_{n})$
(11lu)

satisfies $f(\boldsymbol{x}_{n})+b>0$ if $y_{n}=1$ and
$f(\boldsymbol{x}_{n})+b<0$ otherwise. $f$ is a function that lives in a
function space $\mathcal{H}$ endowed with the norm
$\lVert\cdot\rVert_{\mathcal{H}}$. Then, this function can be used to classify
new data points $\boldsymbol{x}$ according to the sign of
$f(\boldsymbol{x})+b$. For support vector machines, the parameters
$\alpha_{n}$ and $b$ are chosen so as to maximise the margin, i.e. the
distance between the decision boundary given by Equation (11lu) and the point
$\boldsymbol{x}_{n}$ that is closest to the boundary. It can be shown that
this can be achieved by solving the quadratic programming problem (see e.g.
Bishop,, 2006, Rakotomamonjy et al.,, 2008)

(11lva)

However, in real applications, it is usually not possible to separate the two
classes perfectly. Hence, in order to take into account misclassifications, it
is necessary to introduce a penalty term that is linear with respect to the
distance of the misclassified points to the classification boundary (Bennett
and Mangasarian,, 1992). To this end, we define a variable $\xi_{n}$ (known as
a slack variable) for each data point such that

$\xi_{n}=\left\\{\begin{array}[]{l l}0,&\text{if }\boldsymbol{x}_{n}\text{ is
correctly classified,}\\\ \lvert
y_{n}-f(\boldsymbol{x}_{n})\rvert,&\text{otherwise}.\end{array}\right.$
(11lvy)

The optimisation problem of Equation (11lv) then becomes

(11lvza) (11lvzb)

where $\lambda>0$ is a parameter that controls the penalisation of
misclassifications. The objective functions (11lv) and (11lvz) are quadratic,
so any local optimum is also a global optimum. One of the most popular
approaches to solve this type of problems is _sequential minimal optimisation_
(Platt,, 1999). For more details about SVMs see e.g. Bishop, (2006).

##### Multiple kernel learning for SVMs

In the multiple kernel learning framework for SVMs, we consider $M$ different
feature representations, with mapping functions $\phi_{m}$ and corresponding
kernel functions $\delta_{m}$ and feature spaces $\mathcal{X}_{m}$. We
substitute the kernel $\delta$ of Equation (11lu) with a convex combination of
kernels $\delta_{m}$ (Lanckriet et al.,, 2004):

$f(\boldsymbol{x})+b=\sum_{m=1}^{M}\theta_{m}f_{m}(\boldsymbol{x})+b$
(11lvzaa)

where $\theta_{m}\geq 0$, $\sum_{m}\theta_{m}=1$ and
$f_{m}=\sum_{n}\delta_{m}(\boldsymbol{x},\boldsymbol{x}_{n})$. Rakotomamonjy
and Bach, (2007) proposed then to solve the optimisation problem

(11lvzaba) (11lvzabb) (11lvzabc) (11lvzabd)

using the convention that $x/0=0$ if $x=0$ and $\infty$ otherwise. The
algorithm of Rakotomamonjy and Bach takes the name of _simpleMKL_ and is based
on the idea that one can iteratively solve a standard SVM problem (11lvz) for
a fixed value of $\boldsymbol{\theta}$ and then update the vector of weights
$\boldsymbol{\theta}$ using the gradient descent method on the objective
function $J(\boldsymbol{\theta})$. Since the objective function is smooth and
differentiable with Lipschitz gradient, it can be easily optimised with the
reduced gradient algorithm (Luenberger and Ye,, 1984, Chapter 11). If the
standard SVM problem is solved exactly at each iteration, then convergence to
the global optimum is guaranteed (Luenberger and Ye,, 1984).

##### Multiclass multiple kernel learning

SVMs can be used also when the target value $y_{n}$ takes more than two
different values. The most commonly used approaches are called one-versus-one
(Knerr et al.,, 1990) and one-versus-the-rest (Vapnik,, 1998). In the first
one, we consider in turn each class as the “positive” case, and all the others
as the “negative” cases. This way, we construct $K$ different classifiers and
then assign a new observation $\boldsymbol{x}$ using

$y(\boldsymbol{x})=\max\limits_{k\in\\{1,\dots,K\\}}y_{k}(\boldsymbol{x})$
(11lvzabac)

The second approach is to train one SVM for each pair of classes and then
assign a point $\boldsymbol{x}$ to the class to which it is assigned more
often.

Rakotomamonjy et al., (2008) extended the _simpleMKL_ algorithm to the case of
a response with $K>2$ classes. Both these approaches can be used with the
simpleMKL algorithm, defining a new cost function $J(\boldsymbol{\theta})$ as
the sum of all the cost functions of the partial SVMs
$J_{s}(\boldsymbol{\theta})$:

$J(\boldsymbol{\theta})=\sum_{s\in\mathcal{S}}J_{s}({\boldsymbol{\theta}})$
(11lvzabad)

where $\mathcal{S}$ indicates the set of all partial SVMs and each $J_{s}$ is
defined as in Equation (11lvzab).

DatasetsStep 1MCMC, CC, …KernelsResponseStep 2SimpleMKL Weighted kernel Step
3Kernel $k$-meansCluster labels Figure 2: Schematic representation of the MKL-
based outcome-guided integrative clustering approach. Each colour indicates a
different dataset/kernel. First, a mixture model is fit on each dataset
separately. The resulting PSMs are valid kernels that can be used as input to
simpleMKL, if a response variable is available, to find a global clustering of
the data.

It is important to note that none of these approaches explicitly rely on the
fact that the $\Delta_{m}$ are posterior similarity matrices. Hence, any other
type of matrix $\Delta_{m}$ can be used, as long as it is symmetric, positive
semi-definite and the entries $\Delta^{m}_{ij}$ can be interpreted as some
measure of the similarity between $\boldsymbol{x}_{i}$ and
$\boldsymbol{x}_{j}$.

###  3.4 Choosing the number of clusters in the final summary

Many possible approaches have been proposed to choose the number of clusters
(see e.g. Milligan and Cooper,, 1985, Tibshirani et al.,, 2001, Yeung et al.,,
2001, Dudoit and Fridlyand,, 2002). We focus on the so-called silhouette, a
measure of compactness of the clustering structure, proposed by Rousseeuw,
(1987). There are two ways of defining the silhouette of a cluster, based
respectively on the similarities and dissimilarities between the data. Here we
briefly explain the former.

Given some cluster assignment labels $\boldsymbol{c}=[c_{1},\dots,c_{N}]$ and
some measure of the dissimilarity between the data points $\Delta_{ij}$ for
all $i,j=1,\dots,N$, we can define the following quantities: $a_{n}$ is the
average similarity of $\boldsymbol{x}_{n}$ to all objects of cluster $c_{n}$
and, for each $c_{i}\not=c_{n}$, $\Delta_{n,c_{i}}$ is the average similarity
of $n$ to all objects belonging to cluster $c_{i}$. Moreover, let us indicate
by $b_{n}$ the maximum $\Delta_{n,c_{i}}$ over all $i$ such that
$c_{i}\not=c_{n}$. Then, for each observation $n=1,\dots,N$, we can calculate

$s_{n}=\left\\{\begin{array}[]{ll}1-a_{n}/b_{n},&\text{if }a_{n}<b_{n},\\\
0,&\text{if }a_{n}=b_{n},\\\ a_{n}/b_{n}-1,&\text{if
}a_{n}>b_{n}.\end{array}\right.$ (11lvzabah)

This takes values between $-1$ and $1$, with higher values indicating that
$\boldsymbol{x}_{i}$ is well-clustered and negative values suggesting that
$\boldsymbol{x}_{i}$ has been misclassified.

Thus, we run our algorithms for combining the posterior similarity matrices
with different number of clusters from $K_{\text{min}}$ to $K_{\text{max}}$.
We consider the overall average silhouette width
$\bar{s}=\frac{1}{N}\sum_{n=1}^{N}s_{n}$ as a measure of the compactness of
clusters and we choose the value of $K$ that gives the highest value of
$\bar{s}$.

##  4 Simulation examples

Here we show how the methods presented above perform in practice. In Section
4.1 we explain how we generate the synthetic datasets for the simulation
studies. In Section 4.2 we show that the kernel $k$-means approach applied to
a posterior similarity matrix derived from a single dataset performs similarly
to standard clustering methods. Additionally, to assess the MKL-based
integrative clustering approaches described in Sections 3.3.2 and 3.3.3, we
perform a range of simulation studies; the results are presented in Section
4.3.

###  4.1 Synthetic datasets

We generate four synthetic datasets, each composed of data belonging to six
different clusters of equal size. Each observation
$\boldsymbol{x}_{n}^{(k)}\in\\{0,1,2\\}^{10}$ belonging to cluster $k$ is
drawn from a multivariate categorical distribution such that, for each
covariate $j=1,\dots,10$,

$x_{nj}^{(k)}\sim\text{Categorical}(\pi_{1k},\pi_{2k},\pi_{3k}),$ (11lvzabai)

where $\pi_{ik}$, $i=1,2,3$ are such that $\pi_{ik}=w\rho_{ik}+(1-w)/3$, with
$[\rho_{1k},\rho_{2k},\rho_{3k}]\sim\text{Dirichlet}(0.01)$. Each dataset has
a different value of $w\in[0,1]$. Higher values of $w$ give clearer clustering
structures. The response variable is binary, with
$P(y_{n}=1|z_{n}=k)=\theta_{k}$, where
$\theta_{k}\in\left\\{0.01,0.1,0.15,0.85,0.9,0.99\right\\}$. We repeat each
experiment 100 times. For each synthetic dataset, we use the MCMC algorithm
for Dirichlet process mixture models implemented in the R package `PreMiuM` of
Liverani et al., (2015) to obtain the PSMs. We use discrete mixtures (Liverani
et al.,, 2015, Section 3.2) except in one setting (detailed below) where the
profile regression model of Molitor et al., (2010) is used, using a discrete
mixture with categorical response. The idea of profile regression is that, if
a response $y_{n}$ is available for each $n=1,\dots,N$, the observations
$d_{n}=(x_{n},y_{n})$ are jointly modelled as the product of the response
model and a covariate model. In both cases we use the default hyperparameters,
which we found to work well in practice.

We consider four different simulation settings:

  1. (A)

The clustering structure in every dataset is the same and is related to the
outcome of interest.

  2. (B)

As in setting A, the clustering structure in each dataset is the same. In this
case, however, each dataset contains some additional covariates that have no
clustering structure.

  3. (C)

The dataset with highest cluster separability has a clustering structure that
is unrelated to the response variable, all the other datasets are the same as
in setting A.

  4. (D)

This is the same as setting C, but profile regression is used to derive the
PSMs.

One set of PSMs used for setting A is shown in Figure 3.

Figure 3: PSMs of the datasets used for setting A. The rows and columns
correspond to the statistical units. The coloured bar on the right of each PSM
represents the true clusters. Higher probabilities of belonging to the same
cluster are indicated in blue. The values of $w$ used to generate these
matrices are, from left to right, 0.2, 0.4, 0.6, and 0.8.

In this case we know that higher values of $w_{d}$ are associated with higher
levels of cluster separability. However, in general, we do not know how
dispersed the elements in each matrix are. So, we define another way to score
how strong the signal is in each dataset. We use the _cophenetic correlation
coefficient_ , a measure of how faithfully hierarchical clustering would
preserve the pairwise distances between the original data points. Given a
dataset $X=[\boldsymbol{x}_{1},\boldsymbol{x}_{2},\dots,\boldsymbol{x}_{N}]$
and a similarity matrix $\Delta\in\mathbb{R}^{N\times N}$, we define the
_dendrogrammatic distance_ between $\boldsymbol{x}_{i}$ and
$\boldsymbol{x}_{j}$ as the height of dendrogram at which these two points are
first joined together by hierarchical clustering and we denote it by
$\eta_{ij}$. The cophenetic correlation coefficient $\rho$ is calculated as

$\rho=\frac{\sum_{i<j}(\Delta_{ij}-\bar{\Delta})(\eta_{ij}-\bar{\eta})}{\sqrt{\sum_{i<j}(\Delta_{ij}-\bar{\Delta})\sum_{i<j}(\eta_{ij}-\bar{\eta})}},$
(11lvzabaj)

where $\bar{\Delta}$ and $\bar{\eta}$ are the average values of $\Delta_{ij}$
and $\eta_{ij}$ respectively.

###  4.2 Summarising PSMs

Before using the PSMs with the MKL-based integrative clustering methods, we
carry out some simulation studies to ensure that the kernel $k$-means
algorithm performs equally well at summarising the MCMC output as other
methods from the literature. The advantage is that the kernel $k$-means can be
extended to combine multiple datasets.

We use the adjusted Rand index (ARI; Hubert and Arabie,, 1985) as a measure of
the similarity between the output of each clustering method and the true
partition of the data.

To compute the ARI given two partitions $U=\\{U_{1},\dots,U_{K}\\}$ and
$V=\\{V_{1},\dots,V_{L}\\}$, we summarise the overlapping between each pair of
subsets $U_{i}$ and $V_{j}$ by the values $\nu_{ij}=|U_{i}\cap V_{j}|$,
$i=1,\ldots,K$, $j=1,\dots,L$. The ARI is then calculated as

$ARI=\frac{\sum_{k,l}\binom{\nu_{kl}}{2}-\sum_{k}\binom{\nu_{k\cdot}}{2}\sum_{l}\binom{\nu_{\cdot
l}}{2}\big{/}\binom{\nu}{2}}{\frac{1}{2}\big{[}\sum_{k}\binom{\nu_{k\cdot}}{2}+\sum_{l}\binom{\nu_{\cdot
l}}{2}\big{]}-\sum_{k}\binom{\nu_{k\cdot}}{2}\sum_{l}\binom{\nu_{\cdot
l}}{2}\big{/}\frac{\nu}{2}}.$ (11lvzabak)

This is the corrected-for-chance version of the Rand index (Rand,, 1971) that
is simply the number $n_{s}$ of pairs of elements that are in the same subsets
in both partitions $U$ and $V$, plus the number $n_{d}$ of pairs of elements
that are in different subsets in both $U$ and $V$, divided by the total number
of pairs $N^{2}$

$RI=\frac{n_{s}+n_{d}}{N^{2}}.$ (11lvzabal)

We use the synthetic datasets described in setting A of Section 4.1 and
compare the kernel $k$-means algorithm to the methods implemented in the R
package `mcclust` (Fritsch and Ickstadt,, 2009). All these methods take a PSM
$\Delta$ as input and find the clustering $\boldsymbol{c}^{*}$ that maximises
the posterior expected adjusted Rand index (PEAR). This is achieved by
choosing the clustering $\boldsymbol{c}^{*}$ that maximises the following
quantity:

$\frac{\sum_{i<j}\mathbb{I}_{\\{c^{*}_{i}=c^{*}_{j}\\}}\Delta_{ij}-\sum_{i<j}\mathbb{I}_{\\{c_{i}^{*}=c_{j}^{*}\\}}\sum_{i<j}\Delta_{ij}\big{/}\binom{n}{2}}{\frac{1}{2}\big{[}\sum_{i<j}\mathbb{I}_{\\{c_{i}^{*}=c_{j}^{*}\\}}+\sum_{i<j}\Delta_{ij}\big{]}-\sum_{i<j}\mathbb{I}_{\\{c_{i}^{*}=c_{j}^{*}\\}}\sum_{i<j}\Delta_{ij}\big{/}\frac{n}{2}}.$
(11lvzabam)

The clusterings $\boldsymbol{c}^{*}$ taken into consideration by these methods
can be chosen in different ways. We use hierarchical clustering with
$1-\Delta$ as distance matrix with average and complete linkage. The maxpear
function tries all the possible number of clusters between one and a maximum
number of cluster $K_{\max}$ specified by the user. We consider both 6 and 20
as values for $K_{\max}$. We also use the maxpear function to try with all the
clusterings in the MCMC output and take the one that maximises the quantity of
Equation (11lvzabam). We repeat this procedure using the minVI function of the
R package mcclust.ext of Wade and Ghahramani, (2018), where the selected
$\boldsymbol{c^{*}}$ is the one that minimises the lower bound to the
posterior expected variation of information (Equation 6) from Jensen’s
inequality.

Figure 4 shows the boxplots of the adjusted Rand index obtained repeating the
experiment with 100 different sets of synthetic data calculated for each of
the considered clustering algorithms and for all values of $w$. We can see
that the kernel $k$-means on average performs better than the other methods
when the number of clusters is known and that it performs similarly to the
others when the number of clusters is chosen to maximise the average
silhouette.

Figure 4: ARI for the kernel $k$-means applied to one dataset at a time, for
different values of $\rho$ compared to maximising the PEAR as suggested by
Fritsch and Ickstadt, (2009) and to minimising the VI as suggested by Wade and
Ghahramani, (2018). For both methods, we try different settings, namely:
performing hierarchical clustering on the matrix $1-\Delta$ with average (avg)
and complete (comp) linkage, with maximum number of clusters equal to either 6
or 20, as well as considering all the clusterings samples that are appear in
the MCMC output (draws). For kernel $k$-means, the results obtained fixing the
number of clusters to six and choosing it via the silhouette are presented.

###  4.3 Integrative clustering

We assess the MKL-based integrative clustering approaches described in
Sections 3.3.2 and 3.3.3 in the four settings presented in Section 4.1. For
each setting we consider four different subsets of data, each combining three
out of our four synthetic datasets generated with $w=0.2,0.4,0.6,0.8$. In what
follows, we refer to the dataset generated with $w=0.2$ as “dataset 0”, the
one generated with $w=0.4$ as “dataset 1”, and so on, for ease of reference.
Moreover, we indicate by “0+1+2” the integration of the datasets generated
with values of $w$ equal to 0.2, 0.4, and 0.6 respectively, and similarly for
the other combinations of datasets. Here we show the ARI between the
clusterings found via MKL integration and the true cluster labels, the weights
assigned to each dataset in each setting are instead reported in the
Supplementary Material.

##### Setting A

The ARI obtained by combining the datasets in the unsupervised and outcome-
guided frameworks is shown in the first row of Figure 5. The values of the ARI
obtained in the previous section on each dataset separately are also reported.
In all settings we set the number of clusters to the true value, six. The
unsupervised integration performed using localised multiple kernel $k$-means
allows us to reach values of the ARI that are close to those of the “best”
dataset (i.e. the dataset that has the highest value of cluster separability)
among the three datasets in each subset. This is because the unsupervised MKL
approach considered here assign higher weights to the datasets that give rise
to kernels with higher values of $\rho$, which in this case correspond to
higher values of $w$ (Cabassi and Kirk,, 2020). Moreover, even higher values
of the ARI are achieved via outcome-guided integration thanks to the smarter
weighting of the datasets. In this case, the kernels that help separate the
classes in the response have higher weights than the others.

##### Setting B

In the second row of Figure 5 are shown the results obtained for setting B,
where the PSMs are obtained exploiting (an adaptation of) the variable
selection strategy of Chung and Dunson, (2009) implemented in the R package
PReMiuM. Despite the fact that the ARI of dataset 2 is lower than in the
previous case, the integration results are better than in Setting A. Again,
this is due to the fact that most informative kernels are weighted more highly
than the other ones.

##### Setting C

This simulation study helps us to show that the outcome-guided approach
favours the clustering structures that agree with the structure in the
response. For this reason we use a dataset with high cophenetic correlation
coefficient whose clustering structure is not related to the response. The
results are presented in the third row of Figure 5. Again, localised multiple
kernel $k$-means assigns higher weights to the datasets that are more easily
separable, i.e. datasets that give rise to kernels having higher cophenetic
correlation coefficients. Note that here higher values of $w$ correspond to
higher cophenetic correlation. In this situation, this causes the ARI of the
subsets of kernels that include dataset 4 to drop to zero. In the outcome-
guided case, instead, the dataset that has the highest level of cluster
separability but is not related to the outcome of interest has (almost) always
weight equal to zero.

##### Setting D

Lastly, we consider the case where the model used to generate the PSMs is
profile regression (fourth row of Figure 5). We see that, as expected, the ARI
is higher than in the previous cases for the clustering obtained with each
dataset taken separately, except of course for dataset 4, that has a different
clustering structure. This is reflected in an improvement of the ARI of the
unsupervised and outcome-guided integration, for all considered subsets of
data. In particular, the latter almost always allows us to retrieve the true
clustering.

Figure 5: Simulation study, adjusted Rand index obtained by summarising the
PSMs one at a time using kernel $k$-means (left), combining different subsets
of three PSMs in an unsupervised fashion using localised multiple kernel
$k$-means (centre), and combining the same subsets making use of a response
variable and multi-class SVMs to determine each PSM’s weight and using kernel
$k$-means for the final clustering (right).

##  5 Integrative clustering: biological applications

We present results from applying our method to two exemplar integrative
clustering problems from the literature. In Section 5.1 we perform integrative
clustering on the dataset used for the analysis of 3,527 tumour samples from
12 different tumour types of Hoadley et al., (2014). In Section 5.2 we
consider an example of transcriptional module discovery, to which a number of
existing integrative clustering approaches have previously been applied
(Savage et al.,, 2010, Kirk et al.,, 2012, Cabassi and Kirk,, 2020).

###  5.1 Multiplatform analysis of ten cancer types

We analyse the dataset of Hoadley et al., (2014), which contains data
regarding 3,527 tumour samples spanning 12 different tumour types. This
dataset is particularly suitable for two purposes: (i) determining whether
’omic signatures are shared across different tumour types and (ii) discovering
tumour subtypes. The ’omic layers available are: DNA copy number, DNA
methylation, mRNA expression, microRNA expression and protein expression.
Hoadley et al., (2014) used _Cluster-Of-Clusters Analysis_ (COCA; Cabassi and
Kirk,, 2020) to cluster the data and found that the clusters where highly
correlated with the tissue of origin of each tumour sample and were shown to
be of clinical interest.

Here, we combine the data layers both in the unsupervised and outcome-guided
frameworks. We make use of the C implementation (Mason et al.,, 2016) of the
_multiple dataset integration_ (MDI) method of Kirk et al., (2012) to produce
PSMs for each data layer separately. In order to be able to do so, we only
include in our analysis the tumour samples that have no missing values; this
reduces the sample size to 2,421 and the number of tumour types available for
the analysis to ten. A mixture of Gaussians is used for the continuous layers
(DNA copy number, microRNA, and protein expression), while the multinomial
model is used for the methylation data, which are categorical. Due to the high
number of features, it is not possible to produce a PSM for the full mRNA
dataset, so we exclude it from the analysis presented here. In the
Supplementary Material, however, we show how the variable selection method
developed specifically for multi-omic data by Cabassi et al., (2020) can be
employed in this case to reduce the size of each data layer and integrate all
five data types.

####  5.1.1 Unsupervised integration

We combine the PSMs of the four data layers via multiple kernel $k$-means with
number of clusters going from 2 to 50. We choose the number of clusters that
maximises the silhouette, which is 9 (Supplementary Material). The resulting
clusters are shown in Figure 6. Six out of the nine clusters contain almost
exclusively samples from one tissue: most samples of renal cell carcinoma
(KIRC) are in cluster 2, almost all statistical units in clusters 1 and 4 are
breast cancer samples, most serous ovarian carcinoma (OV) samples are in
cluster 3, bladder urothelial adenocarcinoma (BLCA) samples in cluster 8, and
endometrial cancer (UCEC) samples in cluster 9. Cluster 5, instead, is formed
by the colon and rectal adenocarcinoma samples together, and corresponds
exactly to cluster 7 of Hoadley _et al._ (COAD/READ). Moreover, lung squamous
cell carcinoma (LUSC), lung adenocarcinoma (LUAD), and head and neck squamous
cell carcinoma (HNSC) are divided into two clusters (6 and 7). Cluster 8
contains the remaining samples. The average weights assigned to each data
layer are: 6.9% to the copy number data, 7.5% to the methylation data, 7% to
the microRNA expression data, and 78.7% to the protein expression data.

(a) Clusters and weighted kernel.

(b) Coincidence matrix.

Figure 6: Unsupervised multiplatform analysis of ten cancer types. (a) Left:
weighted kernel. The rows and columns correspond to cancer samples. Higher
values of similarity between samples are indicated in blue. Right: final
clusters, tissues of origin, and COCA clusters. (b) Coincidence matrix
comparing the tissue of origin of the tumour samples (rows) with the clusters
(columns).

####  5.1.2 Outcome-guided integration

We obtain the weights for the outcome-guided integration via the SimpleMKL
algorithm, which are as follows: DNA copy number 35.9%, methylation 13.5%,
microRNA expression 33.8%, and protein expression 16.8%. We then cluster the
data using kernel $k$-means with number of clusters going from 2 to 50. The
silhouette is maximised at $K=27$ (Supplementary Material). The clusters
obtained in this way are shown in Figure 7. It is interesting to note that, in
this case, each cluster contains almost exclusively tumour samples from the
same tissue. The only exceptions are clusters 4 and 22, which contain both
lung and head/neck squamous cell carcinoma samples, and clusters 14 and 25 in
which colon and rectal adenocarcinomas are clustered together, like in the
unsupervised case. Each tumour type, except for ovarian and bladder cancers,
is divided into multiple subclusters. Further analysis would be required to
assess whether these clusters are clinically relevant. Interestingly, we
observe a distinction between luminal (i.e. estrogen receptor-positive and
HER2-positive) and basal breast cancer samples (the former are in clusters 8,
9, 17, 19, 23, 24, 27, while the latter are in cluster 13). This was also
observed by Hoadley et al., (2014).

(a) Clusters and weighted kernel.

(b) Coincidence matrix.

Figure 7: Outcome-guided multiplatform analysis of ten cancer types. (a) Left:
weighted kernel. The rows and columns correspond to cancer samples. Higher
values of similarity between samples are indicated in blue. Right: final
clusters, tissues of origin, and COCA clusters. (b) Coincidence matrix
comparing the tissue of origin of the tumour samples (rows) with the clusters
(columns).

###  5.2 Transcriptional module discovery

For this example we consider transcriptional module discovery for yeast
(_Saccharomyces cerevisiae_). The goal is to find clusters of genes that share
a common biological function and are co-regulated by a common set of
transcription factors. Previous studies have demonstrated that combining gene
expression datasets with information about transcription factor binding can
improve detection of meaningful transcriptional modules (Ihmels et al.,, 2002,
Savage et al.,, 2010).

We combine the ChIP-chip dataset of Harbison et al., (2004), which provides
binding information for 117 transcriptional regulators, with the expression
dataset of Granovskaia et al., (2010). The ChIP-chip data are discretised as
in Savage et al., (2010) and Kirk et al., (2012). The measurements in the
dataset of Granovskaia et al., (2010) represent the expression profiles of 551
genes at 41 different time points of the cell cycle.

Since the goal is to find clusters of genes, here the statistical units
correspond to the genes and the covariates to 41 experiments of gene
expression measurement for the expression dataset and to the 117 considered
transcriptional regulators in the ChIP-chip dataset. To produce the posterior
similarity matrices for the two datasets, we use the `DPMSysBio` Matlab
package of Žurauskienė et al., (2016). For each dataset, we run 10.000
iterations of the MCMC algorithm and summarise the output into a PSM. The PSMs
obtained in this way are reported in the Supplementary Material.

(a) Expression data.

(b) ChIP-chip data.

(c) Weighted kernel.

Figure 8: Transcriptional module discovery, integration of the Harbison et
al., (2004) and Granovskaia et al., (2010) datasets. (a) Expression data. Each
row corresponds to a gene and each column to a different time point. (b) ChIP-
chip data. Each row corresponds to a gene and each column to a transcriptional
regulator. (c) Weighted kernel. The rows and columns correspond to the genes.
Higher values of similarity between genes are indicated in blue. To the left
of each plot is shown the final clustering, obtained by integrating the PSMs
of the expression and ChIP-chip data via multiple kernel $k$-means.

We combine the PSMs using our unsupervised MKL approach. The average weights
assigned by the localised multiple kernel $k$-means to each matrix and the
values of the average silhouette for different numbers of clusters are
reported in the Supplementary Material. We set the number of clusters to 25,
which is the value that maximises the silhouette. The final clusters are shown
in Figure 8 next to the two datasets and the combined PSM. In order to
determine whether our clustering is biologically meaningful, we use the Gene
Ontology Term Overlap (GOTO) scores defined by Mistry and Pavlidis, (2008).
Denoting by $\text{annot}g_{i}$ the set of all direct annotations for each
gene and all of their associated parent terms, the GOTO similarity between two
genes $g_{i},g_{j}$ is the number of annotations that the two genes share:

$\text{sim}_{\text{GOTO}}(g_{i},g_{j})=|\text{annot}g_{i}\cap\text{annot}g_{j}|.$
(11lvzaban)

Details on how to compute the average overall GOTO score of a clustering are
given in the Supplementary Material of Kirk et al., (2012). The GOTO scores
are reported in Table 1. The two datasets combined achieve higher GOTO scores
than those of the clusters obtained using each dataset separately. We also
compare these GOTO scores to those obtained with the two methods used in
Cabassi and Kirk, (2020). The first one is COCA, a simple, unweighted
algorithm for integrative clustering that is widely used in practice (The
Cancer Genome Atlas Research Network,, 2012, Hoadley et al.,, 2014). The other
integrative method considered here is KLIC (_Kernel Learning Integrative
Clustering_), by which we mean the integration of multiple kernels via
multiple kernel $k$-means, where the kernels are generated via consensus
clustering (Figure 1). The clusters obtained with these two alternative
methods have lower GOTO scores than the integration of PSMs. For COCA, this is
result not unexpected, since the method is unweighted and has previously been
shown to perform less well than KLIC. The difference between what is referred
to as KLIC here and the unsupervised integration of PSMs, however, only lies
in how the kernels are constructed. These scores therefore suggest that
kernels generated from probabilistic models can lead to more accurate results
than those built using consensus clustering.

Dataset(s) | GOTO BP | GOTO MF | GOTO CC
---|---|---|---
ChIP data (Harbison et al.) | 6.18 | 0.97 | 8.54
Expression data (Granovskaia et al.) | 7.07 | 1.04 | 8.90
ChIP+Expression data: COCA | 5.74 | 0.90 | 8.19
ChIP+Expression data: KLIC | 6.60 | 0.96 | 8.66
ChIP+Expression data: integration of PSMs | 7.15 | 1.05 | 8.93
Table 1: Gene Ontology Term Overlap scores. “BP” stands for Biological Process
ontology, “MF” for Molecular Function, and “CC” for Cellular Component. The
number of clusters used for every method is 25.

##  6 Conclusion

We have presented a novel method for summarising a sample of clusterings from
the posterior distribution of an MCMC algorithm for Bayesian clustering, based
on kernel methods. We have also extended this method to allow us to integrate
multiple PSMs. This can be done either in an unsupervised or in an outcome-
guided way. The former weights each PSM according to how well defined is the
clustering structure that it shows, the latter gives more importance to the
PSMs that better reflect the structure encountered in the response variable of
choice.

We have used simulation examples to show that our method gives comparable
performances in terms of proportion of correct co-clustering as the existing
techniques. We have also demonstrated that the integration of multiple
datasets gives better results than using one dataset at a time. Additionally,
we proved that, if a variable related to the output of interest is available,
our method can assign higher weights to the PSMs that are more closely related
to that. The simulation examples prove that this feature can be extremely
useful when not all the PSMs have the same clustering structure.

Finally, we have applied the novel methods to two real data applications. The
pancancer data analysis shows that the outcome-guided integration of multiple
PSMs can potentially be used in the context of tumour subtype discovery. The
yeast example demonstrates that the proposed method is able to identify groups
of genes that are co-expressed and co-regulated that are more biologically
meaningful than those determined via state-of-the-art integrative algorithms.

## Funding

A. Cabassi and P.D.W. Kirk are supported by the MRC [MC_UU_00002/13]. This
work was supported by the National Institute for Health Research [Cambridge
Biomedical Research Centre at the Cambridge University Hospitals NHS
Foundation Trust]. The views expressed are those of the authors and not
necessarily those of the NHS, the NIHR or the Department of Health and Social
Care. Partly funded by the RESCUER project. RESCUER has received funding from
the European Union’s Horizon 2020 research and innovation programme under
grant agreement No. 847912.

## References

  * Antoniak, (1974) Antoniak, C. E. (1974).  Mixtures of dirichlet processes with applications to bayesian nonparametric problems.  The annals of statistics, pages 1152–1174.
  * Baudat and Anouar, (2000) Baudat, G. and Anouar, F. (2000).  Generalized discriminant analysis using a kernel approach.  Neural Computation, 12(10):2385–2404.
  * Bennett and Mangasarian, (1992) Bennett, K. and Mangasarian, O. L. (1992).  Robust linear programming discrimination of two linearly inseparable sets.  Optimization Methods and Software, 1(1):23–34.
  * Binder, (1978) Binder, D. A. (1978).  Bayesian cluster analysis.  Biometrika, 1:31–38.
  * Bishop, (2006) Bishop, C. M. (2006).  Pattern Recognition and Machine Learning.  Springer, 1st edition edition.
  * Blei et al., (2006) Blei, D. M., Jordan, M. I., et al. (2006).  Variational inference for dirichlet process mixtures.  Bayesian analysis, 1(1):121–143.
  * Boser et al., (1992) Boser, B. E., Guyon, I. M., and Vapnik, V. N. (1992).  A training algorithm for optimal margin classifiers.  Proceedings of the 5th annual workshop on Computational Learning Theory.
  * Cabassi and Kirk, (2020) Cabassi, A. and Kirk, P. D. W. (2020).  Multiple kernel learning for integrative consensus clustering of ’omic datasets.  Bioinformatics.  btaa593.
  * Cabassi et al., (2020) Cabassi, A., Seyres, D., Frontini, M., and Kirk, P. D. W. (2020).  Penalised logistic regression for multi-omic data with an application to cardiometabolic syndrome.  arXiv.  2008.00235.
  * Celeux et al., (2000) Celeux, G., Hurn, M., and Robert, C. P. (2000).  Computational and inferential difficulties with mixture posterior distributions.  Journal of the American Statistical Association, 95(451):957–970.
  * Chung and Dunson, (2009) Chung, Y. and Dunson, D. B. (2009).  Nonparametric bayes conditional distribution modeling with variable selection.  Journal of the American Statistical Association, 104(488):1646–1660.
  * Dahl, (2006) Dahl, D. B. (2006).  Model-based clustering for expression data via a dirichlet process mixture model.  Bayesian inference for gene expression and proteomics, 4:201–218.
  * Dudoit and Fridlyand, (2002) Dudoit, S. and Fridlyand, J. (2002).  A prediction-based resampling method for estimating the number of clusters in a dataset.  Genome Biology, 3(7).
  * Everitt, (1993) Everitt, B. S. (1993).  Cluster analysis.  Hodder Education.
  * Ferguson, (1973) Ferguson, T. S. (1973).  A bayesian analysis of some nonparametric problems.  The Annals of Statistics, pages 209–230.
  * Friedman et al., (2001) Friedman, J., Hastie, T., and Tibshirani, R. (2001).  The elements of statistical learning, volume 1.  Springer Series in Statistics.
  * Fritsch and Ickstadt, (2009) Fritsch, A. and Ickstadt, K. (2009).  Improved Criteria for Clustering Based on the Posterior Similarity Matrix.  Bayesian Analysis, 4(2):367–392.
  * Gelfand and Smith, (1990) Gelfand, A. E. and Smith, A. F. (1990).  Sampling-based approaches to calculating marginal densities.  Journal of the American statistical association, 85(410):398–409.
  * Girolami, (2002) Girolami, M. (2002).  Mercer kernel-based clustering in feature space.  IEEE Transactions on Neural Networks, 13(3):780–784.
  * Gönen and Margolin, (2014) Gönen, M. and Margolin, A. A. (2014).  Localized data fusion for kernel k-means clustering with application to cancer biology.  In Advances in Neural Information Processing Systems, pages 1305–1313.
  * Granovskaia et al., (2010) Granovskaia, M. V., Jensen, L. J., Ritchie, M. E., Toedling, J., Ning, Y., Bork, P., Huber, W., and Steinmetz, L. M. (2010).  High-resolution transcription atlas of the mitotic cell cycle in budding yeast.  Genome Biology, 11(3):R24.
  * Harbison et al., (2004) Harbison, C. T., Gordon, D. B., Lee, T. I., Rinaldi, N. J., Macisaac, K. D., Danford, T. W., Hannett, N. M., Tagne, J.-B., Reynolds, D. B., Yoo, J., et al. (2004).  Transcriptional regulatory code of a eucaryotic genome.  Nature, 431(7004):99–104.
  * Hartigan and Wong, (1979) Hartigan, J. A. and Wong, M. A. (1979).  Algorithm as 136: A k-means clustering algorithm.  Journal of the Royal Statistical Society. Series c (Applied Statistics), 28(1):100–108.
  * Hoadley et al., (2014) Hoadley, K. A., Yau, C., Wolf, D. M., Cherniack, A. D., Tamborero, D., Ng, S., Leiserson, M. D., Niu, B., McLellan, M. D., Uzunangelov, V., et al. (2014).  Multiplatform analysis of 12 cancer types reveals molecular classification within and across tissues of origin.  Cell, 158(4):929–944.
  * Hubert and Arabie, (1985) Hubert, L. and Arabie, P. (1985).  Comparing partitions.  Journal of Classification, 2(1):193–218.
  * Ihmels et al., (2002) Ihmels, J., Friedlander, G., Bergmann, S., Sarig, O., Ziv, Y., and Barkai, N. (2002).  Revealing modular organization in the yeast transcriptional network.  Nature Genetics, 31(4):370.
  * Kaufman and Rousseeuw, (1990) Kaufman, L. and Rousseeuw, P. J. (1990).  Finding groups in data: an introduction to cluster analysis, volume 344.  John Wiley & Sons.
  * Kirk et al., (2012) Kirk, P. D., Griffin, J. E., Savage, R. S., Ghahramani, Z., and Wild, D. L. (2012).  Bayesian correlated clustering to integrate multiple datasets.  Bioinformatics, 28(24):3290–3297.
  * Knerr et al., (1990) Knerr, S., Personnaz, L., and Dreyfus, G. (1990).  Single-layer learning revisited: a stepwise procedure for building and training a neural network.  Neurocomputing: Algorithms, Architectures and Applications, 68(41-50):71.
  * Lanckriet et al., (2004) Lanckriet, G. R., De Bie, T., Cristianini, N., Jordan, M. I., and Noble, W. S. (2004).  A statistical framework for genomic data fusion.  Bioinformatics, 20(16):2626–2635.
  * Liverani et al., (2015) Liverani, S., Hastie, D. I., Azizi, L., Papathomas, M., and Richardson, S. (2015).  PReMiuM: An R package for Profile Regression Mixture Models Using Dirichlet Processes.  Journal of Statistical Software, 31(2).
  * Luenberger and Ye, (1984) Luenberger, D. G. and Ye, Y. (1984).  Linear and nonlinear programming, volume 2.  Springer.
  * Mason et al., (2016) Mason, S. A., Sayyid, F., Kirk, P. D., Starr, C., and Wild, D. L. (2016).  Mdi-gpu: accelerating integrative modelling for genomic-scale data using gp-gpu computing.  Statistical Applications in Genetics and Molecular Biology, 15(1):83–86.
  * McLachlan and Peel, (2004) McLachlan, G. J. and Peel, D. (2004).  Finite mixture models.  John Wiley & Sons.
  * Medvedovic and Sivaganesan, (2002) Medvedovic, M. and Sivaganesan, S. (2002).  Bayesian infinite mixture model based clustering of gene expression profiles.  Bioinformatics, 18(9):1194–1206.
  * Meilă, (2007) Meilă, M. (2007).  Comparing clusterings—an information based distance.  Journal of multivariate analysis, 98(5):873–895.
  * Mika et al., (1999) Mika, S., Rätsch, G., Weston, J., Schölkopf, B., and Müllers, K.-R. (1999).  Fisher discriminant analysis with kernels.  In Proceedings of the 1999 IEEE Signal Processing Society Workshop., pages 41–48. IEEE.
  * Milligan and Cooper, (1985) Milligan, G. W. and Cooper, M. C. (1985).  An examination of procedures for determining the number of clusters in a data set.  Psychometrika, 50(2):159–179.
  * Mistry and Pavlidis, (2008) Mistry, M. and Pavlidis, P. (2008).  Gene ontology term overlap as a measure of gene functional similarity.  BMC bioinformatics, 9:327.
  * Molitor et al., (2010) Molitor, J., Papathomas, M., Jerrett, M., and Richardson, S. (2010).  Bayesian profile regression with an application to the national survey of children’s health.  Biostatistics, 11(3):484–498.
  * Monti et al., (2003) Monti, S., Tamayo, P., Mesirov, J., and Golub, T. (2003).  Consensus clustering: a resampling-based method for class discovery and visualization of gene expression microarray data.  Machine Learning, 52(1-2):91–118.
  * Neal, (2000) Neal, R. M. (2000).  Markov Chain Sampling Methods for Dirichlet Process Mixture Models.  Journal of Computational and Graphical Statistics, 9(2):249–265.
  * Platt, (1999) Platt, J. C. (1999).  Fast training of support vector machines using sequential minimal optimization.  In B. Schölkopf, C. J. C. Burges, A. J. S., editor, Advances in Kernel Methods – Support Vector Learning, chapter 12, pages 185–208. MIT Press.
  * Rakotomamonjy and Bach, (2007) Rakotomamonjy, A. and Bach, F. (2007).  More efficiency in multiple kernel learning.  In Proceedings of the 24th International Conference on Machine Learning, pages 775–782.
  * Rakotomamonjy et al., (2008) Rakotomamonjy, A., Bach, F. R., Canu, S., and Grandvalet, Y. (2008).  SimpleMKL.  Journal of Machine Learning Research, 9:2491–2521.
  * Rand, (1971) Rand, W. M. (1971).  Objective Criteria for the Evaluation of Clustering Methods.  Journal of the American Statistical Association, 66(336):846–850.
  * Rasmussen, (2000) Rasmussen, C. E. (2000).  The infinite gaussian mixture model.  In Advances in neural information processing systems, pages 554–560.
  * Richardson and Green, (1997) Richardson, S. and Green, P. (1997).  On Bayesian Analysis of Mixtures with an Unknown Number of Components.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), (59.4):731–792.
  * Roth and Steinhage, (2000) Roth, V. and Steinhage, V. (2000).  Nonlinear discriminant analysis using kernel functions.  In Advances in neural information processing systems, pages 568–574.
  * Rousseeuw, (1987) Rousseeuw, P. J. (1987).  Silhouettes: a graphical aid to the interpretation and validation of cluster analysis.  Journal of computational and applied mathematics, 20:53–65.
  * Savage et al., (2010) Savage, R. S., Ghahramani, Z., Griffin, J. E., De La Cruz, B. J., and Wild, D. L. (2010).  Discovering transcriptional modules by bayesian data integration.  Bioinformatics, 26(12):i158–i167.
  * Schölkopf et al., (1998) Schölkopf, B., Smola, A., and Müller, K.-R. (1998).  Nonlinear component analysis as a kernel eigenvalue problem.  Neural Computation, 10(5):1299–1319.
  * Schölkopf and Smola, (2001) Schölkopf, B. and Smola, A. J. (2001).  Learning with kernels: support vector machines, regularization, optimization, and beyond.  MIT press.
  * Shawe-Taylor and Cristianini, (2004) Shawe-Taylor, J. and Cristianini, N. (2004).  Kernel methods for pattern analysis.  Cambridge University Press.
  * Steinhaus, (1956) Steinhaus, H. (1956).  Sur la division des corps matériels en parties.  Bulletin de l’Académie Polonaise des Sciences, IV(12):801–804.
  * The Cancer Genome Atlas Research Network, (2012) The Cancer Genome Atlas Research Network (2012).  Comprehensive molecular portraits of human breast tumours.  Nature, 487(7407):61–70.
  * Tibshirani et al., (2001) Tibshirani, R. et al. (2001).  Estimating the number of clusters in a data set via the gap statistic.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63(2):411–423.
  * Vapnik, (1998) Vapnik, V. N. (1998).  Statistical learning theory.  Wiley New York.
  * Wade and Ghahramani, (2018) Wade, S. and Ghahramani, Z. (2018).  Bayesian cluster analysis: Point estimation and credible balls (with discussion).  Bayesian Analysis, 13(2):559–626.
  * Yeung et al., (2001) Yeung, K. Y., Haynor, D. R., and Ruzzo, W. L. (2001).  Validating clustering for gene expression data.  Bioinformatics, 17(4):309–318.
  * Žurauskienė et al., (2016) Žurauskienė, J., Kirk, P. D., and Stumpf, M. P. (2016).  A graph theoretical approach to data fusion.  Statistical Applications in Genetics and Molecular Biology, 15(2):107–122.