Bayesian outcome-guided multi-view mixture models with applications in molecular precision medicine

TL;DR

This paper introduces a multi-view Bayesian mixture model that identifies multiple clustering structures in high-dimensional 'omics data, guided by clinical outcomes, to improve disease subtype discovery in molecular medicine.

Contribution

It proposes a novel semi-supervised, outcome-guided multi-view Bayesian clustering method that captures distinct biological processes and enhances stratified medicine applications.

Findings

Effective in simulated data demonstrating multiple clustering structures.

Applied to pan-cancer proteomics, revealing biologically meaningful subtypes.

Successfully integrated multi-omics data for breast cancer subtyping.

Abstract

Clustering is commonly performed as an initial analysis step for uncovering structure in 'omics datasets, e.g. to discover molecular subtypes of disease. The high-throughput, high-dimensional nature of these datasets means that they provide information on a diverse array of different biomolecular processes and pathways. Different groups of variables (e.g. genes or proteins) will be implicated in different biomolecular processes, and hence undertaking analyses that are limited to identifying just a single clustering partition of the whole dataset is therefore liable to conflate the multiple clustering structures that may arise from these distinct processes. To address this, we propose a multi-view Bayesian mixture model that identifies groups of variables (``views"), each of which defines a distinct clustering structure. We consider applications in stratified medicine, for which our principal goal is to identify clusters of patients that define distinct, clinically actionable disease subtypes. We adopt the semi-supervised, outcome-guided mixture modelling approach of Bayesian profile regression that makes use of a response variable in order to guide inference toward the clusterings that are most relevant in a stratified medicine context. We present the model, together with illustrative simulation examples, and examples from pan-cancer proteomics. We demonstrate how the approach can be used to perform integrative clustering, and consider an example in which different 'omics datasets are integrated in the context of breast cancer subtyping.

Equations

$$p(\mathbf{v}|\mathbf{\rho },\mathbf{\pi })$$

$$=\sum _{k=1}^K{\pi }_k{f}_{\mathbf{y}}(\mathbf{y}|{\mathbf{\theta }}_k,\mathbf{x})f_{\mathbf{x}}(\mathbf{x}|{\mathbf{\varphi }}_k),$$

$$p(\mathbf{x},\mathbf{y}|\mathbf{\Phi },\mathbf{\Theta },\mathbf{\pi },\mathbf{\beta },\mathbf{w})=\sum _{k=1}^K{\pi }_k{f}_{\mathbf{y}}(\mathbf{y}|{\mathbf{\theta }}_k,\mathbf{w},\mathbf{\beta })f_{\mathbf{x}}(\mathbf{x}|{\mathbf{\varphi }}_k).$$

$$p({\mathbf{x}}_i,{\mathbf{y}}_i|z_i=k,\mathbf{\varphi },\mathbf{\theta },\mathbf{\beta },{\mathbf{w}}_i)=f_{\mathbf{y}}({\mathbf{y}}_i|{\mathbf{\theta }}_k,{\mathbf{w}}_i,\mathbf{\beta })f_{\mathbf{x}}({\mathbf{x}}_i|{\mathbf{\varphi }}_k).$$

$$p(z_i=k|D,$$

$$p(z_i=k,\text{mbox}wherek\in {\mathbf{z}}_{-i}|D,$$

$$p(z_i\notin {\mathbf{z}}_{-i}|D,$$

$$p(z_i\notin {\mathbf{z}}_{-i}|D,{\mathbf{z}}_{-i},\dots )+\sum _{k\in {\mathbf{z}}_{-i}}p(z_i=k,\text{mbox}wherek\in {\mathbf{z}}_{-i}|D,{\mathbf{z}}_{-i},\dots )=1.$$

$$f({\mathbf{x}}_i,{\mathbf{y}}_i|z_i=k,\mathbf{\varphi },\mathbf{\theta },\mathbf{\beta },{\mathbf{w}}_i)$$

$$=f_{\mathbf{y}}({\mathbf{y}}_i|{\mathbf{\theta }}_k,\mathbf{\beta },{\mathbf{w}}_i)\prod _{j=1}^P{f}_x(x_{ij}|{\mathbf{\varphi }}_k),\text{vspace}-0.5cm$$

$$f({\mathbf{x}}_i,y_i|z_i=k,\mathbf{\varphi },{\mathbf{\varphi }}_{0\cdot },\mathbf{\theta },{\mathbf{w}}_i,\mathbf{\beta },\mathbf{\gamma })$$

$$=f_{\mathbf{y}}({\mathbf{y}}_i|{\mathbf{\theta }}_k,\mathbf{\beta },{\mathbf{w}}_i)\prod _{j=1}^P{q}_x(x_{ij}|{\mathbf{\varphi }}_{0j}{)}^{\text{mathbbm}{1}_{{\gamma }_j=0}}{f}_x(x_{ij}|{\mathbf{\varphi }}_k{)}^{\text{mathbbm}{1}_{{\gamma }_j=1}}.\text{vspace}-0.5cm$$

$$f({\mathbf{x}}_i,y_i|\{z_{i\ell }=k_{\ell },{\mathbf{\varphi }}_{\ell \cdot }{\}}_{\ell =1}^{L-1},{\mathbf{\varphi }}_{0\cdot },\mathbf{\theta },$$

$$f_{\mathbf{y}}({\mathbf{y}}_i|{\mathbf{\theta }}_{k_1},\mathbf{\beta },{\mathbf{w}}_i)\prod _{j=1}^P{q}_x(x_{ij}|{\mathbf{\varphi }}_{0j}{)}^{\text{mathbbm}{1}_{{\gamma }_j=0}}\prod _{\ell =1}^{L-1}{f}_x(x_{ij}|{\mathbf{\varphi }}_{\ell k_{\ell }}{)}^{\text{mathbbm}{1}_{{\gamma }_j=\ell }}.\text{vspace}-0.5cm$$

$$f({\mathbf{x}}_i,y_i|$$

$$\displaystyle\definecolor{currentcolor}{rgb}{0,0,0}\color[rgb]{0,0,0}\underbrace{\color[rgb]{0,0,0}\left(\prod_{j\in\Gamma_{0}}q_{x}(x_{ij}|\boldsymbol{\phi}_{0j})\right)}_{\color[rgb]{0,0,0}\mathclap{\text{View 0}}}\color[rgb]{0,0,0}\definecolor{currentcolor}{rgb}{0,0,0}\color[rgb]{0,0,0}\underbrace{\color[rgb]{0,0,0}\left(f_{\bf y}({\bf y}_{i}|\boldsymbol{\theta}_{k_{1}},\boldsymbol{\beta},\mathbf{w}_{i})\prod_{j\in\Gamma_{1}}f_{x}(x_{ij}|\boldsymbol{\phi}_{1k_{1}})\right)}_{\color[rgb]{0,0,0}\mathclap{\text{View 1}}}\color[rgb]{0,0,0}\prod_{\ell=2}^{L-1}\definecolor{currentcolor}{rgb}{0,0,0}\color[rgb]{0,0,0}\underbrace{\color[rgb]{0,0,0}\left(\prod_{j\in\Gamma_{\ell}}f_{x}(x_{ij}|\boldsymbol{\phi}_{\ell k_{\ell}})\right)}_{\color[rgb]{0,0,0}\mathclap{\text{View $\ell$}}}\color[rgb]{0,0,0},\vspace{-0.5cm}$$

$$\displaystyle p(\gamma_{j}=\ell|z_{\cdot\ell},\boldsymbol{\phi}_{\ell\cdot})=\frac{1}{Z}\nu_{\ell}\prod_{i=1}^{n}f_{x}(x_{ij}|\boldsymbol{\phi}_{\ell z_{i\ell}})\mbox{ for $\ell=1,\ldots,L-1$},$$

$$p({\gamma }_j=0|{\mathbf{\varphi }}_{0j})=\frac{1}{Z}{\nu }_0\prod _{i=1}^n{q}_x(x_{ij}|{\mathbf{\varphi }}_{0j}),$$

$$Z={\nu }_0\prod _{i=1}^n{q}_x(x_{ij}|{\mathbf{\varphi }}_{0j})+\sum _{\ell =1}^{L-1}\left({\nu }_{\ell }\prod _{i=1}^n{f}_x(x_{ij}|{\mathbf{\varphi }}_{\ell z_{i\ell }})\right).$$

$$x_{ij}|z_i=k,{\gamma }_j=\ell$$

$$\text{mbox}where{\mathbf{\varphi }}_{\ell kj}=[{\varphi }_{\ell kj1},{\varphi }_{\ell kj2},{\varphi }_{\ell kj3}]$$

$$p(\mathbf{x}|\mathbf{\varphi },\mathbf{\pi })$$

$$p({\mathbf{x}}_i|z_i,\mathbf{\varphi })=f_{\mathbf{x}}({\mathbf{x}}_i|{\mathbf{\varphi }}_{z_i}),$$

$$p({\mathbf{x}}_i,z_i=k|\mathbf{\varphi },\pi )$$

$$=f_{\mathbf{x}}({\mathbf{x}}_i|{\mathbf{\varphi }}_k){\pi }_k.$$

$$p({\mathbf{x}}_i|\mathbf{\varphi },\pi )$$

$$p({\mathbf{x}}_i,z_i,\mathbf{\varphi },\pi )$$

$$=f_{\mathbf{x}}({\mathbf{x}}_i|{\mathbf{\varphi }}_{z_i})p(z_i|\mathbf{\pi })p(\mathbf{\pi })\prod _{k=1}^{K}p({\mathbf{\varphi }}_k),$$

$$p({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n,z_1,\dots ,z_n,\mathbf{\varphi },\mathbf{\pi })$$

$$=\left(\prod _{i=1}^n{f}_{\mathbf{x}}({\mathbf{x}}_i|{\mathbf{\varphi }}_{z_i})p(z_i|\mathbf{\pi })\right)p(\mathbf{\pi })\prod _{k=1}^{K}p({\mathbf{\varphi }}_k).$$

$$p({\mathbf{\varphi }}_k|{\mathbf{x}}_1,\dots ,{\mathbf{x}}_n,z_1,\dots ,z_n,{\mathbf{\varphi }}_{-k},{\pi }_1,\dots ,{\pi }_K)\propto p({\mathbf{\varphi }}_k)\prod _{i:z_i=k}{f}_{\mathbf{x}}({\mathbf{x}}_i|{\mathbf{\varphi }}_{z_i}),$$

$$p({\pi }_1,\dots ,{\pi }_K|{\mathbf{x}}_1,\dots ,{\mathbf{x}}_n,z_1,\dots ,z_n,{\mathbf{\varphi }}_1,\dots ,{\mathbf{\varphi }}_K)$$

$$p(z_i=k|{\mathbf{x}}_1,\dots ,{\mathbf{x}}_n,{\mathbf{\varphi }}_1,\dots ,{\mathbf{\varphi }}_K,{\pi }_1,\dots ,{\pi }_K,z_{-i})$$

$$={\pi }_k{f}_{\mathbf{x}}({\mathbf{x}}_i|{\varphi }_k),$$

$$p(z_i=k|{\mathbf{x}}_1,\dots ,{\mathbf{x}}_n,\mathbf{\varphi },\pi ,z_{-i})$$

Full text

1 Introduction

Clustering is ubiquoitously used in the analysis of omics data as a means to uncover structure and patterns in these large, high-dimensional datasets (e.g. Eisen et al., 1998; Heyer et al., 1999; Alon et al., 1999; Ben-Dor et al., 1999; Son et al., 2005). Here we are particularly interested in molecular precision medicine applications, in which the analysis objective is to identify molecular subtypes of disease (e.g. Golub et al., 1999; Perou et al., 2000; Sørlie et al., 2001; Cancer Genome Atlas Network, 2012; Kuijjer et al., 2018). In this context, the aim is to identify clusters of patients on the basis of a diverse range of molecular variables, measurements of which are obtained using high-throughput ’omics technologies. One feature of these datasets is that they are typically high-dimensional, which frequently necessitates the use of variable screening or selection strategies (e.g. Witten and Tibshirani, 2010; Fop and Murphy, 2018; Crook et al., 2018), or other dimension reduction techniques (e.g. Yeung and Ruzzo, 2001; McLachlan et al., 2002; Taschler et al., 2019). However, a commonly overlooked challenge is that ’omics datasets often define multiple clustering structures in the patient population, as a consequence of different subsets of variables (e.g. those corresponding to functional groups of genes or proteins) being implicated in a variety of different biomolecular processes. Thus, depending on the subset of variables we consider, we can identify different patient clusters.

A number of papers have proposed methods for identifying multiple clustering structures (e.g. Cui et al., 2007; Niu et al., 2010; Guan et al., 2010; Li and Shafto, 2011; Niu et al., 2014). In the literature, a set of variables that define the same clustering structure has been termed a view, while the task of identifying views and their associated clustering structures has been referred to as either multi-view clustering (Cui et al., 2007) or cross-clustering (Li and Shafto, 2011). One potential challenge faced by these approaches is how to decide which of the identified clustering structures is the most useful or relevant for a given task. This is particularly important in stratified medicine and disease subtyping applications, where we seek clusters (strata) that define groups of patients who have, for example, similar prognoses or respond similarly to treatment (e.g. Perou et al., 2000; Sørlie et al., 2001). In practice, to determine if a given clustering structure is “relevant”, it is common to make use of a left out outcome variable $y$ (e.g. survival data) and to assess whether or not different clusters are associated with different distributions of $y$ (e.g. Curtis et al., 2012). An alternative approach is to adopt a semi-supervised (outcome guided) approach that makes use of $y$ when performing the clustering analysis.

Bayesian profile regression is one such semi-supervised mixture modelling approach that makes use of an outcome/response in order to guide inference toward relevant clustering structures (Molitor et al., 2010). Informally, a clustering is said to be relevant (for a given response) if individuals allocated to the same cluster tend to have similar values for the response; or, more generally, if the responses of individuals in the same cluster can be accurately described by a common model. More precisely, a clustering is defined to be relevant (for a given response) if the value taken by an individual’s response is not independent of their cluster allocation. It is clear from this definition that the relevance of a clustering can only be specified relative to a given response – and, in particular, that for different responses, different clusterings might be relevant.

For example, we could use Bayesian profile regression to retrospectively cluster patients on the basis of genetic or genomics data, using their survival times as a response to guide the clustering toward prognostically relevant disease subtypes. If we were to use a different response (e.g. height), we might end up with a completely different clustering structure. Neither one of these clustering structures would necessarily be “wrong” – they are just relevant with respect to different responses. Crucially, if we were to adopt an unsupervised clustering approach (which is commonly the default analysis choice), there is no guarantee that this would identify a clustering structure that was relevant for either response.

As we demonstrate in Section 5.1, a limitation of the Bayesian profile regression model is that, with increasing data dimension, the influence exerted by the (typically low dimensional) response on the inference of the clustering structure grows weaker. Here we propose a semi-supervised multi-view Bayesian clustering model that simultaneously addresses both this limitation, as well as the challenge faced by existing multi-view approaches of picking out the (most) relevant clustering structure.

2 Profile regression

We suppose that we have data comprising observations on a vector of clustering variables, ${\bf x}$, and responses, ${\bf y}$. We denote the concatenated vector of clustering variables and response by ${\bf v}=[{\bf x};{\bf y}]$, and model the data using a mixture model with $K$ components (where $K$ could be finite or infinite), as follows:

[TABLE]

where $\pi_{k}$ is the mixture weight associated with the $k$-th component, $\boldsymbol{\rho}_{k}$ denotes the parameters associated with the $k$-th component, and we write $\boldsymbol{\rho}$ and $\boldsymbol{\pi}$ to denote $\{\boldsymbol{\rho}_{k}\}_{k=1}^{K}$ and $\{\pi_{k}\}_{k=1}^{K}$ respectively. We assume that the joint density $f_{\bf v}$ can be factorised into $f_{\bf x}$ and $f_{{\bf y}}$ as shown in Equation (2), with $\boldsymbol{\phi}_{k}$ denoting the parameters of the model for ${\bf x}$, and $\boldsymbol{\theta}_{k}$ denoting the parameters of the model for ${\bf y}$. We will write $\boldsymbol{\Phi}$ and $\boldsymbol{\Theta}$ to denote $\{\boldsymbol{\phi}_{k}\}_{k=1}^{K}$ and $\{\boldsymbol{\theta}_{k}\}_{k=1}^{K}$ respectively. In the profile regression model, it is further assumed that $f_{\bf y}({\bf y}|\boldsymbol{\theta}_{k},{\bf x})=f_{\bf y}({\bf y}|\boldsymbol{\theta}_{k})$; i.e. that ${\bf y}$ is conditionally independent of ${\bf x}$ given $\boldsymbol{\theta}_{k}$ (Molitor et al., 2010). A related model, in which this conditional independence assumption is not made, is given in Shahbaba and Neal (2009).

In Molitor et al. (2010), the authors allow the model for ${\bf y}$ to include a dependence upon additional adjustment covariates, ${\bf w}$, that may be predictive of the response but that we do not wish to contribute to the clustering (e.g. confounders that we wish to control for, such as age or sex), together with associated “global” (i.e. not component-specific) parameters $\boldsymbol{\beta}$ (see Appendix for details). The general profile regression model is then:

[TABLE]

The original formulation of the profile regression model, which we also adopt here, is specifically in terms of infinite mixture models using Dirichlet process priors (Molitor et al., 2010); however, we note that the model is equally applicable in the case of finite $K$.

2.1 Dirichlet process formulation

In Molitor et al. (2010), a finite approximation to the Bayesian nonparametric case was considered, using a truncated stick breaking construction of the Dirichlet process (Ishwaran and James, 2001) to define the prior on the $\pi_{k}$’s in Equation (3). Subsequent Bayesian profile regression papers (Hastie et al., 2014) and implementations (Liverani et al., 2015) also considered stick breaking constructions, with inference performed via slice sampling (Walker, 2007; Kalli et al., 2011). Following the derivations of Neal (2000) and Rasmussen (2000) in the unsupervised case, here we instead consider the Dirichlet process mixture model as a limiting case of a (finite) $K$ component mixture model, in which a symmetric Dirichlet prior with parameter $\alpha/K$ is placed on the mixture weights, $\pi_{k}$ (see also Ishwaran and Zarepour, 2002). This is closely related to the Pólya urn (Blackwell and MacQueen, 1973) and Chinese restaurant process (Aldous et al., 1985) constructions for the Dirichlet process, and permits inference via a collapsed Gibbs sampler in which the mixture weights are marginalised. Details are provided in the Appendix, with sampling performed as in Neal (2000) – although we note that alternative approaches for performing inference in Bayesian mixture models could also be employed in this context (e.g. Richardson and Green, 1997; Green and Richardson, 2001; Jain and Neal, 2004, 2007; Walker, 2007; Kalli et al., 2011; Miller and Harrison, 2017)

To provide a very brief overview, let $D=\{({\bf x}_{i},{\bf y}_{i},{\bf w}_{i})\}_{i=1}^{n}$ denote a dataset comprising $({\bf x},{\bf y},{\bf w})$ triples for $n$ individuals, such that $({\bf x}_{i},{\bf y}_{i},{\bf w}_{i})$ corresponds to the $i$-th individual. As is common for mixture models, we introduce latent component allocation variables, $z_{i}$, where $z_{i}=k$ if the $i$-th individual is associated with the $k$-th component, and $p(z_{i}=k|\boldsymbol{\pi})=\pi_{k}$. We define ${\bf z}=\{z_{1},\ldots,z_{n}\}$ to be the multiset of all $n$ component allocations. The component-conditional likelihood associated with the $i$-th individual is then:

[TABLE]

Independent priors, $p(\boldsymbol{\theta}),p(\boldsymbol{\phi})$ and $p(\boldsymbol{\beta})$, are taken for the component-specific and global parameters. The collapsed Gibbs sampler then iterates between:

• •

Step (i): updating the component allocations, $z_{i}$, given the data $D$, the most recently sampled parameters $\boldsymbol{\Theta},\boldsymbol{\Phi},\boldsymbol{\beta}$, and the most recently sampled values for the other allocation variables, ${\bf z}_{-i}={\bf z}\backslash\{z_{i}\}$; and

• •

Step (ii): updating the parameters $\boldsymbol{\Theta},\boldsymbol{\Phi},\boldsymbol{\beta}$, given the data $D$ and the most recently sampled component allocations ${\bf z}$.

In the finite $K$ case, if we take a symmetric Dirichlet prior with parameter $\alpha/K$ for the mixture weights, $\pi_{1},\ldots,\pi_{K}\sim\mbox{Dir}(\alpha/K)$, then the conditional posterior probability of allocating individual $i$ to the $k$-th component – required for Step (i) above – is given by:

[TABLE]

where $b$ is a normalising constant that ensures that $\sum_{k=1}^{K}p(z_{i}=k|D,{\bf z}_{-i},\ldots)=1$, and $n_{-i,k}$ is the number of individuals currently allocated to component $k$, excluding the $i$-th individual. That is, if we let $\mathbbm{1}_{a}$ denote the indicator function (which is equal to 1 if $a$ is true and zero otherwise), then $n_{-i,k}=\sum_{z_{j}\in{\bf z}_{-i}}\mathbbm{1}_{z_{j}=k}$.

The Dirichlet process (DP) mixture model may be derived by considering the limit $K\rightarrow\infty$, in which case the conditional posterior probability of allocating individual $i$ to an existing component (i.e. one to which other individuals are currently allocated) is given by:

[TABLE]

where $n_{-i,k}$ is as before and $b$ is again a normalising constant (see Equation (8) below). For the DP, we also require the conditional posterior probability of allocating individual $i$ to a new component, which is:

[TABLE]

The normalising constant, $b$, in Equations (6) and (7) is chosen to ensure that:

[TABLE]

Given values for the component allocation variables, ${\bf z}$, it is relatively straightforward to perform Step (ii); i.e. to sample the parameters $\boldsymbol{\Theta},\boldsymbol{\Phi},\boldsymbol{\beta}$. See Appendix for full details, where we also describe how to sample the DP hyperparameter, $\alpha$, according to the method described in Escobar and West (1995).

3 Variable selection in profile regression

The multi-view model that we propose in Section 4 may be regarded as an extended form of variable selection for clustering. We therefore start by presenting a version of the profile regression model that permits variable selection in the same manner as in the models proposed by Law et al. (2003, 2004); Tadesse et al. (2005); Papathomas et al. (2012).

We again consider the component-conditional likelihood associated with the $i$-th individual, making the additional assumption of conditionally independent clustering variables:

[TABLE]

where $P$ is the dimension of ${\bf x}$, so ${\bf x}_{i}=[x_{i1},\ldots,x_{iP}]^{\top}$. In the following, we assume throughout a common parametric form, $f_{x}$, for all clustering variables; however, we could straightforwardly extend to allow different parametric forms for different clustering variables in order to allow, for example, modelling of mixed (continuous and categorical) data.

In order to perform variable selection, we follow the same approach as in Law et al. (2003, 2004) and Tadesse et al. (2005), and introduce latent binary indicator variables $\gamma_{j}$, such that $\gamma_{j}=1$ if the $j$-th clustering variable contributes to the clustering structure, and 0 otherwise. We model the clustering variables for which $\gamma_{j}=0$ as having density $q_{x}(x_{j}|\boldsymbol{\phi}_{0j})$, where $\boldsymbol{\phi}_{0j}$ is a “global” (i.e. not component-specific) parameter and $q_{x}$ need not be of the same parametric form as $f_{x}$. We denote the collection of all $\boldsymbol{\phi}_{0j}$ parameters by $\boldsymbol{\phi}_{0\cdot}=\{\boldsymbol{\phi}_{0j}\}_{j=1}^{P}$. The allocation-conditional likelihood is then:

[TABLE]

The introduction of the latent indicator variables in Equation (10) results in the grouping of clustering variables into two disjoint sets. Adopting the terminology of Cui et al. (2007), we will refer to these sets as views. We define View 1 to be the set of clustering variables (for which $\gamma_{j}=1$) that define a clustering structure, and View 0 to be the set (for which $\gamma_{j}=0$) that do not define a clustering structure. We will refer to this latter set of clustering variables as the null view. The variable selection model may therefore be considered to be one which not only groups individuals together (into clusters), but also groups clustering variables together (into 2 views). Inference of the view allocation variables, $\gamma_{j}$, can be performed via Gibbs sampling, and amounts to sampling according to the posterior probabilities associated with each of the models (i.e. the “no clustering structure” model, $q_{x}$, and the “clustering structure” model, $f_{x}$) given the observed data for the $j$-th clustering variable, $\{x_{ij}\}_{i=1}^{n}$. See Appendix for details.

4 Multi-view Bayesian profile regression

A natural extension to the variable selection model is to allow there to be $L$ views (with $L\in\mathbb{Z}^{+}$ now allowed to be more than 2) by allowing the indicator variables, $\gamma_{j}$, to be categorical with $L$ categories, $0,1,2,\ldots,L-1$. It is then clear that the $\gamma_{j}$ variables are view allocation variables, which serve an analogous role to the component allocation variables, $z_{i}$, but which act to group together the clustering variables rather than the individuals.

As in the variable selection model, the introduction of multiple views necessitates the introduction of view-specific models. In the variable selection case, the clustering variables in View 1 contribute to a mixture model (which defines clusters among the individuals), while the clustering variables in the null view contribute to a single density (corresponding to all individuals being in a single cluster). In general, different choices are possible for the view-specific models. Here we consider an $L$-view case in which View 0 is a null view, while Views $1,\ldots,L-1$ are each associated with a mixture model, such that each of these views defines a different clustering structure among the individuals. We associate the response with View 1 only, and refer to this view as the relevant view. Thus, the View 1 model is a (semi-supervised) Bayesian profile regression model, while the models for Views $2,\ldots,L-1$ are all unsupervised mixture models.

4.1 Allocations-conditional likelihood

For each of the non-null views, we introduce component allocation variables $z_{i\ell}$, such that $z_{i\ell}$ denotes the component that is responsible for the $i$-th individual in the $\ell$-th view’s mixture model. We denote the number of mixture components in the $\ell$-th view’s mixture model by $K_{\ell}$, where – as previously – $K_{\ell}$ may be finite or infinite. We moreover denote the parameters associated with the $k$-th component in the $\ell$-th view by $\boldsymbol{\phi}_{\ell k}$ and define $\boldsymbol{\phi}_{\ell\cdot}=\{\boldsymbol{\phi}_{k\ell}\}_{k=1}^{K_{\ell}}$ to be the full complement of component-specific parameters associated with the $\ell$-th view. The $\boldsymbol{\theta}$ parameters of the model for ${\bf y}$ are associated with view $\ell=1$ only.

The likelihood associated with the $i$-th individual, conditioned on these component allocations, is then:

[TABLE]

Defining $\Gamma_{\ell}=\{j:\gamma_{j}=\ell\}$ to be the index set of clustering variables allocated to the $\ell$-th view, the above may alternatively be written as follows:

[TABLE]

where each of the bracketed terms corresponds to a different view, as shown. This expression makes clear that (conditioned on the allocation of clustering variables to views) each of the views is modelled independently, as also illustrated in Figure 1, with View 0 being modelled as a single cluster, View 1 being modelled with a profile regression model (compare to Equation (9)), and the remaining views each modelled by their own (unsupervised) mixture model.

4.2 Inference

Conditioned on the view allocations, the models describing each of the views are independent, and hence we can perform inference for the component allocations and parameters within each view using existing approaches either for Bayesian profile regression models (Molitor et al., 2010) in the case of View 1, or for unsupervised Bayesian mixture models (e.g. Richardson and Green, 1997; Neal, 2000) for Views 2, …, $L-1$. We therefore adopt a Gibbs sampling approach, in which we iterate between sampling the view allocation variables and performing inference for the models in each view. We provide details of the update for the view allocation indicators below. Updates for the within-view parameters and latent variables are performed as in Neal (2000); see Appendix for further details.

4.2.1 Updating the view allocation indicators

Define $\nu_{\ell}$ to be the prior probability that a clustering variable is allocated to the $\ell$-th view. Given the component allocations $z_{\cdot\ell}=\{z_{i\ell}\}_{i=1}^{n}$ for the $\ell$-th view (for $\ell=1,\ldots,L-1$), the posterior probability that the $j$-th clustering variable is allocated to the $\ell$-th view is then:

[TABLE]

where $Z$ is a normalising constant that ensures that the posterior view allocation probabilities sum to 1, and, as before, $\boldsymbol{\phi}_{\ell z_{i\ell}}$ denotes the parameter associated with the $z_{i\ell}$-th component in the $\ell$-th view, and $\boldsymbol{\phi}_{\ell\cdot}$ is the full complement of component-specific parameters associated with the $\ell$-th view.

For the null view, we have:

[TABLE]

and it is now clear that

[TABLE]

In the case where conjugate priors are taken for the component-specific parameters, $\boldsymbol{\phi}_{\ell k}$, and/or the parameters of the null view model, $\boldsymbol{\phi}_{0.}$, these parameters may be integrated out and the likelihood functions $f_{x}$ and $q_{x}$ may be replaced with marginal likelihoods; see Appendix.

4.3 Choice of $L$, the number of views

It may be noted that inference for the view allocations is somewhat analogous to inference for the mixture component allocations within each (non-null) view. In the latter case, individuals are allocated with a higher probability to components in which individuals with similar clustering variable profiles are allocated; whereas in the former case, clustering variables are allocated with a higher probability to views in which clustering variables defining a similar clustering structure are allocated.

Similarly, the choice of the number of views in a multi-view model, $L$, is analogous to the choice of the number of components, $K$, in a conventional mixture model. Moreover, in much the same way that the prior component allocation probabilities in a mixture model may be treated as parameters, $\pi_{k}$, to be inferred (e.g. taking a Dirichlet or Dirichlet process prior), it is also possible to treat the prior view allocation probabilities in the multi-view model, $\nu_{\ell}$, as parameters. By adopting a Dirichlet (or Dirichlet process) prior, and taking $L$ to be large (or infinite), the number of views may, in principle, be inferred automatically. However, this comes with an associated computational expense, since each additional view brings with it a mixture model whose parameters and latent component allocations must be inferred. Approximate inference procedures for these models, such as variational Bayes (previously considered in the context of multi-view clustering by Guan et al., 2010) will be an important direction for future research, but in the present work we focus upon small values of $L\geq 2$, for which inference via Markov chain Monte Carlo (MCMC) is feasible.

4.4 Initialisation

We have found that a good initialisation strategy is to start with all variables in the relevant view, so that irrelevant and null variables are “selected out” at subsequent iterations. We adopt this initialisation strategy in all examples

5 Examples

Although until now we have deliberately kept the exposition general, in the examples that follow we restrict our attention to Dirichlet process mixture models for the non-null views. In Section 5.1 we present simulation examples that allow us to illustrate how clustering approaches that do not model multiple clustering structures can fail, even if they exploit response information and employ variable selection. In Section 5.2 we consider an integrative clustering example in the context of breast cancer subtype characterisation, in which we fit a multi-view model with 3 views (including one null view).

5.1 Simulation study to illustrate the limitations of methods that ignore multiple views

We construct a simulation example to demonstrate the limitations of existing semi-supervised clustering approaches that fail to account for multiple clustering structures. We consider simulated datasets with $n=300$ individuals, $p=10$ categorical clustering variables (each of which has 3 categories), and a univariate binary response, $y$. The clustering variables define 2 views, with the first $q$ clustering variables ($q\in\{3,5,7,9\}$) defining a relevant clustering structure (which relates to the response) and the remaining $p-q$ defining an irrelevant clustering structure (unrelated to the response). The clustering variables in the relevant view (View 1) define 6 equally-sized clusters, which are related to the response according to $P(y_{i}=1|z_{i1}=k)=\theta_{k}$, where $\theta_{k}\in\{0.01,0.15,0.40,0.60,0.85,0.90\}$. The clustering variables in the irrelevant view (View 2) also define 6 clusters, but there is no link between these clusters and the response; i.e. $P(y_{i}=1|z_{i2}=k)=P(y_{i}=1)=0.485$, regardless of the component allocation in View 2. An illustration of a dataset for $q=5$ is provided in Figure 2, with Figure 2(a) showing the clustering structure defined by the clustering variables in View 1, and Figure 2(b) showing the clustering structure defined by the clustering variables in View 2 (which is irrelevant for the response shown).

Within each view, the categorical data are simulated according to a mixture distribution. With probability $w\in[0,1]$ we simulate $x_{ij}$ according to:

[TABLE]

and with probability $1-w$ we simulate $x_{ij}$ according to a discrete uniform distribution on the three categories. The parameter $w$ thereby allows us to control cluster separability. The clusters are easily separable when $w$ is close to 1, but the clusters are increasingly noisy and hard to separate as $w$ approaches 0.

5.1.1 Results

We applied an existing semi-supervised clustering approach to each dataset, as implemented in the PReMiuM R package for Bayesian profile regression (Liverani et al., 2015). For each dataset, we ran PReMiuM both with and without variable selection. To run with variable selection, we specified the “varSelectType” option in PReMiuM to be “continuous”, which performs variable selection with latent selection weights as described in Papathomas et al. (2012). To obtain a “gold standard” that reflects the performance that could be achieved if only the relevant clustering variables were selected, we additionally applied PReMiuM to each dataset after having removed all of the irrelevant clustering variables. In practical examples, identifying the relevant clustering variables a priori in this way will not be possible; however, here we suppose we have access to an oracle that can provide this information. In all cases, we used PReMiuM to perform 10,000 Gibbs sampling iterations, discarding the first 1,000 as burn in and thinning to retain every 5-th draw.

To summarise the MCMC output, we calculated the adjusted Rand index (ARI) between the true clustering structure in the relevant view and each retained clustering structure (partition) sampled from the posterior. Thus, for each distinct PReMiuM run, we obtained a distribution of ARI values, which assesses how well the inferred clustering structure matches the clustering structure in the relevant view. The ARI can take a value of at most 1 (indicating a perfect match), while a value of 0 indicates that the match is no better than we would expect by random chance.

Figure 3 illustrates typical output for 4 datasets, each having a different value for $q$ (the number of relevant clustering variables). Corresponding figures for other simulated datasets are provided in Supplementary Results, and are qualitatively similar. As we might expect, as $q$ diminishes, so too does our ability to infer the correct clustering structure. Perhaps less intuitively, we also see that – as a consequence of failing to model the multiple views – running PReMiuM with variable selection does not necessarily help to improve inference of the relevant clustering structure, and can actually be damaging. Figure 3 shows that variable selection can be useful when the number of relevant clustering variables is large relative to the number of irrelevant clustering variables (e.g. if 7 out of 10 are relevant). In these cases, the irrelevant clustering variables are discarded (i.e. strongly down-weighted), resulting in performance that is comparable with that achieved when the irrelevant clustering variables are artificially removed (e.g. consider $q=7$ in Figure 3). However, when the number of relevant clustering variables is small relative to the number of irrelevant clustering variables (e.g. if only 3 out of 10 are relevant), variable selection can diminish performance; as, in these cases, it is the relevant clustering variables that are discarded. Since PReMiuM models the data as possessing only one (non-null) clustering structure, it tends to home in on a single “dominant” clustering structure (here, the one that is defined by the majority of clustering variables). Variable selection reinforces this by removing or down-weighting any clustering variables that define a different clustering structure, resulting in performance that is better if the dominant clustering also happens to be the relevant one, and worse if not.

Crucially, we stress that the results seen here are not specific either to the inference or variable selection procedures implemented in PReMiuM, but are a consequence of failing to model multiple (non-null) views in the data, when they exist.

5.2 Breast cancer subtyping

In recent years, many authors (e.g. Shen et al., 2009; Kirk et al., 2012; Lock and Dunson, 2013; Savage et al., 2013) have considered the problem of how to perform integrative clustering, using multiple ’omics datasets in order to better characterise cancer subtypes at the molecular level (see also Kristensen et al., 2014, for a reveiw). Although multi-view approaches have not previously been applied for this purpose, they straightforwardly permit integrative clustering. In the multi-view model, clustering variables are allocated to the same view if they define the same clustering structure. Whether or not these clustering variables come from the same dataset or are of the same data type is irrelevant; all that matters is the clustering structure that they define.

Here we apply semi-supervised multiview modelling to identify clusters among breast cancer tumour samples on the basis of reverse phase protein array (RPPA) and micro-RNA (miRNA) data from TCGA (The Cancer Genome Atlas, 2012). A great deal of existing work has considered the use of mRNA expression data for identifying cancer subtypes, including the PAM50 predictive model for classifying breast cancer tumour samples on the basis of the expression of 50 genes (Parker et al., 2009). We consider $n=108$ breast cancer tumour samples, 66 of which have been classified on the basis of mRNA expression data as basal-like and 42 as Luminal A. We use this classification as a binary response, to guide the clustering on the basis of the RPPA and miRNA data from TCGA. The RPPA data comprise measurements on 171 proteins, while the miRNA data comprise measurements on 423 miRNAs.

Before clustering, we process the miRNA and RPPA datasets as in Lock and Dunson (2013). To robustify against misspecification of models for continuous data, we take the additional pre-processing step of using tertile discretisation within each tumour sample, and treat the resulting data as categorical. After pre-processing, the datasets are concatenated, so that our final working dataset has $p=594$ clustering variables (corresponding to 171 proteins and 423 miRNAs).

5.2.1 Results

We fitted our semi-supervised multi-view model to the concatenated TCGA miRNA and RPPA data, assuming $L=3$ views (1 relevant, 1 irrelevant, and 1 null view). We performed 10,000 Gibbs sampling iterations, removing the first 5,000 as burn in, and then thinning to retain every 5-th draw. For each clustering variable, we calculated Monte Carlo estimates of the probability of being selected into each view. Bar plots of these probabilities are shown in Figure 4. The majority of clustering variables are selected with high probability into the null view, with relatively few selected with high probability into the other two views.

To summarise the clustering structure within each non-null view, we first calculated the posterior similarity matrix (PSM) for each view. A PSM is an $n\times n$ matrix whose $i,j$-entry indicates the proportion of clusterings sampled from the posterior in which tumour sample $i$ and $j$ had the same cluster label (Fritsch and Ickstadt, 2009). Since PSMs summarise pairwise co-clustering probabilities, they provide a summary of the MCMC output that avoids challenges associated with label-switching. The PSMs for the non-null views are provided in Figure 5. It is clear from Figure 5a that, as desired, the clustering structure in the relevant view has a strong association with the subtype label. As shown in Figure 5b, the irrelevant view also possesses a strong clustering structure, but one which is not associated with the subtype label.

To visualise the clusters in the relevant view, we thresholded the selection probabilities to retain only those clustering variables selected with probability at least 0.90 of being selected into the relevant view. This left only 53 clustering variables, of which 32 were miRNAs and 21 were proteins. A heatmap representation of the data for these 53 clustering variables is provided in Figure 6a, from which both the clustering structure and its association with the response (subtype) is clear. We similarly visualised the clusters in the irrelevant view by retaining only those clustering variables selected with probability at least 0.90 of being selected into the irrelevant view. This left 76 clustering variables, all of which were proteins. The resulting heatmap representation is provided in Figure 6b. In this case, while there is an evident clustering structure, it is not associated with tumour subtype.

6 Discussion

We have demonstrated that, when there are multiple clustering structures present in data, existing (single view) clustering approaches can fail to recover the most relevant clustering structure, even when guided by an appropriate response (Section 5.1.1). Moreover, traditional variable selection approaches for clustering do not necessarily improve matters, since they tend to select variables that define the dominant clustering structure, regardless of whether or not it is associated with a response of interest. In Section 5.2.1, we have shown that real molecular datasets can and do possess multiple clustering structures, and that our semi-supervised multi-view model can allow both relevant and irrelevant structures to be identified.

While multi-view approaches clearly provide advantages relative to existing (single view) alternatives, computational considerations are a notable challenge. Dirichlet process mixture models are already computationally costly, and the multi-view approach proposed here introduces an additional Dirichlet process mixture model for each additional view. While restricting the number of views provided adequate results in the examples considered here, this need not be the case if the “true” number of views in the data is greater than the selected value for $L$. Both computational approaches (such as parallelisation, as in Suchard et al., 2010) and fast approximate inference procedures (such as variational Bayes, as in Guan et al., 2010) will be important considerations for future work.

Funding

PDWK and SR were supported by the MRC (MC_UU_00002/13 and MC_UU_00002/10) respectively. FP and SR were supported by EPSRC project EP/R018561/1.

7 Appendix

7.1 Unsupervised Bayesian mixture models

We start by considering unsupervised mixture models of the following general form:

[TABLE]

where $\boldsymbol{\pi}=[\pi_{1},\ldots,\pi_{K}]$ is the vector of mixture weights and $\boldsymbol{\phi}=\{\boldsymbol{\phi}_{1},\ldots,\boldsymbol{\phi}_{K}\}$ is the collection of all component-specific parameters. We shall initially consider the case of finite $K$, and return later to the infinite mixture model.

As is common for mixture models, we introduce latent component allocation variables, $z_{i}$, where $z_{i}=k$ if the $i$-th observation ${\bf x}_{i}$ is associated with the $k$-th component, and $p(z_{i}=k|\pi)=\pi_{k}$. Then,

[TABLE]

and hence

[TABLE]

Integrating out $z_{i}$ by summing over all $K$ possible values, we obtain (as we would hope):

[TABLE]

Making the usual conditional independence assumptions, the full joint model for ${\bf x}_{i},z_{i},\boldsymbol{\phi},\boldsymbol{\pi}$ is:

[TABLE]

where we assume independent priors for the component-specific parameters, $\phi_{k}$, and a symmetric Dirichlet prior for the mixture weights, $\pi_{1},\ldots,\pi_{K}\sim\mbox{Dir}(\alpha/K)$.

For the full dataset, we have:

[TABLE]

7.1.1 Inference via Gibbs sampling (finite $K$ case)

Given Equation (24), it is straightforward to write down the conditionals for Gibbs sampling. For the time being, we assume finite $K$ (from which we will later derive the infinite limit).

Conditional for $\boldsymbol{\phi}_{k}$

By examination of the RHS of Equation 24, we have:

[TABLE]

where $\boldsymbol{\phi}_{-k}$ denotes the set comprising all $\boldsymbol{\phi}_{j}$ for which $j\neq k$. Thus the conditional for $\boldsymbol{\phi}_{k}$ is the posterior density for $\boldsymbol{\phi}_{k}$ given all ${\bf x}_{i}$ for which $z_{i}=k$. If conjugate priors are taken, this posterior is available analytically, otherwise samples may be drawn by, for example, the Metropolis-Hastings algorithm. Note that if there are no ${\bf x}_{i}$ for which $z_{i}=k$ (i.e. if the $k$-th component has no observations associated with it), then $\boldsymbol{\phi}_{k}$ is simply sampled from the prior, $p(\boldsymbol{\phi}_{k})$.

Conditional for $\pi$

By examination of the RHS of Equation 24, we have:

[TABLE]

Hence, the conditional for $\pi$ is the posterior for $\pi$ given the values taken by the categorical latent allocation variables, $z_{i}$, $i=1,\ldots,n$. If we take a conjugate Dirichlet prior, this posterior is available in closed form.

Conditional for $z_{i}$

By examination of the RHS of Equation 24, we have:

[TABLE]

where $z_{-i}$ denotes the set comprising all $z_{j}$ for which $j\neq i$. Since $\sum_{k=1}^{K}p(z_{i}=k|{\bf x}_{1},\ldots,{\bf x}_{n},\boldsymbol{\phi},\pi,z_{-i})=1$, it follows that the conditional is:

[TABLE]

which may be straightforwardly evaluated for finite $K$.

Marginalising $\pi$

Taking a conjugate Dirichlet prior for $\pi$, an alternative strategy is to marginalise $\boldsymbol{\pi}$ rather than to sample it. We assume a symmetirc Dirichlet prior with concentration parameter $\alpha/K$.

Note that the $z_{i}$’s are only conditionally independent of one another given $\boldsymbol{\pi}$, so if we marginalise $\boldsymbol{\pi}$ then we must be careful to model the dependence of $z_{i}$ on $z_{-i}$ in our conditional for $z_{i}$.

We have

[TABLE]

To marginalise $\boldsymbol{\pi}$, we must therefore evaluate $\int_{\pi}p(z_{i}=k|z_{-i},\boldsymbol{\pi},\alpha)p(\boldsymbol{\pi}|\alpha)d\boldsymbol{\pi}=p(z_{i}=k|z_{-i},\alpha)$, which is the conditional prior for $z_{i}$ given the values for the other latent allocation variables, $z_{-i}$.

We have,

[TABLE]

where in the final line we exploit the fact that the $z_{i}$’s are conditionally independent of $\alpha$, given $\pi$.

In order to proceed, we must evaluate this fraction. To do this we require a standard result about Dirichlet distributions, which says that moments of random variables distributed according to a symmetric Dirichlet distribution with parameter $\alpha/K$ can be expressed as follows:

[TABLE]

where the $m_{k}$’s are any natural numbers.

Moreover, we note the following two equalities:

[TABLE]

and

[TABLE]

where $n_{-i,c}$ is the number of $z_{j}$’s with $j\neq i$ for which $z_{j}=c$. It then follows that we may use the result given in Equation (34) in order to evaluate the numerator and denominator in the RHS of Equation (33). After some algebra, and exploiting the property of Gamma functions that $\Gamma(t+1)=t\Gamma(t)$, we obtain:

[TABLE]

Hence,

[TABLE]

Moreover, since $K$ is finite, we may straightforwardly evaluate the equality:

[TABLE]

where

[TABLE]

Marginalising $\phi$

Similarly, if a conjugate prior is available for the $\phi_{k}$’s, then these may be marginalised too. Note that (similar to the case with the $z_{i}$’s when we marginalised $\pi$) the ${\bf x}_{i}$’s are only conditionally independent of one another given the $\phi_{k}$’s and the $z_{i}$’s, so if we marginalise $\phi$ then we must be careful to model the dependence of $v_{i}$ on $v_{-i}$ in our conditional for $z_{i}$.

After some algebra, it is straightforward to show that marginalising $\phi$ gives the following for the conditional for $z_{i}$:

[TABLE]

where ${\bf x}_{-i,k}$ denotes all observations ${\bf x}_{j}$ for which $j\neq i$ and $z_{j}=k$, and hence $p(\phi_{k}|{\bf x}_{-i,k})$ is the posterior for $\phi_{k}$ given all of the observations currently associated with component $k$ (excluding ${\bf x}_{i}$). If there are no ${\bf x}_{j}$ for which $j\neq i$ and $z_{j}=k$ (i.e. if $k$ is a component to which no other observations have been allocated), then we say that the $k$-th component is empty and define $p(\phi_{k}|{\bf x}_{-i,k}):=p(\phi_{k})$ to be the prior for $\phi_{k}$.

When implementing the sampler, it is useful to observe that

[TABLE]

Hence, still assuming that the $k$-th component is not empty, the integral in Equation (39) is

[TABLE]

which is a ratio of marginal likelihoods: one in which we include ${\bf x}_{i}$ amongst the observations associated with component $k$, and one in which we exclude ${\bf x}_{i}$ from the observations associated with component $k$.

This expression aids the interpretation of the sampler: at each iteration, and for each component, we weigh the evidence that ${\bf x}_{i}$ is associated with component $k$ against the evidence that ${\bf x}_{i}$ is not associated with component $k$ (given the other observations currently associated with that component, ${\bf x}_{-i,k}$). Intuitively, this expression ensures that we are more likely to allocate ${\bf x}_{i}$ to a component to which similar observations have previously been allocated.

Note also that the term $\frac{n_{-i,k}+\alpha/K}{n-1+\alpha}$ in Equation (39) represents the conditional prior probability that ${\bf x}_{i}$ should be allocated to component $k$ represents our prior belief that we should allocate ${\bf x}_{i}$ to component $k$, given the allocation of all of the other observations. Since $n_{-i,k}$ is in the numerator, this expresses a “rich-get-richer” prior belief; i.e. that, a priori, we are more likely to assign ${\bf x}_{i}$ to a component that already has many observations assigned to it, rather than to one with fewer.

Final note

Note that, having marginalised the $\pi_{k}$’s and $\phi_{k}$’s, we may use Equation (39) to sample just the $z_{i}$’s, without having to sample any other parameters. The one exception is the $\alpha$ hyperparameter, which we may either fix or sample (using, for example, the approach described in Escobar and West, 1995).

7.1.2 Conditionals for Gibbs sampling ($K$ infinite)

The infinite case can be derived from the finite case by considering the limit $K\rightarrow\infty$. The key change is to the conditional priors for the component allocation variables, $p(z_{i}=k|z_{-i},\alpha)$. There are two cases to consider:

1. $k$ is the label for an existing (non-empty) component:

This is the case where $k=z_{j}$ for some $j\neq i$. By considering the limit as $K\rightarrow\infty$ of Equation (35), we obtain the following:

[TABLE]

2. $k$ is a new label:

This is the case where $k\neq z_{j}$ for all $j\neq i$. Hence, the probability we seek is $p(z_{i}\neq z_{j}\mbox{ for all }j\neq i|z_{-i},\alpha)$. Calculation of this probability is aided by using the following identity:

[TABLE]

where $\mathcal{C}_{-i}=\{k:z_{j}=k\mbox{ for some }j\neq i\}$ is the set of all current component labels (excluding $z_{i}$). Using this identity, we obtain (after a little algebra) the following:

[TABLE]

Conditional for $z_{i}$

Using Equations (41) and (42), it is straightforward to show (cf. Equation (39), and making use of Equation (40)) that we have the following conditionals for $z_{i}$:

[TABLE]

where $b$ is a normalising constant that ensures that

[TABLE]

Up to changes in notation and expansion of some terms, these expressions are identical to those given in Equation (3.7) of Neal (2000).

7.2 Profile regression for categorical data

Recall that our data comprises input-output (i.e. clustering variable vector-response) pairs, $D=\{{\bf v}_{i}\}_{i=1}^{n}$, where ${\bf v}_{i}=[{\bf x}_{i},y_{i}]^{\top}$, with ${\bf x}_{i}$ and $y_{i}$ denoting the $i$-th clustering variable vector and response, respectively. We now consider the case in which the clustering variables and response are all categorical.

7.2.1 Modelling the clustering variables

We now assume that each clustering variable (i.e. each element of the vector ${\bf x}_{i}$) is categorical, with the $j$-th clustering variable having $R_{j}$ categories, which we label as $1,2,\ldots,R_{j}$. We model the data using categorical distributions. We define $\phi_{k,j,r}$ to be the probability that the $j$-th clustering variable takes value $r$ in the $k$-th component, and write $\Phi_{k,j}=[\phi_{k,j,1},\phi_{k,j,2},\ldots,\phi_{k,j,R_{j}}]$ for the collection of probabilities associated with the $j$-th clustering variable in the $k$-th component. We further define $\Phi_{k}=\{\Phi_{k,1},\Phi_{k,2},\ldots,\Phi_{k,J}\}$ to be the collection of all probabilities (over all $J$ clustering variables) associated with the $k$-th component, and $\Phi=\{\Phi_{k}\}_{k\in\mathcal{C}}$ to be the set of all $\Phi_{k}$’s that are associated with non-empty components (here, $\mathcal{C}=\{k:z_{i}=k\mbox{ for some }i\in\{1,\ldots,n\}\}$).

We assume that the clustering variables are conditionally independent, given their component allocation, so that

[TABLE]

Conditional for $\Phi_{k,j}$

From Equation (25), the conditional that we require for Gibbs sampling is the posterior for $\Phi_{k,j}$, given the observations associated with the $k$-th component. For each $j$, we adopt a conjugate Dirichlet prior for $\Phi_{k,j}$,

[TABLE]

where ${\bf a}_{j}=[a_{j,1},\ldots,a_{j,R_{j}}]$ is the vector of concentration parameters. The posterior is then:

[TABLE]

where ${\bf x}_{i_{1}},{\bf x}_{i_{2}},\ldots,{\bf x}_{i_{n_{k}}}$ are the observations associated with component $k$, and $s_{k,j,r}$ is defined to be the number of observations associated with component $k$ for which the $j$-th clustering variable is in category $r$.

Marginalising $\Phi_{k,j}$

We may also integrate out $\Phi_{k,j}$ in order to write down the marginal likelihood associated with ${x}_{i_{1},j},{x}_{i_{2},j},\ldots,{x}_{i_{n_{k}},j}$. Note that the marginal likelihood is (by definition) the prior expectation of the product $\phi_{k,j,1}^{s_{k,j,1}}\ldots\phi_{k,j,{R_{j}}}^{s_{k,j,{R_{j}}}}$. We may therefore use the same standard result that was used to derive Equation (34) in order to immediately write down the marginal likelihood. Still assuming that ${\bf x}_{i_{1}},{\bf x}_{i_{2}},\ldots,{\bf x}_{i_{n_{k}}}$ are the observations associated with component $k$, we have:

[TABLE]

To shorten notation, define $X_{k}=\{{\bf x}_{i_{1}},{\bf x}_{i_{2}},\ldots,{\bf x}_{i_{n_{k}}}\}$ to be the set of observations associated with component $k$, and $X_{k,j}=\{{x}_{i_{1},j},{x}_{i_{2},j},\ldots,{x}_{i_{n_{k}},j}\}$ to be the set containing the $j$-th elements of the vectors in $X_{k}$. Since we assume that the clustering variables are conditionally independent, given their component allocation, it follows that:

[TABLE]

where $p(X_{k,j}|{\bf a}_{j})=p({x}_{i_{1},j},{x}_{i_{2},j},\ldots,{x}_{i_{n_{k}},j}|{\bf a}_{j})$ is as given in Equation (48).

7.2.2 Modelling the response

We assume that the response, $y_{i}$, is also categorical (this time with $R_{y}$ categories), and model using a categorical distribution. Similar to before, we define $\theta_{k,r}$ to be the probability that $y$ takes value $r$ in the $k$-th component, and write $\Theta_{k}=[\theta_{k,1},\theta_{k,1},\ldots,\theta_{k,{R_{y}}}]$ to be the collection of all probabilities associated with the $k$-th component. Also define $Y_{k}=\{y_{i_{1}},y_{i_{2}},\ldots,y_{i_{n_{k}}}\}$ to be the set of $y$’s allocated to component $k$.

Conditional for $\theta_{k}$

We adopt a conjugate Dirichlet prior for $\theta_{k}$,

[TABLE]

where ${\bf a}_{y}=[a_{y,1},\ldots,a_{j,R_{y}}]$ is the vector of concentration parameters. Similar to previously, the posterior for $\theta_{k}$ is then

[TABLE]

where $s_{k,y,r}$ is the number of observations in $Y_{k}$ that are in the $r$-th category.

Marginalising $\theta_{k}$

As before, we may also calculate the marginal likelihood, $p(Y_{k}|{\bf a}_{j})$. We have

[TABLE]

7.2.3 Joint marginal likelihood

In order to proceed, we need an expression for the marginal likelihood associated with $V_{k}=\{X_{k},Y_{k}\}$. We assume that $y$ and ${\bf x}$ are conditionally independent, given their component allocation. It follows that

[TABLE]

where the expressions for $p(Y_{k}|{\bf a}_{y})$ and $p(X_{k,j}|{\bf a}_{j})$ are as given previously. Note that:

Setting $V_{k}=\{{\bf x}_{i},{\bf x}_{-i,k}\}$, we can evaluate the marginal likelihood in the numerator in Equation (43). 2. 2.

Setting $V_{k}=\{{\bf x}_{-i,k}\}$, we can evaluate the marginal likelihood in the denominator in Equation (43). 3. 3.

Setting $V_{k}=\{{\bf x}_{i}\}$, we can evaluate the marginal likelihood in Equation (44).

Thus, we may evaluate all of the terms required for the conditionals for $z_{i}$, and hence (leaving aside, for the time being, the problem of sampling $\alpha$) we have everything we need in order to perform inference for our model.

7.3 Variable selection

7.3.1 A null model for the clustering variable data

We introduce a null model for the ${\bf x}$ data, under the assumption that there is no clustering structure. Under the null model, we again model the clustering variable data using categorical distributions, and – similar to previously – define $\Phi_{0,j}=[\phi_{0,j,1},\phi_{0,j,2},\ldots,\phi_{0,j,R_{j}}]$ and $\Phi_{0}=\{\Phi_{0,1},\Phi_{0,2},\ldots,\Phi_{0,J}\}$. However, under the null model, we assume that there is no clustering structure (or, equivalently, that the data form a single large cluster). The likelihood associated with ${\bf x}_{i}$ therefore does not involve a component allocation variable, $z_{i}$, and is instead:

[TABLE]

In our null model, the observations ${\bf x}_{1},\ldots,{\bf x}_{n}$ are assumed to be conditionally independent given $\boldsymbol{\Phi}_{0}$. It follows that, under the null model, the likelihood associated with the full clustering variable dataset is:

[TABLE]

7.3.2 Variable selection

We consider the possibility that only some of the clustering variables are relevant for determining the clustering structure, while the others are irrelevant. We introduce binary indicators $\gamma_{j}$, such that $\gamma_{j}$ = 1 means that the $j$-th clustering variable is relevant for the clustering structure, while $\gamma_{j}$ = 0 means that the variable is irrelevant.

Informally, we wish to use the full mixture model for the relevant variables, and the null model for the irrelevant variables. More precisely, we consider the following model for ${\bf x}$:

[TABLE]

7.3.3 Inference

We consider how the introduction of the $\gamma_{j}$’s affects how we infer the model parameters. As previously, we start by writing down the full joint model:

[TABLE]

Inferring $\Phi_{k,j}$

As before, we adopt independent Dirichlet priors for $\Phi_{k,j}$,

[TABLE]

where ${\bf a}_{j}=[a_{j,1},\ldots,a_{j,R_{j}}]$ is the vector of concentration parameters.

It is clear from Equation (56) that the conditional for $\Phi_{k,j}$ is

[TABLE]

If $\gamma_{j}=1$, then – as before – the conditional for $\Phi_{k,j}$ is the posterior for $\Phi_{k,j}$ given all observations currently associated with component $k$, i.e.

[TABLE]

where the $s_{k,j,r}$’s are as previously defined

If $\gamma_{j}=0$, then the conditional for $\Phi_{k,j}$ is just the prior, $p(\Phi_{k,j})$, as provided in Equation (57). [But note that, in practice, we would never need to sample from the conditional in this case].

Inferring $\Phi_{0,j}$

As we did for $\Phi_{k,j}$, we adopt a Dirichlet prior for $\Phi_{0,j}$,

[TABLE]

where ${\bf a}_{j}^{(0)}=[a_{j,1}^{(0)},\ldots,a_{j,R_{j}}^{(0)}]$ is the vector of concentration parameters. It is clear from Equation (56) that the conditional for $\Phi_{0,j}$ is

[TABLE]

If $\gamma_{j}=1$, then the conditional for $\Phi_{0,j}$ is just the prior, $p(\Phi_{0,j})$, as provided in Equation (60). [But note that, in practice, we would never need to sample from the conditional in this case].

If $\gamma_{j}=0$, then the conditional for $\Phi_{0,j}$ is the posterior for $\Phi_{0,j}$ given all observations, i.e.

[TABLE]

where $s_{j,r}$ is defined to be the number of observations for which the $j$-th variable is in category $r$.

7.3.4 Conditional for $z_{i}$

Let us temporarily return to the case of finite $K$, and recall Equation (37), which provides the conditional for $z_{i}$ (assuming that the $\theta_{k}$’s are sampled, rather than integrated out). Introducing the $\gamma_{j}$’s (to allow for variable selection), and exploiting the conditional independence of the variables given $z_{i}$, Equation (37) becomes:

[TABLE]

Marginalising $\theta_{k}$

After some algebra, it is straightforward to show that marginalising the component-specific parameters, $\theta_{k}$, gives the following for the conditional for $z_{i}$:

[TABLE]

where

• •

$\mathcal{M}_{-i}=\{m:m\in\{1,\ldots,n\}\mbox{ and }m\neq i\}$; and

• •

$\mathcal{M}_{-i,k}=\{m:m\in\{1,\ldots,n\}\mbox{ and }m\neq i\mbox{ and }z_{m}=k\}$.

Thus, as previously, $p(\theta_{k}|\{{y}_{m}\}_{m\in\mathcal{M}_{-i,k}})$ is the posterior for $\theta_{k}$ given all observations (excluding $y_{i}$) currently associated with the $k$-th component (and similarly for $\Phi_{k,j}$). Also as before, if the $k$-th component is empty we define $p(\theta_{k}|\{{y}_{m}\}_{m\in\mathcal{M}_{-i,k}}):=p(\theta_{k})$ to be the prior for $\theta_{k}$ (and similarly for $\Phi_{k,j}$).

Note that the expression involving $\Phi_{0,j}$ is slightly different: here, we have the posterior for $\Phi_{0,j}$ given all observations except the $i$-th. Since this expression does not involve $k$, we may absorb this term into the normalising constant, $Z$, so that:

[TABLE]

Conditional for $z_{i}$ ($K$ infinite)

Proceeding as previously, we consider the limit as $K\rightarrow\infty$. In this case, the conditionals for $z_{i}$ are as follows:

If $k=z_{j}$ for some $j\neq i$, then

[TABLE]

Moreover,

[TABLE]

where $b$ is a normalising constant that ensures that

[TABLE]

Note that these expressions are essentially the same as those provided in Equations (43) and (44), except that the irrelevant variables do not contribute to the conditionals (as we would hope). The marginal likelihoods may be evaluated as previously described (see Sections 7.2.1 and 7.2.2).

7.3.5 Sampling $\gamma_{j}$

By examination of Equation (56), the conditional for $\gamma_{j}$ is immediately

[TABLE]

Since $\gamma_{j}\in\{0,1\}$, we have:

[TABLE]

and

[TABLE]

where

[TABLE]

Marginalising $\Phi_{k,j}$ and $\Phi_{0,j}$

As previously, we may also marginalise the $\Phi_{k,j}$ and $\Phi_{0,j}$ parameters. Define $X_{k}=\{{\bf x}_{i_{1}},{\bf x}_{i_{2}},\ldots,{\bf x}_{i_{n_{k}}}\}$ to be the set of observations associated with component $k$, and $X_{k,j}=\{{x}_{i_{1},j},{x}_{i_{2},j},\ldots,{x}_{i_{n_{k}},j}\}$ to be the set containing the $j$-th elements of the vectors in $X_{k}$ (as before). Also define $\mathcal{C}=\{k:z_{i}=k\mbox{ for some }i\}$ to be the set of all current component labels, and $X_{\cdot,j}=\cup_{k\in\mathcal{C}}\{X_{k,j}\}$ to be the set containing the $j$-th elements of all observations.

We then obtain

[TABLE]

and

[TABLE]

where $z^{\prime}$ is a normalising constant that ensures that the probabilities sum to 1. The marginal likelihoods may be evaluated as previously described (see Section 7.2.1).

7.4 Multi-view clustering

The method for variable selection described in the previous section may be considered to be a special case of multi-view clustering, in which we allow the data to possess multiple clustering structures (depending upon which group of variables we consider to be “relevant”). In the variable selection case, we model the data as possessing 2 clustering structures: one non-trivial clustering structure (to which the “relevant” variables contribute), and one trivial clustering structure (comprising just a single cluster, to which the “irrelevant” variables contribute). This may be straightforwardly extended to $L$ clustering structures, by allowing $\gamma_{j}$ to be a categorical variable whose value indicates the clustering structure to which the $j$-th variable contributes.

7.4.1 Notational conventions

We assume that $\gamma_{j}\in\{0,1,2,\ldots,V-1\}$ is a categorical variable with $L$ categories. As previously, we assume that if $\gamma_{j}=0$ then the $j$-th variable is completely irrelevant, and does not contribute to any clustering structure (except the trivial structure comprising one cluster).

[TABLE]

[TABLE]

For $\ell=1,\ldots,L-1$, we assume that there are non-overlapping variable sets $\mathcal{S}_{\ell}$, such that each $\mathcal{S}_{\ell}$ possesses its own clustering structure. Within each $\mathcal{S}_{\ell}$, we model the data using a mixture model. To this end, we introduce latent component allocation variables, $z_{i}^{(\ell)}$, which indicate the component in the model for $\mathcal{S}_{\ell}$ that is responsible for generating the $i$-th observation. Similarly, we define $\Phi_{k}^{(\ell)}$ to be the parameters associated with the $k$-th component in the $\ell$-th mixture model.

We moreover define $\gamma_{j}=\ell$ if the $j$-th variable is in $\mathcal{S}_{\ell}$. Note that, a priori, we do not know which variables belong to which $\mathcal{S}_{\ell}$, so we must perform inference for the $\gamma_{j}$’s.

Finally, we introduce another categorical variable, $\nu\in\{1,2,\ldots,L-1\}$, such that the clustering structure present in $\mathcal{S}_{\nu}$ is the one that is relevant for profile regression, i.e. the likelihood for $y_{i}$,

[TABLE]

depends only on the $z_{i}^{(\ell)}$ component allocation variable.

Our conditional model for ${\bf x}_{i}=[{\bf x}_{i},y_{i}]^{\top}$ (given the $z_{i}^{(\ell)}$’s, $\gamma_{j}$’s, and $\nu$) is then

[TABLE]

7.4.2 Conditionals for Gibbs sampling

We may follow the same method of argument presented in Section 3 in order to derive the required conditionals. For brevity, we provide the conditionals obtained after integrating out the component-specific parameters.

Conditionals for $z_{i}^{(\ell)}$: $\nu=\ell$ case

In this case, we are dealing with the clustering structure that is relevant for profile regression, so must include the contribution of the responses variable, $y$. We have (cf. Equations (67) and (68)): If $k=z_{j}^{(\ell)}$ for some $j\neq i$, then

[TABLE]

where $n_{-i,k}^{(\ell)}$ is the number of observations, ${\bf x}_{m}$, for which $z_{m}^{(\ell)}=k$ (excluding the $i$-th observation), and $\mathcal{M}_{-i,k}^{(\ell)}=\{m:m\in\{1,\ldots,n\}\mbox{ and }m\neq i\mbox{ and }z_{m}^{(\ell)}=k\}$.

Moreover,

[TABLE]

where $b^{(\ell)}$ is a constant that ensures that the probabilities sum to 1, and $\alpha^{(\ell)}$ is the concentration parameter of the Dirichlet process prior on the mixture weights for the $\nu$-th mixture model.

Conditionals for $z_{i}^{(\ell)}$: $\nu\neq\ell$ case

In this case, the response variable, $y$, does not contribute to the conditionals. We have (cf. Equations (67) and (68)): If $k=z_{j}^{(\ell)}$ for some $j\neq i$, then

[TABLE]

where $n_{-i,k}^{(\ell)}$ is the number of observations, ${\bf x}_{m}$, for which $z_{m}^{(\ell)}=k$ (excluding the $i$-th observation). Moreover,

[TABLE]

where $b^{(\ell)}$ is a constant that ensures that the probabilities sum to 1, and $\alpha^{(\ell)}$ is the concentration parameter of the Dirichlet process prior on the mixture weights for the $\ell$-the mixture model.

7.4.3 Conditional for $\gamma_{j}$

Define $X_{k}^{(\ell)}=\{{\bf x}_{i_{1}},{\bf x}_{i_{2}},\ldots,{\bf x}_{i_{n_{k}^{(\ell)}}}\}$ to be the set of observations associated with component $k$ of the $\ell$-th mixture model, and $X_{k,j}^{(\ell)}=\{{x}_{i_{1},j},{x}_{i_{2},j},\ldots,{x}_{i_{n_{k}^{(\ell)}},j}\}$ to be the set containing the $j$-th elements of the vectors in $X_{k}^{(\ell)}$. Also define $\mathcal{C}^{(\ell)}=\{k:z_{i}^{(\ell)}=k\mbox{ for some }i\}$ to be the set of all current component labels for the $\ell$-th mixture model.

We have

[TABLE]

and, for all other $\ell\neq 0$,

[TABLE]

where $z^{\prime}$ is a normalising constant that ensures that the probabilities sum to 1. The marginal likelihoods may be evaluated as previously described (see Section 7.2.1).

7.4.4 Conditional for $\nu$

The conditional for $\nu$ is rather similar to the conditional for $\gamma_{j}$. Define $Y_{k}^{(\ell)}=\{{y}_{i_{1}},{y}_{i_{2}},\ldots,{y}_{i_{n_{k}^{(\ell)}}}\}$. Then, for $\ell\in\{1,\ldots,L-1\}$, we have:

[TABLE]

where $z^{\prime\prime}$ is a normalising constant to ensure that probabilities sum to 1.

Note that, when sampling $\nu$, we are effectively comparing different models for partitioning the $y$’s, and selecting from amongst these according to probabilities given by Equation (81) above.