Are Clusterings of Multiple Data Views Independent?

TL;DR

This paper investigates whether clusters derived from different data types in a multi-view dataset are dependent or independent, developing a new statistical test and applying it to wellness data to assess the validity of combined clustering approaches.

Contribution

The paper introduces a novel test for independence of clusterings across multiple data views and applies it to real-world wellness data to evaluate clustering assumptions.

Findings

Clusterings within each data view are dependent across time.

Clusterings across different data types are largely independent.

Implication: combining data views for clustering may not always be appropriate.

Abstract

In the Pioneer 100 (P100) Wellness Project (Price and others, 2017), multiple types of data are collected on a single set of healthy participants at multiple timepoints in order to characterize and optimize wellness. One way to do this is to identify clusters, or subgroups, among the participants, and then to tailor personalized health recommendations to each subgroup. It is tempting to cluster the participants using all of the data types and timepoints, in order to fully exploit the available information. However, clustering the participants based on multiple data views implicitly assumes that a single underlying clustering of the participants is shared across all data views. If this assumption does not hold, then clustering the participants using multiple data views may lead to spurious results. In this paper, we seek to evaluate the assumption that there is some underlying relationship among the clusterings from the different data views, by asking the question: are the clusters within each data view dependent or independent? We develop a new test for answering this question, which we then apply to clinical, proteomic, and metabolomic data, across two distinct timepoints, from the P100 study. We find that while the subgroups of the participants defined with respect to any single data type seem to be dependent across time, the clustering among the participants based on one data type (e.g. proteomic data) appears not to be associated with the clustering based on another data type (e.g. clinical data).

Tables

Table 1. Table 1: Results from the test of H 0 : C = 1 K ( 1 ) ​ 1 K ( 2 ) T : subscript 𝐻 0 𝐶 subscript 1 superscript 𝐾 1 superscript subscript 1 superscript 𝐾 2 𝑇 H_{0}:C=1_{K^{(1)}}1_{K^{(2)}}^{T} developed in Section 3.1 applied to clinical, proteomic, and metabolomic data at the first and third timepoints, and applied to pairs of data views defined by different data types. Sample sizes n 𝑛 n , dimensions in each view p 1 subscript 𝑝 1 p_{1} and p 2 subscript 𝑝 2 p_{2} , and p-values obtained using the permutation approximation from Section 3.3 are reported.

View 1	View 2	$n$	$p_1$	$p_2$	p-value
Clinical at Timepoint 1	Clinical at Timepoint 3	83	204	198
Proteomic at Timepoint 1	Proteomic at Timepoint 3	66	249	257
Metabolomic at Timepoint 1	Metabolomic at Timepoint 3	88	641	640
Clinical at Timepoint 1	Proteomic at Timepoint 1	70	204	249	0.236
Clinical at Timepoint 2	Proteomic at Timepoint 2	60	205	254	0.091
Clinical at Timepoint 3	Proteomic at Timepoint 3	66	198	257	0.950
Clinical at Timepoint 1	Metabolomic at Timepoint 1	98	204	641	0.034
Clinical at Timepoint 2	Metabolomic at Timepoint 2	89	205	641	0.073
Clinical at Timepoint 3	Metabolomic at Timepoint 3	81	198	640	0.328
Proteomic at Timepoint 1	Metabolomic at Timepoint 1	72	249	641	0.402
Proteomic at Timepoint 2	Metabolomic at Timepoint 2	67	254	641	0.004
Proteomic at Timepoint 3	Metabolomic at Timepoint 3	73	257	640	0.020

Equations

$$f(X^{(l)}\mid Z^{(l)}=k)={\varphi }^{(l)}\left(X^{(l)};{\theta }_k^{(l)}\right)$$

$$P(Z^{(1)}=k,Z^{(2)}=k^{\prime })={\Pi }_{k{k}^{\prime }}$$

$$f(X^{(1)},X^{(2)};{\theta }^{(1)},{\theta }^{(2)},\Pi )$$

$$=\sum _{k,k^{\prime }}{\Pi }_{k{k}^{\prime }}f(X^{(1)}\mid Z^{(1)}=k)f(X^{(2)}\mid Z^{(2)}=k^{\prime })$$

$$=\sum _{k,k^{\prime }}{\Pi }_{k{k}^{\prime }}{\varphi }^{(1)}\left(X^{(1)};{\theta }_k^{(1)}\right){\varphi }^{(2)}\left(X^{(2)};{\theta }_{k^{\prime }}^{(2)}\right),$$

$$\left\{\Pi \in {\Delta }^{K^{(1)}\times K^{(2)}}:\Pi 1_{K^{(2)}}={\pi }^{(1)},{\Pi }^T{1}_{K^{(1)}}={\pi }^{(2)}\right\}=\left\{\mathrm{diag}({\pi }^{(1)})C\mathrm{diag}({\pi }^{(2)}):C\in {\mathcal{C}}_{{\pi }^{(1)},{\pi }^{(2)}}\right\},$$

$$f(X^{(1)},X^{(2)};{\theta }^{(1)},{\theta }^{(2)},\Pi )$$

$$\equiv f(X^{(1)},X^{(2)};{\theta }^{(1)},{\theta }^{(2)},{\pi }^{(1)},{\pi }^{(2)},C).$$

$$f(X^{(l)};{\theta }^{(l)},{\pi }^{(l)})$$

$$\ell ({\theta }^{(1)},{\theta }^{(2)},{\pi }^{(1)},{\pi }^{(2)},C)=\sum _{i=1}^n\log f(X_i^{(1)},X_i^{(2)};{\theta }^{(1)},{\theta }^{(2)},{\pi }^{(1)},{\pi }^{(2)},C),$$

$$\ell ({\theta }^{(l)},{\pi }^{(l)})=\sum _{i=1}^n\log f(X_i^{(l)};{\theta }^{(l)},{\pi }^{(l)}),$$

$$\hat{C}$$

$${\hat{\varphi }}_{ik}^{(1)}={\varphi }^{(1)}\left(X_i^{(1)};{\hat{\theta }}_k^{(1)}\right)\text{and}{\hat{\varphi }}_{i{k}^{\prime }}^{(2)}={\varphi }^{(2)}\left(X_i^{(2)};{\hat{\theta }}_{k^{\prime }}^{(2)}\right).$$

$$f(X^{(1)},X^{(2)};{\theta }^{(1)},{\theta }^{(2)},{\pi }^{(1)},{\pi }^{(2)},C)=f(X^{(1)}\mid X^{(2)};{\theta }^{(1)},{\theta }^{(2)},{\pi }^{(1)},{\pi }^{(2)},C)f(X^{(2)};{\theta }^{(2)},{\pi }^{(2)}).$$

$$\log \Lambda$$

$$=\underset{{\theta }^{(1)},{\theta }^{(2)},{\pi }^{(1)},{\pi }^{(2)},C}{\sup }\ell ({\theta }^{(1)},{\theta }^{(2)},{\pi }^{(1)},{\pi }^{(2)},C)-[\ell ({\hat{\theta }}^{(1)},{\hat{\pi }}^{(1)})+\ell ({\hat{\theta }}^{(2)},{\hat{\pi }}^{(2)})],$$

$$\ell ({\theta }^{(1)},{\theta }^{(2)},{\pi }^{(1)},{\pi }^{(2)},C)=\ell ({\theta }^{(1)},{\pi }^{(1)})+\ell ({\theta }^{(2)},{\pi }^{(2)}),$$

$$\log \tilde{\Lambda }$$

$$=\ell ({\hat{\theta }}^{(1)},{\hat{\theta }}^{(2)},{\hat{\pi }}^{(1)},{\hat{\pi }}^{(2)},\hat{C})-[\ell ({\hat{\theta }}^{(1)},{\hat{\pi }}^{(1)})+\ell ({\hat{\theta }}^{(2)},{\hat{\pi }}^{(2)})]$$

$$=\sum _{i=1}^n\log \left[\frac{({\hat{\varphi }}_i^{(1)}{)}^T\mathrm{diag}({\hat{\pi }}^{(1)})\hat{C}\mathrm{diag}({\hat{\pi }}^{(2)}){\hat{\varphi }}_i^{(2)}}{({\hat{\varphi }}_i^{(1)}{)}^T{\hat{\pi }}^{(1)}({\hat{\pi }}^{(2)}{)}^T{\hat{\varphi }}_i^{(2)}}\right].$$

$$\Pi =\frac{1-\delta }{K^2}{1}_K{1}_K^T+\frac{\delta }{K}{I}_K,$$

$$G^2({\hat{M}}^{(1)},{\hat{M}}^{(2)})=2\sum _{k=1}^{K^{(1)}}\sum _{k^{\prime }=1}^{K^{(2)}}{\hat{N}}_{k{k}^{\prime }}\log \left[(n{\hat{N}}_{k{k}^{\prime }})/({\hat{N}}_{k.}{\hat{N}}_{.k^{\prime }})\right],$$

$$G^2({\hat{M}}^{(1)},{\hat{M}}^{(2)})\stackrel{d}{\to }{\chi }_{(K^{(1)}-1)(K^{(2)}-1)}^2.$$

$$I({\hat{M}}^{(1)},{\hat{M}}^{(2)})=G^2({\hat{M}}^{(1)},{\hat{M}}^{(2)})/2n.$$

$$\log \tilde{\Lambda }\left(\hat{\Pi },\{{\hat{r}}_i^{(1)},{\hat{r}}_i^{(2)}{\}}_{i=1}^n\right)=\sum _{i=1}^n\log \left[\frac{({\hat{r}}_i^{(1)}{)}^T\hat{\Pi }{\hat{r}}_i^{(2)}}{({\hat{r}}_i^{(1)}{)}^T(\hat{\Pi }{1}_{K^{(2)}})(1_{K^{(1)}}^T\hat{\Pi }){\hat{r}}_i^{(2)}}\right],$$

$${\hat{M}}_i^{(1)}\equiv \underset{k\in \{1,2,\dots ,K^{(1)}\}}{\text{arg max}}{\hat{r}}_{ik}^{(1)},{\hat{M}}_i^{(2)}\equiv \underset{k^{\prime }\in \{1,2,\dots ,K^{(2)}\}}{\text{arg max}}{\hat{r}}_{i{k}^{\prime }}^{(2)}.$$

$$\log \tilde{\Lambda }(\hat{N}/n,\{e_{{\hat{M}}_i^{(1)}},e_{{\hat{M}}_i^{(2)}}{\}}_{i=1}^n)=nI({\hat{M}}^{(1)},{\hat{M}}^{(2)})=G^2({\hat{M}}^{(1)},{\hat{M}}^{(2)})/2,$$

$$\left\{\Pi \in {\Delta }^{K^{(1)}\times K^{(2)}}:\Pi 1_{K^{(2)}}={\pi }^{(1)},{\Pi }^T{1}_{K^{(1)}}={\pi }^{(2)}\right\}\subseteq \left\{\mathrm{diag}({\pi }^{(1)})C\mathrm{diag}({\pi }^{(2)}):C\in {\mathcal{C}}_{{\pi }^{(1)},{\pi }^{(2)}}\right\}.$$

$$C{\pi }^{(2)}=\left[\begin{array}{c}\sum _{k^{\prime }}{C}_{1k^{\prime }}{\pi }_{k^{\prime }}^{(2)}\\ ⋮\\ \sum _{k^{\prime }}{C}_{K^{(1)}{k}^{\prime }}{\pi }_{k^{\prime }}^{(2)}\end{array}\right]=\left[\begin{array}{c}\sum _{k^{\prime }}\frac{{\Pi }_{1k^{\prime }}}{{\pi }_1^{(1)}}\\ ⋮\\ \sum _{k^{\prime }}\frac{{\Pi }_{K^{(1)}{k}^{\prime }}}{{\pi }_{K^{(1)}}^{(1)}}\end{array}\right]=1_{K^{(1)}}$$

$$\left\{\Pi \in {\Delta }^{K^{(1)}\times K^{(2)}}:\Pi 1_{K^{(2)}}={\pi }^{(1)},{\Pi }^T{1}_{K^{(1)}}={\pi }^{(2)}\right\}\supseteq \left\{\mathrm{diag}({\pi }^{(1)})C\mathrm{diag}({\pi }^{(2)}):C\in {\mathcal{C}}_{{\pi }^{(1)},{\pi }^{(2)}}\right\}.$$

$$\Pi =\mathrm{diag}({\pi }^{(1)})C\mathrm{diag}({\pi }^{(2)}).$$

$$\Pi 1_{K^{(2)}}$$

$$=\mathrm{diag}({\pi }^{(1)})C{\pi }^{(2)}$$

$$=\mathrm{diag}({\pi }^{(1)})1_{K^{(1)}}={\pi }^{(1)}.$$

$${\hat{r}}_i^{(1)}=\frac{{\hat{\varphi }}_i^{(1)}}{1^T{\hat{\varphi }}_i^{(1)}},{\hat{r}}_i^{(2)}=\frac{{\hat{\varphi }}_i^{(2)}}{1^T{\hat{\varphi }}_i^{(2)}},$$

Full text

1 Introduction

Complex biological systems consist of diverse components with dynamics that may vary over time, and so these sytems often cannot be fully characterized by any single type of data, or at any single snapshot in time. Consequently, it has become increasingly common for researchers to collect multiple data sets, or views, for a single set of observations. In the machine learning literature, this is known as the multiple-view or multi-view data setting.

Multiple-view data has been applied extensively to characterize disease, such as in The Cancer Genome Atlas Project (Cancer Genome Atlas Research Network, 2008). In contrast, The Pioneer 100 (P100) Wellness Project (Price and others, 2017) collected multiple-view data from healthy participants to characterize wellness, and to optimize wellness of the participants through personalized healthcare recommendations. One way to do this is to identify subgroups of similar participants using cluster analysis, and then tailor recommendations to each subgroup.

In recent years, many papers have proposed clustering methods in the multiple-view data setting (Bickel and Scheffer, 2004; Shen and others, 2009; Kumar and others, 2011; Kirk and others, 2012; Lock and Dunson, 2013; Gabasova and others, 2017). The vast majority of these methods “borrow strength” across the data views to obtain a more accurate clustering of the observations than would be possible based on a single data view. Implicitly, these methods assume that there is a single consensus clustering shared by all data views.

The P100 data contains many data views; multiple data types (e.g. clinical data and proteomic data) are available at multiple timepoints. Thus, it is tempting to apply consensus clustering methods to identify subgroups of the P100 participants. However, before doing so, it is important to check the assumption that there exists a single consensus clustering. If instead different views reflect unrelated aspects of the participants, then there is no “strength to be borrowed” across the views, and it would be better to perform a separate clustering of the observations in each view. Before attempting cluster analysis of the P100 data, it is critical that we determine which combinations of views have “strength to be borrowed”, and which combinations do not.

This raises the natural question of how associated the underlying clusterings are in each view. Suppose we cluster the P100 participants twice, once using their baseline clinical data, and once using their baseline proteomic data. Can we tell from the data whether the two views’ underlying clusterings are related or unrelated? Answering this question provides useful information:

• Case 1: If the underlying clusterings appear related, then this increases confidence that the clusterings are scientifically meaningful, and offers some support for performing a consensus clustering of the P100 participants that integrates baseline clinical and proteomic views.
• Case 2: If the underlying clusterings appear unrelated, we must consider two explanations.

Perhaps clinical and proteomic views measure different properties about the participants, and therefore identify complementary (or “orthogonal”) clusterings. If so, then a consensus clustering is unlikely to provide meaningful results, and may cause us to lose valuable information about the subgroups underlying the individual data views. 2. 2.

Perhaps the subgroups underlying the data views are indeed related, but they appear unrelated due to noise. If so, then we might be skeptical of any results obtained on these very noisy data, whether from consensus clustering or another approach.

In Case 2, it would not be appropriate to perform consensus clustering.

To determine from the data whether the two views’ clusterings are related or unrelated, it is tempting to apply a clustering procedure (e.g. k-means) to each view, then apply well-studied tests of independence of categorical variables (e.g. the $\chi^{2}$-test for independence, the $G$-test for independence, or Fisher’s exact test) to the estimated cluster assignments. However, such an approach relies on an assumption that the estimated cluster assignments are independent and identically distributed samples from the joint distribution of the cluster membership variables, which is not satisfied in practice. Thus, there is a need for an approach which takes into account the fact that the clusterings are estimated from the data.

The rest of this paper is organized as follows. In Section 2, we propose a mixture model for two-view data. In Section 3, we use this model to develop a test of the null hypothesis that clusterings on two views of a single set of observations are independent. We explore the performance of our proposed hypothesis test via numerical simulation in Section 4. In Section 5, we connect and compare our proposed hypothesis test to the aforementioned approach of applying the $G$-test for independence to the estimated cluster assignments, and draw connections between this approach and the mutual information statistic (Meilă, 2007). In Section 6, we apply our method to the clinical, proteomic, and metabolomic datasets from the P100 study. In Section 7, we provide a discussion, which includes the extension to more than two views.

2 A Mixture Model for Multiple-View Data

2.1 Model Specification

In what follows, we consider the case of two data views. We will discuss the extension to more than two views in Section 7.

Suppose we have $p_{1}$ and $p_{2}$ features in the first and second data view, respectively. For a single observation, let $X^{(1)}\in\mathbb{R}^{p_{1}}$ and $X^{(2)}\in\mathbb{R}^{p_{2}}$ denote the random vectors corresponding to the two data views and let $Z^{(1)}\in\{1,\ldots,K^{(1)}\}$ and $Z^{(2)}\in\{1,\ldots,K^{(2)}\}$ be unobserved random variables, indicating the latent group memberships of this observation in the two data views. Here, $K^{(1)}$ and $K^{(2)}$ represent the number of clusters in the two data views, which we assume for now to be known (we will consider the case in which they are unknown in Section 2.4). We assume that $X^{(1)}$ and $X^{(2)}$ are conditionally independent given the pair of cluster memberships, $(Z^{(1)},Z^{(2)})$; this assumption is common in the multi-view clustering literature (see e.g. Bickel and Scheffer 2004, Rogers and others 2008, Kumar and others 2011, Lock and Dunson 2013, Gabasova and others 2017). Further, suppose that

[TABLE]

where $\phi^{(l)}(\cdot;\theta)$ denotes a density function with parameter $\theta$, and $\Pi\in\Delta^{K^{(1)}\times K^{(2)}}\equiv\{S\in\mathbb{R}^{K^{(1)}\times K^{(2)}}:S_{kk^{\prime}}\geq 0,\sum\limits_{k=1}^{K^{(1)}}\sum\limits_{k^{\prime}=1}^{K^{(2)}}S_{kk^{\prime}}=1\}$. Equations (2.1)–(2.2) are an extension of the finite mixture model (McLachlan and Peel, 2000) to the case of two data views. We further assume that each cluster has positive probability, i.e. $P(Z^{(1)}=k)>0$ and $P(Z^{(2)}=k^{\prime})>0$, and so $\Pi 1_{K^{(2)}}\in\Delta_{+}^{K^{(1)}}$ and $\Pi^{T}1_{K^{(1)}}\in\Delta_{+}^{K^{(2)}}$, where $\Delta^{K}_{+}\equiv\{s\in\mathbb{R}^{K}:s_{k}>0,\sum\limits_{k=1}^{K}s_{k}=1\}.$

Let $\theta^{(1)}\equiv(\theta^{(1)}_{1},\ldots,\theta^{(1)}_{K^{(1)}})$ and $\theta^{(2)}\equiv(\theta^{(2)}_{1},\ldots,\theta^{(2)}_{K^{(2)}})$. The joint density of $X^{(1)}$ and $X^{(2)}$ is

[TABLE]

where the second equality follows from conditional independence of $X^{(1)}$ and $X^{(2)}$ given $Z^{(1)}$ and $Z^{(2)}$, and the last equality follows from (2.1).

The matrix $\Pi$ governs the statistical dependence between the two data views. It will be useful for us to parameterize $\Pi$ in terms of a triplet $(\pi^{(1)},\pi^{(2)},C)$ that separates the single-view information from the cross-view information.

Proposition 1**.**

Suppose $\pi^{(1)}\in\Delta^{K^{(1)}}_{+}$ and $\pi^{(2)}\in\Delta^{K^{(2)}}_{+}$. Then,

[TABLE]

where $\mathcal{C}_{\pi^{(1)},\pi^{(2)}}=\{C\in\mathbb{R}^{K^{(1)}\times K^{(2)}}:C_{kk^{\prime}}\geq 0,~{}C\pi^{(2)}=1_{K^{(1)}},~{}C^{T}\pi^{(1)}=1_{K^{(2)}}\}.$

A proof of Proposition 1 is given in Appendix A.1 of the Supplementary Materials.

Proposition 1 indicates that any matrix $\Pi\in\Delta^{K^{(1)}\times K^{(2)}}$ with $\Pi 1_{K^{(2)}}\in\Delta_{+}^{K^{(1)}}$ and $\Pi^{T}1_{K^{(1)}}\in\Delta_{+}^{K^{(2)}}$ can be written as the product of its row sums $\pi^{(1)}$, its column sums $\pi^{(2)}$, and a matrix $C$. Therefore, we can rewrite the joint probability density (2.3) as follows:

[TABLE]

In what follows, we will parametrize the density of $X^{(1)}$ and $X^{(2)}$ in terms of $\theta^{(1)},\theta^{(2)},\pi^{(1)},\pi^{(2)}$, and $C$, rather than in terms of $\theta^{(1)},\theta^{(2)}$, and $\Pi$.

The following proposition characterizes the marginal distributions of $X^{(1)}$ and $X^{(2)}$.

Proposition 2**.**

Suppose $X^{(1)}$ and $X^{(2)}$ have joint distribution (2.1). Then for $l=1,2$, $X^{(l)}$ has marginal density given by

[TABLE]

Proposition 2 follows from (2.1) – (2.2). Proposition 2 shows that for $l=1,2$, $X^{(l)}$ marginally follows a mixture model with parameters $\theta^{(l)}$ and cluster membership probabilities $\pi^{(l)}$. Note that the marginal density of $X^{(1)}$ does not depend on $\theta^{(2)},\pi^{(2)}$, and $C$, and similarly, the marginal density of $X^{(2)}$ does not depend on $\theta^{(1)},\pi^{(1)}$, and $C$; this fact will be critical to our approach to parameter estimation in Section 2.3.

The model described in this section is closely related to several multiple-view mixture models proposed in the literature: see e.g. Rogers and others (2008), Kirk and others (2012), Lock and Dunson (2013), and Gabasova and others (2017). However, the focus of those papers is cluster estimation: they do not provide a statistical test of association, and for the most part, impose additional structure on the probability matrix $\Pi$ in order to encourage similarity between the clusters estimated in each data view. By constrast, the focus of this paper is inference: testing for dependence between the clusterings in different data views. The model described in this section is a step towards that goal.

2.2 Interpreting $\Pi$

In Figures 1(i)–(iii), $n=15$ independent pairs $\{(X^{(1)}_{i},X^{(2)}_{i})\}_{i=1}^{n}$ are drawn from the model (2.1)–(2.2), for three choices of $\Pi$. The left-hand panel represents the $p_{1}=2$ features in the first data view, and the right-hand panel represents the $p_{2}=2$ features in the second data view. For $l=1,2$, the observations $\{X^{(l)}_{i}\}_{i=1}^{n}$ in the $l$th data view belong to two clusters, where the latent variables $\{Z^{(l)}_{i}\}_{i=1}^{n}$ characterize cluster membership in the $l$th data view. Light and dark gray represent the clusters in the first view, and circles and triangles represent the clusters in the second view.

Figures 1(i)–(ii) correspond to two special cases of $\Pi$ that are easily interpretable. In Figure 1(i), $\Pi$ has rank one, i.e. $\Pi=\pi^{(1)}[\pi^{(2)}]^{T}$, so that the clusterings in the two data views are independent. Thus, whether an observation is light or dark appears to be roughly independent of whether it is a circle or a triangle. In Figure 1(ii), $K^{(1)}=K^{(2)}$ and $\Pi$ is diagonal (up to a permutation of the rows), so that the clusterings in the two data views are identical. Thus, all of the circles are light and all of the triangles are dark. Another special case is when $\Pi$ is block diagonal (up to a permutation) with $K_{B}$ blocks. Then, the clusterings of the two data views agree about the presence of $K_{B}$ “meta-clusters” in the data. For example, one clustering might be a refinement of the other, or if one view has clusters $A,B,C,D$, and the other has clusters $a,b,c,d$, it could be that $A\cup B=a$ and $C=b\cup c$ and $D=d$.

In general, $\Pi$ will be neither exactly rank one nor exactly (block) diagonal; Figure 1(iii) provides such an example. Furthermore, $\hat{\Pi}$ (an estimator for $\Pi$) almost certainly will be neither. Nonetheless, examination of $\hat{\Pi}$ can provide insight into the relationships between the two clusterings. For example, if $\hat{\Pi}$ is far from rank one, then this suggests that the clusterings in the two data views may be dependent. We will formalize this intuition in Section 3.

2.3 Estimation

2.3.1 Estimation Procedure and Algorithm

Given $n$ independent pairs $(X^{(1)}_{1},X^{(2)}_{1}),\ldots,(X^{(1)}_{n},X^{(2)}_{n})$ drawn from the model (2.1)–(2.2), the log-likelihood takes the form

[TABLE]

where $f(\cdot,\cdot;\theta^{(1)},\theta^{(2)},\pi^{(1)},\pi^{(2)},C)$ is defined in (2.1). A custom expectation-maximization (EM; Dempster and others 1977; McLachlan and Krishnan 2007) algorithm could be developed to solve (2.6) for a local optimum (a global optimum is typically unattainable, as (2.6) is non-concave). We instead take a simpler approach. Proposition 2 implies that for $l=1,2$, we can estimate $\theta^{(l)}$ and $\pi^{(l)}$ by maximizing the marginal likelihood for the $l$th data view, given by

[TABLE]

where $f(\cdot;\theta^{(l)},\pi^{(l)})$ is defined in (2.5). Each of these maximizations can be performed using standard EM-based software for model-based clustering of a single data view. Let $\hat{\theta}^{(1)},\hat{\pi}^{(1)},\hat{\theta}^{(2)},$ and $\hat{\pi}^{(2)}$ denote the maximizers of (2.7). Next, to estimate $C$, we maximize the joint log-likelihood (2.6) evaluated at $\hat{\theta}^{(1)},\hat{\pi}^{(1)},\hat{\theta}^{(2)},$ and $\hat{\pi}^{(2)}$, subject to the constraints imposed by Proposition 1:

[TABLE]

where $\mathcal{C}_{\hat{\pi}^{(1)},\hat{\pi}^{(2)}}=\{C\in\mathbb{R}^{K^{(1)}\times K^{(2)}}:C_{kk^{\prime}}\geq 0,~{}C\hat{\pi}^{(2)}=1_{K^{(1)}},~{}C^{T}\hat{\pi}^{(1)}=1_{K^{(2)}}\}$. Equation 2.8 is a convex optimization problem, which we solve using a combination of exponentiated gradient descent (Kivinen and Warmuth, 1997) and the Sinkhorn-Knopp algorithm (Franklin and Lorenz, 1989), as detailed in Appendix B of the Supplementary Materials. Details of our approach for fitting the model (2.1)–(2.2) are given in Algorithm 1.

2.3.2 Justification of Estimation Procedure

The estimation procedure in Section 2.3.1 does not maximize the joint likelihood (2.6); nonetheless, we will argue that it is an attractive approach.

To begin, in Step 1 of Algorithm 1, we estimate $\theta^{(1)}$ and $\pi^{(1)}$ by maximizing the marginal likelihood (2.7). This decision leads to computational advantages, as it enables us to make use of efficient software for clustering a single data view, such as the mclust package (Scrucca and others, 2016) in R. We can further justify this decision using conditional inference theory. Equation 3.6 in Reid (1995) extends the definition of ancillary statistics to a setting with nuisance parameters. We show that $X^{(2)}$ is ancillary (in the extended sense of Reid 1995) for $\theta^{(1)},\pi^{(1)}$, and $C$ by using the definition of conditional densities, and Proposition 2, to rewrite (2.1) as

[TABLE]

Thus, Reid (1995) argues that we should use only $X^{(2)}$, and not $X^{(1)}$, to estimate $\theta^{(2)}$ and $\pi^{(2)}$. In Step 1 of Algorithm 1, we are doing exactly this.

In Steps 3–5 of Algorithm 1, we maximize $\ell(\hat{\theta}^{(1)},\hat{\theta}^{(2)},\hat{\pi}^{(1)},\hat{\pi}^{(2)},\cdot)$, giving $\hat{C}$, which is a pseudo maximum likelihood estimator for $C$ in the sense of Gong and Samaniego (1981). This decision also leads to computational advantages, as it enables us to make use of efficient convex optimization algorithms in estimating $C$. Results in Gong and Samaniego (1981) suggest that when $\hat{\theta}^{(1)}$, $\hat{\theta}^{(2)}$, $\hat{\pi}^{(1)}$, and $\hat{\pi}^{(2)}$ are good estimates, $\hat{C}$ is so as well.

2.4 Selection of the Number of Clusters

In Sections 2 and 3, our discussion assumed that $K^{(1)}$ and $K^{(2)}$ are known. However, this is rarely the case in practice. Recall that we estimate $\theta^{(1)}$ and $\pi^{(1)}$ by maximizing the marginal likelihood (2.7), which amounts to performing model-based clustering of $X^{(1)}$ only. Thus, to select the number of clusters $K^{(1)}$, we can make use of an extensive literature (reviewed in e.g. Mirkin 2011) on choosing the number of clusters when clustering a single data view. For example, we can use AIC or BIC to select $K^{(1)}$ and $K^{(2)}$.

3 Testing Whether Two Clusterings are Independent

3.1 A Brief Review of Pseudo Likelihood Ratio Tests

Let $\ell(\alpha,\beta,\gamma)$ be the log-likelihood function for a random sample, where $\mathcal{A}$ is the parameter space of $\alpha$. Given a null hypothesis $H_{0}:\alpha=\alpha_{0}$ for some $\alpha_{0}\in\mathcal{A}$, an alternative hypothesis $H_{1}:\alpha\neq\alpha_{0}$, and an estimator $\hat{\gamma}$, the pseudo likelihood ratio statistic (Self and Liang, 1987) is defined to be $\log\tilde{\Lambda}\equiv\underset{\alpha,\beta}{\sup}~{}\ell(\alpha,\beta,\hat{\gamma})-\underset{\beta}{\sup}~{}\ell(\alpha_{0},\beta,\hat{\gamma})$. Let $\alpha^{*}$ be the true parameter value for $\alpha$. If $\alpha^{*}$ is an interior point of $\mathcal{A}$, then under some regularity conditions, if $H_{0}$ holds, then $2\log\tilde{\Lambda}\overset{d}{\longrightarrow}\chi^{2}_{r},$ where $r$ is the dimension of $\mathcal{A}$ (Chen and Liang, 2010).

3.2 A Pseudo Likelihood Ratio Test for Independence

In this subsection, we develop a test for the null hypothesis that $H_{0}:C=1_{K^{(1)}}1_{K^{(2)}}^{T}$, or equivalently, that $H_{0}:\Pi=\pi^{(1)}(\pi^{(2)})^{T}$: that is, we test whether $Z^{(1)}$ and $Z^{(2)}$ are independent, i.e. whether the cluster memberships in the two data views are independent. We could use a likelihood ratio test statistic to test $H_{0}$,

[TABLE]

where the second equality follows from noticing that substituting $C=1_{K^{(1)}}1_{K^{(2)}}^{T}$ into (2.6) yields

[TABLE]

where $\ell(\theta^{(l)},\pi^{(l)})$ for $l=1,2$ are defined in (2.7), and recalling the definition of $\hat{\theta}^{(1)}$, $\hat{\pi}^{(1)}$, $\hat{\theta}^{(2)}$, and $\hat{\pi}^{(2)}$ as the maximizers of (2.7). However, (3.2) requires maximizing $\ell(\theta^{(1)},\theta^{(2)},\pi^{(1)},\pi^{(2)},C)$, which would require a custom EM algorithm; furthermore, the resulting test statistic will typically involve the difference between two local maxima (since each term in (3.2) requires fitting an EM algorithm). This leads to erratic behavior, such as negative values of $\log\Lambda$.

Therefore, instead of taking the approach in (3.2), we develop a pseudo likelihood ratio test, as in Section 3.1. We use the marginal MLEs, $\hat{\theta}^{(1)},\hat{\pi}^{(1)}$, and $\hat{\theta}^{(2)},\hat{\pi}^{(2)}$, instead of performing the joint optimization in (3.2). This leads to the test statistic

[TABLE]

where $\hat{C}$ in (3.12) is defined in (2.8), $\mathcal{C}_{\hat{\pi}^{(1)},\hat{\pi}^{(2)}}$ is defined in Proposition 1, $\hat{\phi}^{(1)}$ and $\hat{\phi}^{(2)}$ are defined in (2.9), and the last equality follows from (2.6), (2.7), and (2.9). In addition to taking advantage of the computationally efficient estimation procedure described in Section 2.3.1, the pseudo likelihood ratio test statistic does not exhibit the erratic behavior exhibited by the likelihood ratio test statistic. This stability comes from all three terms in (3.12) involving the same local maxima (as opposed to different local maxima).

3.3 Approximating the Null Distribution of $\log\tilde{\Lambda}$

The discussion in Section 3.1 suggests that under $H_{0}:C=1_{K^{(1)}}1_{K^{(2)}}^{T}$, one might expect that $2\log\tilde{\Lambda}\overset{d}{\longrightarrow}\chi^{2}_{r},$ where $r=(K^{(1)}-1)(K^{(2)}-1)$ is the dimension of $\mathcal{C}_{\pi^{(1)},\pi^{(2)}}$. However, this approximation performs poorly in practice, due to violations of the regularity conditions in Chen and Liang (2010). Furthermore, we will often be interested in data applications in which $n$ is relatively small. Hence, we propose a permutation approach. We observe from (3.11) that under $H_{0}$, the log-likelihood is identical under any permutation of the order of the samples in each view. Hence, we take $B$ random permutations of the samples $X^{(2)}_{i}$ from the second view, and compare the observed value of $\log\tilde{\Lambda}$ to its empirical distribution in these permutation samples. Details are given in Algorithm 2. Since $\hat{\phi}^{(1)}$, $\hat{\phi}^{(2)}$, $\hat{\pi}^{(1)}$, and $\hat{\pi}^{(2)}$ are invariant to permutation, for each permutation we need only to estimate $C$. This is another advantage of our test over the likelihood ratio test discussed in Section 3.2, which would require repeating the EM algorithm in every permutation.

Even when we reject the null hypothesis, the clusters could be only weakly dependent; thus, it is helpful to measure the strength of association between the views. Recalling from Section 2.2 that $\text{rank}(\Pi)=1$ implies independence of the clusterings in the two data views, we propose to calculate the effective rank (Vershynin, 2012) of $\hat{\Pi}$, defined in Algorithm 1 – the ratio of the sum of the singular values of $\hat{\Pi}$, and the largest singular value of $\hat{\Pi}$. The effective rank of a matrix is bounded between 1 and its rank, and the matrix is far from rank-1 when its effective rank is far from 1. For example, in Figure 1(iii), the effective rank of $\Pi$ is 1.5, and is upper bounded by 2. Thus, the effective rank of $\hat{\Pi}$ is bounded between 1 and $\min\{K^{(1)},K^{(2)}\}$, and $\hat{\Pi}$ is far from rank-1 when its effective rank is far from 1.

4 Simulation Results

To investigate the Type I error and power of our test, we generate data from (2.1)–(2.2), with

[TABLE]

for $K=6$ and for a range of values of $\delta\in[0,1]$, where $\delta=0$ corresponds to independent clusterings, and $\delta=1$ corresponds to identical clusterings. We draw the observations in the $l$th data view from a Gaussian mixture model, for which the $k$th mixture component is a $N_{p}(\mu^{(l)}_{k},\Sigma^{(l)})$ distribution, with $p=10$, and with $\mu^{(l)}_{k}$ given in Appendix C.1 of the Supplementary Materials.

We simulate 2000 data sets for $\Sigma^{(1)}=\Sigma^{(2)}=\sigma^{2}I_{p}$ for a range of values of $\sigma$ and $n$, and evaluate the power of the pseudo likelihood ratio test of $H_{0}:C=1_{K_{1}}1_{K_{2}}^{T}$ described in Section 3.2 at nominal significance level $\alpha=0.05$, when the number of clusters is correctly and incorrectly specified. To perform Step 1 of Algorithm 1, we use the package mclust in R to fit Gaussian mixture models with a common $\sigma^{2}I_{p}$ covariance matrix (the “EII” covariance structure in mclust). We use $B=200$ permutation samples in Step 2 of Algorithm 2. Simulations in this paper were conducted using the simulator package (Bien, 2016) in R. Results are shown in Figure 2.

The pseudo likelihood ratio test controls the Type I error close to the nominal $\alpha=0.05$ level, even when the number of clusters is misspecified. Power tends to increase as $\delta$ (defined in (4.14)) increases, and tends to decrease as $\sigma$ increases. Compared to using the correct number of clusters, using too many clusters yields lower power, but using too few clusters can sometimes yield higher power (e.g. in the middle panel of Figure 2). This is because, when the signal-to-noise ratio is low, the true clusters are not accurately estimated; thus, combining several true clusters into a single “meta-cluster” can sometimes, but not always, lead to improved agreement between clusterings across the two data views. We explore the impact of the choice of $K$ on the performance of the pseudo likelihood ratio test in Appendix C.2.1 of the Supplementary Materials.

Additional values of $K$ and $p$ are investigated in Appendix C.2.2.

5 Connection to the G-test for Independence and Mutual Information

Let $\hat{M}^{(1)}=(\hat{M}^{(1)}_{1},\ldots,\hat{M}^{(1)}_{n})$ and $\hat{M}^{(2)}=(\hat{M}^{(2)}_{1},\ldots,\hat{M}^{(2)}_{n})$ denote the results of applying a clustering procedure to $X^{(1)}$ and $X^{(2)}$ respectively. In this notation, $\hat{M}^{(1)}_{i}\in\{1,\ldots,K^{(1)}\}$ and $\hat{M}^{(2)}_{i}\in\{1,\ldots,K^{(2)}\}$ denote the estimated cluster assignment for the $i$th observation in the two views. To test whether $Z^{(1)}$ and $Z^{(2)}$ are independent, we could naively apply tests on $\hat{M}^{(1)}$ and $\hat{M}^{(2)}$ for whether two categorical variables are independent. For instance, we could use the $G$-test statistic for independence (Chapter 3.2, Agresti 2003), given by

[TABLE]

where $\hat{N}_{kk^{\prime}}=|\{i\in\{1,\ldots,n\}:\hat{M}^{(1)}_{i}=k,\hat{M}^{(2)}_{i}=k^{\prime}\}|$, $\hat{N}_{.k^{\prime}}=\sum\limits_{k}\hat{N}_{kk^{\prime}}$, and $\hat{N}_{k.}=\sum\limits_{k^{\prime}}\hat{N}_{kk^{\prime}}$. Under the model $\hat{N}_{kk^{\prime}}~{}\overset{\mathrm{ind}}{\sim}\text{ Poisson }(\mu_{kk^{\prime}}),\log\mu_{kk^{\prime}}=\alpha_{k}+\beta_{k^{\prime}}+\gamma_{kk^{\prime}},$ the $G$-test statistic for independence (5.15) is a likelihood ratio test statistic for testing the null hypothesis of independence, i.e. for testing $H_{0}$: $\gamma_{kk^{\prime}}=0\text{ for all }k,k^{\prime}$. Thus, under $H_{0}$: $\gamma_{kk^{\prime}}=0$,

[TABLE]

The $G$-test statistic for independence (5.15) relies on an assumption which is violated in our setting, namely that $\{(\hat{M}^{(1)}_{i},\hat{M}^{(2)}_{i})\}_{i=1}^{n}$ are independent and identically distributed samples from the distribution of $(Z^{(1)},Z^{(2)})$. It is nonetheless a natural approach to the problem of comparing two views’ clusterings. In fact, the mutual information of Meilă (2007) for measuring the similarity between two clusterings of a single $n\times p$ dataset can be written as a scaled version of the $G$-test statistic; when applied to instead measure the similarity between $\hat{M}^{(1)}$ and $\hat{M}^{(2)}$, the mutual information $I(\hat{M}^{(1)},\hat{M}^{(2)})$ is given by

[TABLE]

While the proposed pseudo likelihood ratio test statistic (3.13) for testing independence of $Z^{(1)}$ and $Z^{(2)}$ does not resemble the simple $G$-test statistic for independence in (5.15), we show here that they are in fact quite related.

Let $\hat{r}_{i}^{(1)}=\frac{\hat{\phi}^{(1)}_{i}}{1_{K^{(1)}}^{T}\hat{\phi}^{(1)}_{i}}$ and $\hat{r}_{i}^{(2)}=\frac{\hat{\phi}_{i}^{(2)}}{1_{K^{(2)}}^{T}\hat{\phi}^{(2)}_{i}}$ be the vectors giving the soft-clustering assignment weights (or “responsibilities”) for the $i$th observation in the two views, where $\hat{\phi}_{i}$ is defined in (2.9). We rewrite the pseudo likelihood ratio test statistic (3.13) as

[TABLE]

where $\hat{\Pi}$ is defined in Algorithm 1. In the following proposition, we consider replacing the “soft” cluster assignments $\hat{r}^{(1)}_{i}$ and $\hat{r}^{(2)}_{i}$ with “hard” cluster assignments, and replacing the estimate $\hat{\Pi}$ derived from the “soft” cluster assignments with an estimate derived from “hard” cluster assignments, in (5.18). In what follows,

[TABLE]

Proposition 3**.**

Let $\hat{M}^{(1)}$ and $\hat{M}^{(2)}$ be the estimated model-based cluster assignments in each data view defined by (5.19). Let $\hat{N}$ be the matrix with entries $\hat{N}_{kk^{\prime}}$ containing the number of observations assigned to cluster $k$ in view 1 and cluster $k^{\prime}$ in view 2. Then,

[TABLE]

where $\log\tilde{\Lambda}(\cdot,\cdot)$ is defined in (5.18), and $e_{t}$ is the unit vector that contains a 1 in the $t$th element.

Proposition 3 follows by algebra, and says that replacing the soft cluster assignments in the pseudo likelihood ratio test statistic of Section 3 with hard cluster assignments yields exactly the $G$-test statistic for independence (5.15) (and the mutual information given in (5.17))! In fact, in the special case of fitting multiple-view Gaussian mixtures with common covariance matrix $\sigma^{2}I_{p_{1}}$ in the first view and $\sigma^{2}I_{p_{2}}$ in the second view, we will show that as $\sigma\rightarrow 0$, and the soft cluster assignments converge to hard cluster assignments, the pseudo likelihood ratio test statistic converges to the $G$-test for independence. In what follows, $\log\tilde{\Lambda}\equiv\log\tilde{\Lambda}\left(\hat{\Pi},\{\hat{r}^{(1)}_{i},\hat{r}^{(2)}_{i}\}_{i=1}^{n}\right)$, as in (3.13) and (5.18).

Proposition 4**.**

Let $\sigma^{2}>0$. Suppose that to compute $\log\tilde{\Lambda}$, we fit the model (2.1)–(2.2), for $\phi^{(1)}$ and $\phi^{(2)}$ densities of Gaussian distributions with covariance matrices $\sigma^{2}I_{p_{1}}$ and $\sigma^{2}I_{p_{2}}$ respectively. Let $\tilde{M}^{(1)}$ and $\tilde{M}^{(2)}$ denote the results of applying k-means clustering on the two data views. Then, as $\sigma^{2}\rightarrow 0$, $\log\tilde{\Lambda}\rightarrow nI(\tilde{M}^{(1)},\tilde{M}^{(2)})=G^{2}(\tilde{M}^{(1)},\tilde{M}^{(2)})/2$.

Proposition 4 is proven in Appendix A.2 of the Supplementary Materials. When $\sigma^{2}>0$, the pseudo likelihood ratio test statistic, the $G$-test statistic, and the mutual information are not equivalent. We can thus think of the pseudo likelihood ratio test statistic as reflecting the uncertainty associated with the clusterings obtained on the two views, and the $G$-test statistic and the mutual information as ignoring the uncertainty associated with the clusterings. This suggests that the pseudo likelihood ratio test of Section 3.2 outperforms the $G$-test for independence when the sample size is small and/or there is little separation between the clusters.

To confirm this intuition, we return to the simulation set-up described in Section 4, and compare the performances of the pseudo likelihood ratio test (3.13) and the G-test for independence (5.15) for testing $H_{0}:C=1_{K^{(1)}}1_{K^{(2)}}^{T}$. We obtain p-values for (5.15) using the $\chi^{2}$ approximation from (5.16), and using a permutation approach, where we take $B$ permutations of the elements of $\hat{M}^{(2)}$, and compare the observed value of (5.15) to its empirical distribution in these permutation samples. The results are shown in Figure 3; we see that the two tests yield similar power when the sample size is larger and/or the value of $\sigma$ is smaller, and that the pseudo likelihood ratio test yields higher power than the G-test for independence when the sample size is smaller and/or the value of $\sigma$ is larger. We note that the $\chi^{2}$ approximation for the G-test from (5.16) does not control the Type I error. Additional values of $p$ and $K$, additional values of $\Sigma^{(l)}$, and non-Gaussian finite mixture models are investigated in Appendices C.3.1, C.3.2, and C.3.3 of the Supplementary Materials, respectively; the results are similar to those described in this section.

6 Application to the Pioneer 100 Wellness Project (Price and others, 2017)

6.1 Introduction to the Scientific Problem

In the P100 Wellness Project (Price and others, 2017), multiple biological data types were collected at multiple timepoints for 108 healthy participants. For each participant, whole genome sequences were measured, activity tracking data were collected daily over nine months, and clinical laboratory tests, metabolomes, proteomes, and microbiomes were measured at three-month, six-month, and nine-month timepoints. The P100 study aims to optimize wellness of the participants through personalized healthcare recommendations. In particular, clinical biomarkers measured at baseline were used to make personalized health recommendations.

As an alternative approach, we could identify subgroups of individuals with similar clinical profiles using cluster analysis, and then develop interventions tailored to each subgroup. It is tempting to identify these subgroups using not just clinical data at baseline, but also other types of data (e.g. proteomic data) at other timepoints. We could do this by applying a multi-view consensus clustering method (e.g. Shen and others 2009). However, such an approach assumes that there is a single true clustering underlying all data types at all timepoints. Therefore, before applying a consensus clustering approach, we should determine whether there is any evidence that the clusterings underlying the data types and/or timepoints are at all related (in which case consensus clustering may lead to improved estimation of the clusters) or whether the clusterings are completely unrelated (in which case one would be better off simply performing a separate clustering of the observations in each view). In what follows, we will use the hypothesis test developed in Section 3 to determine whether clusterings of P100 participants based on clinical, proteomic, and genomic data are dependent across timepoints, and across data types.

6.2 Data Analysis

At each of the three timepoints, 207 clinical measurements, 268 proteomic measurements, and 642 metabolomic measurements were available for $n=108$ observations. In the following, we define a data view to be a single data type at a single timepoint. In each view, we removed features missing in more than 25% of participants, and removed participants missing more than 25% of features. Next, features in each view with standard deviation 0 were removed. The remaining missing data were imputed using nearest neighbors imputation in the impute package in R (Hastie and others, 2017). Features in each view were then adjusted for gender using linear regression. Finally, the remaining features were scaled to have standard deviation 1. As in Section 4, we consider the model (2.1)–(2.2) under the assumption that each component in the mixture is drawn from a Gaussian distribution. For each data view, we fit the model using the mclust package in R, with a common $\sigma^{2}I$ covariance matrix (the “EII” covariance structure in mclust). To test $H_{0}:C=1_{K^{(1)}}1_{K^{(2)}}^{T}$, we compute p-values using the permutation approximation discussed in Section 3.3 with $B=10^{5}$. Based on the results in Appendix C.2.1 of the Supplementary Materials, we choose the number of clusters in each view by BIC under the constraint that the number of clusters is greater than one.

We now compare the clusterings in the clinical data at the first and third timepoints, the clustering in the proteomic data at the first and third timepoints, and the clusterings in the metabolomic data at the first and third timepoints. The sample sizes and results are reported in Table 1. For each data type, the clusters found at each timepoint are displayed in Figure 4.

We find strong evidence that for each data type, the clusterings at the first and third timepoints are not independent. We further measure the strength of dependence through the effective rank of $\hat{\Pi}$, as described in Section 3.3. For the clusterings in the clinical data, the effective rank of $\hat{\Pi}$ is 1.63, and is upper bounded by 2. For the clusterings in the proteomic data, the effective rank of $\hat{\Pi}$ is 1.90, and is upper bounded by 5. For the clusterings in the metabolomic data, the effective rank of $\hat{\Pi}$ is 1.2, and is upper bounded by 3. These results suggest that the strengths of association for the clusterings estimated on the clinical data, the proteomic data, and the metabolomic data, are strong, moderate, and weak respectively. The fact that the clusterings estimated on some data types are strongly dependent over time provides evidence that they are scientifically meaningful. Furthermore, it suggests that performing consensus clustering on some data types (e.g. clinical data and proteomic data) across timepoints may be reasonable.

We now focus on comparing clusterings in the clinical, proteomic, and metabolomic data at a single timepoint. The sample sizes and results are reported in Table 1.

The results provide modest evidence that proteomic and metabolomic data at a given timepoint are dependent, and provide weak evidence that clinical and metabolomic data are dependent. However, on balance, the evidence that the clusterings are dependent across data types is weaker than we might expect. This suggests to us that the underlying subgroups defined by the three data types are in fact quite different, and that we should be very wary of performing a consensus clustering type approach across data types, or any analysis strategy that assumes that all three data types are getting at the same set of underlying clusters.

7 Discussion

Most existing work on multiple-view clustering has focused on the problem of estimation: namely, on exploiting the availability of multiple data views in order to cluster the observations more accurately. In this paper, we have instead focused on the relatively unexplored problem of inference: we have proposed a hypothesis test to determine whether clusterings based on multiple data views are independent or associated.

In Section 6, we applied our test to the P100 Wellness Study (Price and others, 2017). We found strong evidence that clusterings based on clinical data and proteomic data persist over time, i.e. that the subgroups defined by the clinical data and the proteomic data are similar at different timepoints. This suggests that if we wish to identify participant subgroups based on (say) clinical data, then it may be worthwhile to apply a consensus clustering approach to the clinical data from multiple timepoints. However, we found only modest evidence that clusterings based on different data types are dependent! This suggests that we should be cautious about identifying participant subgroups by applying consensus clustering across multiple data types, as the clusterings underlying the distinct data types may be quite different.

Throughout this paper, we compared clusterings on $L=2$ data views. We may also wish to compare clusterings across $L>2$ views. Let $X^{(l)}\in\mathbb{R}^{p_{l}}$ for $1\leq l\leq L$ be the random vectors corresponding to the $L$ views. Suppose $X^{(l)}$ are generated according to (2.1) for $1\leq l\leq L$, where $(Z^{(1)},\ldots,Z^{(L)})$ are unobserved multinomial random variables with probabilities given by $P(Z^{(1)}=k_{1},\ldots,Z^{(L)}=k_{L})=\Pi_{k_{1}\ldots k_{L}},$ for $1\leq k_{l}\leq K^{(l)}$ and $1\leq l\leq L$, where the sum of $\Pi$ over all indices is 1 and $\Pi_{k_{1}\ldots k_{L}}\geq 0$. Results analogous to Propositions 1 and 2 hold in this setting. Thus, we can estimate the parameters in the extended model much as we did in Section 2.3.1, replacing the Sinkhorn-Knopp algorithm for matrix balancing with a tensor balancing algorithm (see e.g. Sugiyama and others 2017). To test the null hypothesis that $Z^{(1)},\ldots,Z^{(L)}$ are mutually independent, we can develop a pseudo likelihood ratio test much as we did in Section 3, where instead of permuting the observations in $X^{(2)}$ in Step 2(a) of Algorithm 2, we permute the observations in $X^{(2)},\ldots,X^{(L)}$. Alternatively, one can simply test for pairwise independence between clusterings, instead of testing for mutual independence between clusterings on all views, as we did in Section 6.

An R package titled multiviewtest is available online at https://github.com/lucylgao/multiviewtest and is forthcoming on CRAN. Code to reproduce the data analysis in Section 6, and to reproduce the simulations in Sections 4 and 5 and in Appendix C, are available online at https://github.com/lucylgao/independent-clusterings-code.

Acknowledgments

We thank Nathan Price and John Earls for responding to inquiries about the P100 data, and Will Fithian for a useful conversation. Lucy L. Gao received funding from the Natural Sciences and Engineering Research Council of Canada. Daniela Witten and Jacob Bien were supported by NIH Grant R01GM123993. Jacob Bien was supported by NSF CAREER Award DMS-1653017. Daniela Witten was supported by NIH Grant DP5OD009145, NSF CAREER Award DMS-1252624, and Simons Investigator Award No. 560585. Conflict of Interest: None declared.

Appendix A Proofs

A.1 Proof of Proposition 1

Suppose $\pi^{(1)}\in\Delta^{K^{(1)}}_{+}$ and $\pi^{(2)}\in\Delta^{K^{(2)}}_{+}$, where $\Delta^{K}_{+}\equiv\{s\in\mathbb{R}^{K}:s_{k}>0,\sum\limits_{k=1}^{K}s_{k}=1\}$. Let $\mathcal{C}_{\pi^{(1)},\pi^{(2)}}=\{C\in\mathbb{R}^{K^{(1)}\times K^{(2)}}:C_{kk^{\prime}}\geq 0,~{}C\pi^{(2)}=1_{K^{(1)}},~{}C^{T}\pi^{(1)}=1_{K^{(2)}}\}$. First, we show that

[TABLE]

Suppose $\Pi\in\Delta^{K^{(1)}\times K^{(2)}}=\{S\in\mathbb{R}^{K^{(1)}\times K^{(2)}}:S_{kk^{\prime}}\geq 0,\sum\limits_{k=1}^{K^{(1)}}\sum\limits_{k^{\prime}=1}^{K^{(2)}}S_{kk^{\prime}}=1\}$ such that $\Pi\cdot 1_{K^{(2)}}=\pi^{(1)}$ and $\Pi^{T}\cdot 1_{K^{(1)}}$. Define $C\in\mathbb{R}^{K^{(1)}\times K^{(2)}}_{+}$ by $C_{kk^{\prime}}=\frac{\Pi_{kk^{\prime}}}{\pi^{(1)}_{k}\pi^{(2)}_{k^{\prime}}}.$ The denominator is nonzero because $\pi^{(1)}\in\Delta^{K^{(1)}}_{+}$ and $\pi^{(2)}\in\Delta^{K^{(2)}}_{+}$. Further, $C_{kk^{\prime}}\geq 0$ and $\Pi=\mathrm{diag}(\pi^{(1)})C\mathrm{diag}(\pi^{(2)})$. Now,

[TABLE]

since $\Pi 1_{K^{(2)}}=\pi^{(1)}$ implies that $\sum\limits_{k^{\prime}}\Pi_{kk^{\prime}}=\pi^{(1)}_{k}$. Similarly, we can show $C^{T}\pi^{(1)}=1_{K^{(2)}}$.

Next, we show that

[TABLE]

Suppose $C\in\mathbb{R}^{K^{(1)}\times K^{(2)}}$ such that $C_{kk^{\prime}}\geq 0$, $C\pi^{(2)}=1_{K^{(1)}}$ and $C^{T}\pi^{(1)}=1_{K^{(2)}}$. Define

[TABLE]

Then $\Pi_{kk^{\prime}}=\pi^{(1)}_{k}C_{kk^{\prime}}\pi^{(2)}_{k^{\prime}}$. Since $\pi^{(1)}_{k}>0$, $C_{kk^{\prime}}\geq 0$ and $\pi^{(2)}_{k^{\prime}}>0$, it follows that $\Pi_{kk^{\prime}}\geq 0$. Now,

[TABLE]

Similarly, we can show $\Pi^{T}1_{K^{(1)}}=\pi^{(2)}$. Further, since $\pi^{(1)}\in\Delta_{+}^{K^{(1)}}$, $\sum\limits_{kk^{\prime}}\Pi_{kk^{\prime}}=1_{K^{(1)}}^{T}\Pi 1_{K^{(2)}}=1_{K^{(1)}}^{T}\pi^{(1)}=1,$ so $\Pi\in\Delta^{K^{(1)}\times K^{(2)}}$.

Hence, we have proved Proposition 1.

A.2 Proof of Proposition 4

Let $\sigma^{2}>0$. Suppose throughout that we fit the model (2.1)–(2.2), for $\phi^{(1)}(\cdot;\theta)$ the density of a $N_{p_{1}}(\theta,\sigma^{2}I_{p_{1}})$ random variable and $\phi^{(2)}(\cdot;\theta)$ the density of a $N_{p_{2}}(\theta,\sigma^{2}I_{p_{2}})$ random variable. This amounts to applying Gaussian mixture model-based clustering with common covariance matrix $\sigma^{2}I$ to each view.

Let $\hat{\pi}^{(1)}$, $\hat{\pi}^{(2)}$, $\hat{\theta}^{(1)}$, and $\hat{\theta}^{(2)}$ denote the maximizers of (2.7), and let

[TABLE]

where $\hat{\phi}^{(1)}$ and $\hat{\phi}^{(2)}$ are defined in (2.9).

Since Gaussian mixture model-based clustering with common covariance matrix $\nu^{2}I$ converges to k-means clustering in each view as $\nu^{2}\rightarrow 0$ (see Section 14.3.7 in Friedman and others (2001) for details), as $\sigma^{2}\rightarrow 0$,

[TABLE]

where $\tilde{M}_{i}^{(1)}$ and $\tilde{M}_{i}^{(2)}$ are the estimated k-means cluster assignments of the $i$th observation in each view, $\tilde{N}_{k.}=\sum\limits_{k^{\prime}}\tilde{N}_{kk^{\prime}}$, and $\tilde{N}_{.k^{\prime}}=\sum\limits_{k}\tilde{N}_{kk^{\prime}}$, for $\tilde{N}_{kk^{\prime}}=|\{i\in\{1,\ldots,n\}:\tilde{M}_{i}^{(1)}=k,\tilde{M}_{i}^{(2)}=k^{\prime}\}|.$

We now rewrite (2.8) as

[TABLE]

where the second equality holds because the quantities $\hat{\phi}^{(1)}_{i}$ and $\hat{\phi}^{(2)}_{i}$ defined in (2.9) do not depend on $C$, and the third equality follows from the definition of $\hat{r}^{(1)}_{i}$ and $\hat{r}^{(2)}_{i}$ in (A.1). It follows that

[TABLE]

where

[TABLE]

By (A.2), as $\sigma^{2}\rightarrow 0$, $g^{\sigma}$ converges pointwise to $g$, where

[TABLE]

Applying the method of Lagrange multipliers, we find that $\tilde{C}\equiv\underset{C\in\mathcal{C}_{\hat{\pi}^{(1)},\hat{\pi}^{(2)}}}{\text{arg min}}~{}g(C)$ satisfies

[TABLE]

By Exercise 7.23(c) in Rockafellar and Wets (1998), $\{g^{\sigma}(C)\}_{\{\sigma>0\}}$ is essentially bounded. By Theorem 7.17 in Rockafellar and Wets (1998), the epigraphical limit of $g^{\sigma}$ is $g$. Finally, $g^{\sigma}$ and $g$ are continuous and proper. Hence, by Theorem 7.33 in Rockafellar and Wets (1998),

[TABLE]

By (A.3), (A.4) and (A.5), as $\sigma^{2}\rightarrow 0$ we have

[TABLE]

By (A.2) and (A.6), and the definition of $\hat{\Pi}$ in Algorithm 1,

[TABLE]

Applying (A.2) and (A.7) to (5.18) yields the result.

Appendix B Exponentiated gradient descent for solving (2.8)

After a transformation of the optimization problem, (2.8) can be efficiently solved using exponentiated gradient descent (Kivinen and Warmuth, 1997), a first-order method specially designed for optimization over the simplex. This is a form of mirror descent, with provable convergence results (Beck and Teboulle, 2003). While there is no analytic solution for the update performed at each iteration of the exponentiated gradient descent, each update can be performed by applying the Sinkhorn-Knopp algorithm, a matrix balancing algorithm with provable convergence results and linear convergence rates (Franklin and Lorenz, 1989).

Define

[TABLE]

where $f(\cdot,\cdot;\theta^{(1)},\theta^{(2)},\Pi)$ is defined in (2.3). Consider the following optimization problem:

[TABLE]

By Proposition 1, we can equivalently write (B.1) as follows:

[TABLE]

which is equivalent to

[TABLE]

which is in turn equivalent to

[TABLE]

which is (2.8), the optimization problem we must solve to estimate $\hat{C}$. Hence, to find $\hat{C}$, we can solve (B.1); let $\hat{\Pi}$ be the minimizer of (B.1). Then, $\hat{C}$ can be found by

[TABLE]

The motivation for this transformation of (2.8) is that the transformed problem (B.1) can be efficiently solved using an algorithm described in Cuturi (2013).

We will describe the exponentiated gradient algorithm to solve the general problem

[TABLE]

To solve (B.3), we apply the update

[TABLE]

This is similar to the proximal gradient method, but instead of $\|\Pi-\hat{\Pi}^{t}\|_{F}^{2}/(2s)$ we use the Bregman divergence,

[TABLE]

An advantage of this choice is that the positivity constraint is automatically enforced. For more on this, see Beck and Teboulle (2003). The optimality conditions for the problem (B.4) are

[TABLE]

for $1\leq k\leq K^{(1)}$, $1\leq k^{\prime}\leq K^{(2)}$, where $\lambda_{k}$ and $\eta_{k^{\prime}}$ are Lagrange multipliers for the row sum and column sum constraints, respectively. This implies that

[TABLE]

The gradient is given by

[TABLE]

In the special case of (B.1), writing $G_{kk^{\prime}}=[\nabla g(\hat{\Pi}^{t})]_{kk^{\prime}},$ the update

[TABLE]

can be written as

[TABLE]

where $u_{k}=\exp\{-s\lambda_{k}\}$, $v_{k^{\prime}}=\exp\{-s\eta_{k^{\prime}}\}$, and $M_{kk^{\prime}}=\hat{\Pi}^{t}_{kk^{\prime}}\exp\{sG_{kk^{\prime}}-1\}$.

Since $\hat{\Pi}^{t+1}$ must satisfy the row and column sum constraints, $u$ and $v$ must be chosen accordingly. As in Cuturi (2013), we can apply the Sinkhorn Theorem and Sinkhorn-Knopp algorithm to find $u$ and $v$. By the Sinkhorn Theorem, $u$ and $v$ are unique modulo scalar multiplication of $u$ with a positive number and scalar division of $v$ by that same positive number. The Sinkhorn-Knopp algorithm (alternatively rescaling the columns so that the rows sum to $\hat{\pi}^{(1)}$ then rescaling the rows so that the columns sum to $\hat{\pi}^{(2)}$) can be applied to $M_{kk^{\prime}}$ to find $u$ and $v$. Hence, to perform the update, we simply multiply $\hat{\Pi}^{t}_{kk^{\prime}}$ by $\exp\{sG_{kk^{\prime}}-1\}$ and then apply the Sinkhorn-Knopp algorithm to the updated matrix so that the row and column sum constraints are satisfied. Algorithm S1 provides the details. Using Proposition (1), we can then use our bijection between $\Pi$ to $C$ (B.2) to obtain Algorithm 1 from Algorithm S1.

Appendix C Simulations

C.1 Mean matrices for simulations in Section 4

In the simulations described in Section 4, data are generated from (2.1)–(2.2) with

[TABLE]

for $K=6$. In the $l$th data view, the observations are drawn from a multivariate Gaussian mixture model, for which the $k$th component in the mixture (corresponding to the $k$th cluster) is a $N_{p}(\mu^{(l)}_{k},\sigma^{2}I_{p})$ distribution, with $p=10$. The $p\times K$ mean matrices for the two data views are of the form

[TABLE]

C.2 Supplementary simulations to Section 4

C.2.1 Selection of the number of clusters

Recall from Section 4 that using too few clusters in the pseudo likelihood ratio test sometimes yields better power than using the correct number of clusters. In this section, we demonstrate that using too few clusters can either increase or decrease the power, depending on the situation.

We generate data from (2.1)–(2.2) with $\Pi=\frac{1-\delta}{K^{2}}1_{K}1_{K}^{T}+\frac{\delta}{K}I_{K}$ for $K=6$. In the $l$th data view, the observations are drawn from a bivariate Gaussian mixture model, for which the $k$th component in the mixture (corresponding to the $k$th cluster) is a $N_{2}(\mu^{(l)}_{k},0.4^{2}I_{p})$ distribution. We simulate 2000 datasets for a range of values of $n$, and for two choices of $\mu^{(l)}$:

• Choice 1:

[TABLE]

• Choice 2:

[TABLE]

We evaluate the power of the pseudo likelihood ratio test of $H_{0}:C=1_{K_{1}}1_{K_{2}}^{T}$ at nominal significance level $\alpha=0.05$, when the number of clusters is incorrectly and correctly specified. Results are displayed in Figure S1.

Observe from the left panel of Figure S1 that when $\mu^{(l)}$ is given in (C.2), using too few clusters $(K=2)$ yields higher power than the correct number of clusters $(K=4$). This is because under (C.2), the clusterings on each view contain two natural “meta-clusters”, formed by combining clusters whose means are close; see Figure S2(i) for an illustration. Because the clusters on each view are not well-separated, it is easier to instead cluster the data into the two “meta-clusters”, which are highly in agreement when the four clusters are in agreement. For example, when $\Pi=I_{4}/4$, the $\Pi$ matrix corresponding to the “meta-clustering” is given by $I_{2}/2$. Thus, testing for independence assuming just two clusters yields better power than correctly assuming four clusters.

By contrast, observe from the right panel of Figure S1 that when $\mu^{(l)}$ is given in (C.3), using too few clusters $(K=2)$ yields much lower power than the correct number of clusters ($K=4)$. Again, under (C.3), the clusterings on each view contain two natural “meta-clusters” (see Figure S2(ii)), and it is easier to cluster the data into the two “meta-clusters”. However, under (C.3), even when the four clusters are highly in agreement, the meta-clusters are not highly in agreement. For example, when $\Pi=I_{4}/4$, the $\Pi$ matrix corresponding to the “meta-clustering” is given by $1_{2}1_{2}^{T}/4$. Thus, testing for independence assuming just two clusters yields worse power than correctly assuming four clusters.

We also observe in both panels of Figure S1 that using too many clusters yields slightly lower power than using the correct number of clusters under both (C.2) and (C.3); this result is similar to the results from the simulation setting described in Section 4 and Appendix C.2.2. Furthermore, in both panels of Figure S1, corresponding to the two choices of $\mu^{(l)}$, choosing the clusters using BIC on each view yields slightly lower power than using the correct number of clusters in small and moderate sample sizes, and performs as well as using the correct number of clusters when the sample size is large.

C.2.2 Additional values of $K$ and $p$

We simulate two data views with $p_{1}=p_{2}=p=10$, and $K^{(1)}=K^{(2)}=K=3$. In the $l$th data view, the observations are drawn from a multivariate Gaussian mixture model, for which the $k$th component in the mixture (corresponding to the $k$th cluster) is a $N_{p}(\mu^{(l)}_{k},\sigma^{2}I_{p})$ distribution.

The means of the components in the mixture model, written as a $p\times K$ matrix, are

[TABLE]

Additionally, we simulate two data views with $p_{1}=p_{2}=p=100$, and $K^{(1)}=K^{(2)}=K$, for $K=3$ and $K=6$. In the $l$th data view, the observations are drawn from a multivariate Gaussian mixture model, for which the $k$th component in the mixture (corresponding to the $k$th cluster) is a $N_{p}(\mu^{(l)}_{k},\sigma^{2}I_{p})$ distribution. For $K=3$, the means of the components in the mixture model, written as a $p\times K$ matrix, are

[TABLE]

For $K=6$, the means are

[TABLE]

To investigate the type I error and power of our test, we generate data according to (2.1)–(2.2), with a range of $\Pi$ defined in (C.1).

We simulate 2000 datasets for a range of values of $n$ and $\sigma$, and evaluate the power of the pseudo likelihood ratio test of $H_{0}:C=1_{K_{1}}1_{K_{2}}^{T}$ at nominal significance level $\alpha=0.05$, when the number of clusters is correctly and incorrectly specified. Results are shown in Figures S3, S4, and S5.

Results are similar to the simulation study described in Section 4.

C.3 Supplementary simulations to Section 5

In this section, we consider the $G$-test for independence in Section 5 using $\hat{M}^{(1)}$ and $\hat{M}^{(2)}$ defined in (5.19). As in Section 5, we perform a simulation study in order to compare the performances of the pseudo likelihood ratio test for testing the null hypothesis $H_{0}:C=1_{K^{(1)}}1_{K^{(2)}}^{T}$ and the $G$-test for independence. We obtain p-values for $G^{2}(\hat{M}^{(1)},\hat{M}^{(2)})$ (5.15) using the $\chi^{2}$ approximation from (5.16), as well as a permutation approach, where we take $B$ permutations of the elements of $\hat{M}^{(2)}$, and compare the observed value of $G^{2}(\hat{M}^{(1)},\hat{M}^{(2)})$ to its empirical distribution in these permutation samples.

C.3.1 Additional values of $p$ and $K$

We return to the simulation set-up described in Section 4 and Appendix C.1, with $\Sigma^{(1)}=\Sigma^{(2)}=\sigma^{2}I_{p}$, and investigate a range of values of $p$ and $K$. In addition to the $G$-test for independence and the pseudo likelihood ratio test, we also compute the adjusted Rand Index (ARI) of Hubert and Arabie (1985) in order to compare the results of model-based clustering (implemented as in Sections 2.3 and 3) on each view; p-values for the ARI are obtained using a permutation approach. We compare the performance of the pseudo likelihood ratio test, the $G$-test for independence, and the ARI for testing the null hypothesis $H_{0}:C=1_{K^{(1)}}1_{K^{(2)}}^{T}$. The results are in Figures S6, S7, and S8. Results are similar to results from the $K=6$ and $p=10$ setting in Figure 3; we note that the ARI performs similarly to the $G$-test for independence in all cases.

C.3.2 Additional values of $\Sigma^{(1)}$ and $\Sigma^{(2)}$

We return to the simulation set-up described in Section 4 with $K=3$, $\mu^{(l)}$ given by

[TABLE]

and for two choices of $\Sigma^{(l)}$:

• Choice 1: $\Sigma^{(1)}=\Sigma^{(2)}=\begin{pmatrix}2.25&0.5\\ 0.5&2.25\end{pmatrix}$,
• Choice 2: $\Sigma^{(1)}=\begin{pmatrix}2.25&0.5\\ 0.5&2.25\end{pmatrix}$ and $\Sigma^{(2)}=\mathrm{diag}(2.25,4)$.

To perform Step 1 of Algorithm 1, we use the package mclust (Scrucca and others, 2016) to fit Gaussian mixture models with a common dense covariance matrix (the “EEE” covariance structure in mclust) for $\Sigma^{(l)}$ given by choice 1 above, and to fit Gaussian mixture models with a common diagonal covariance matrix (the “EEI” covariance structure in mclust) for $\Sigma^{(l)}$ given by choice 2 above. We compare the performance of the pseudo likelihood ratio test of $H_{0}:C=1_{K^{(1)}}1_{K^{(2)}}^{T}$ at nominal significance level $\alpha$ = 0.05 to the $G$-test for independence for testing the null hypothesis $H_{0}:C=1_{K^{(1)}}1_{K^{(2)}}^{T}$. The results are in Figure S9(i) and Figure S9(ii). Results are similar to results from the $K=6$ and $p=10$ setting with $\sigma=4.8$ in Figure 3.

C.3.3 Model misspecification

We consider a simulation set-up which compares the performance of the pseudo likelihood ratio test under model misspecification. We generate data from model (2.1) – (2.2), with $\Pi=\frac{1-\delta}{K^{2}}1_{K}1_{K}^{T}+\frac{\delta}{K}I_{K}$ for $K=3$. In the $l$th data view, the observations are drawn from a finite mixture model for which the $k$th component in the mixture (corresponding to the $k$th cluster) is a bivariate Student’s $t$-distribution with location parameter $\mu^{(l)}_{k}$ and scale matrix $\Sigma=\begin{pmatrix}2.25&0.5\\ 0.5&2.25\end{pmatrix}$, where the mean matrices for the two data views are of the form (C.4). We fit Gaussian mixture models with a common covariance matrix (the “EEE” covariance structure in mclust), and again use $B=200$ permutation samples. The results are in Figure S9(iii). We compare the performance of the pseudo likelihood ratio test of $H_{0}:C=1_{K^{(1)}}1_{K^{(2)}}^{T}$ at nominal significance level $\alpha$ = 0.05 to the $G$-test for independence for testing the null hypothesis $H_{0}:C=1_{K^{(1)}}1_{K^{(2)}}^{T}$, with p-values obtained with a permutation approach. Results remain similar to results from the $K=6$ and $p=10$ setting with $\sigma=4.8$ in Figure 3.