arXiv-1501.05194
companion code for manuscript http://arxiv.org/abs/1501.05194
view repo
The use of mutual information as a similarity measure in agglomerative hierarchical clustering (AHC) raises an important issue: some correction needs to be applied for the dimensionality of variables. In this work, we formulate the decision of merging dependent multivariate normal variables in an AHC procedure as a Bayesian model comparison. We found that the Bayesian formulation naturally shrinks the empirical covariance matrix towards a matrix set a priori (e.g., the identity), provides an automated stopping rule, and corrects for dimensionality using a term that scales up the measure as a function of the dimensionality of the variables. Also, the resulting log Bayes factor is asymptotically proportional to the plug-in estimate of mutual information, with an additive correction for dimensionality in agreement with the Bayesian information criterion. We investigated the behavior of these Bayesian alternatives (in exact and asymptotic forms) to mutual information on simulated and real data. An encouraging result was first derived on simulations: the hierarchical clustering based on the log Bayes factor outperformed off-the-shelf clustering techniques as well as raw and normalized mutual information in terms of classification accuracy. On a toy example, we found that the Bayesian approaches led to results that were similar to those of mutual information clustering techniques, with the advantage of an automated thresholding. On real functional magnetic resonance imaging (fMRI) datasets measuring brain activity, it identified clusters consistent with the established outcome of standard procedures. On this application, normalized mutual information had a highly atypical behavior, in the sense that it systematically favored very large clusters. These initial experiments suggest that the proposed Bayesian alternatives to mutual information are a useful new tool for hierarchical clustering.
READ FULL TEXT VIEW PDFcompanion code for manuscript http://arxiv.org/abs/1501.05194
Cluster analysis aims at uncovering natural groups of objects in a multivariate dataset [see Jain (2010) for a review]. In the vast variety of methods used in cluster analysis, an agglomerative hierarchical clustering (AHC) is a generic procedure that sequentially merges pairs of clusters that are most similar according to an arbitrary function called similarity measure, thereby generating a nested set of partitions, also called hierarchy (Duda et al., 2000). The choice of the similarity measure indirectly defines the shape of the clusters, and thus plays a critical role in the clustering process. While this choice is guided by the features of the problem at hand, it is also often restricted to a limited number of commonly used measures, such as the Euclidean distance or Pearson correlation coefficient (D’haeseleer, 2005). In the present work, we focus on the clustering of random variables based on their mutual information, which has recently gained in popularity in cluster analysis, notably in the field of genomics (Butte and Kohane, 2000; Zhou et al., 2004; Dawy et al., 2006; Priness et al., 2007) and in functional magnetic resonance imaging (fMRI) data analysis (Stausberg et al., 2009; Benjaminsson et al., 2010; Kolchinsky et al., 2014). Mutual information is a general measure of statistical dependency derived from information theory (Shannon, 1948; Kullback, 1968; Cover and Thomas, 1991). A key feature of mutual information is its ability to capture nonlinear interactions for any type of random variables (Steuer et al., 2002); also of interest, it indifferently applies to univariate or multivariate variables and can thus be applied to clusters of arbitrary size. Yet, mutual information is an extensive measure that increases with variable dimensionality. In addition, , the finite-sample estimator of mutual information, suffers from a dimensionality-dependent bias (see Appendix A). Several authors have proposed to correct mutual information for dimensionality by using a “normalized” version of mutual information (Li et al., 2001; Kraskov et al., 2005; Kraskov and Grassberger, 2009). In the clustering literature, normalized mutual information is routinely used. However, the impact of such correction procedure has not been extensively evaluated so far.
In the present paper, we consider Bayesian model-based clustering (Scott and Symons, 1971; Binder, 1981; Heller and Ghahramani, 2005; Jain, 2010) as an alternative to mutual information for the hierarchical clustering of dependent multivariate normal variables. Specifically, we derive a similarity measure by comparing two models: where and are independent (i.e., the covariance between any element of and any element of is equal to zero), against where the covariance between and can be set to any admissible value. The proposed similarity measure is then the log Bayes factor in favor of against (Kass and Raftery, 1995).With appropriate priors on the model parameters, we show that the similarity measure between and can be expressed in closed form. As will be shown below, the Bayesian formulation naturally (1) allows for clustering even when the sample covariance matrix is ill-defined; (2) provides for an automated stopping rule when the clustering reached has for any pair of remaining clusters; (3) corrects for dimensionality using a term that scales up the measure as a function of the dimensionality of the variables; and (4) provides for a local and global measure of similarity, in that it can be used to decide which pair of variables to cluster at each step (local level) as well as to compare different levels of the resulting hierarchy (global level). Asymptotically (i.e., when the number of samples ), the similarity measure is a linear function of mutual information, with a penalization factor that is in agreement with the Bayesian information criterion (BIC) (Schwarz, 1978). In this sense, the present paper makes an explicit connection between Bayesian model comparison for the clustering of dependent random variables and mutual information. The code corresponding to the Bayesian approach is freely available online^{1}^{1}1https://github.com/SIMEXP/arXiv-1501.05194/releases/tag/1.0.
We evaluated an AHC procedure based on this approach with synthetic datasets. The experiment aimed to evaluate how it behaved under both its exact and asymptotic forms compared to other approaches, including raw and normalized mutual information. We finally tested the new measures on two real datasets: a toy dataset and functional magnetic resonance imaging (fMRI) data.
In the following, we develop a Bayesian solution to the problem of clustering detailed above. We first introduce the model together with the Bayesian framework and a general expression for the similarity measure. In subsequent subsections, we derive a closed-form expression for the marginal model likelihoods under both assumptions of dependence and independence as well as exact and asymptotic expressions for the similarity measure. We then provide a description of the hierarchical agglomerative clustering algorithm resulting from the present development. We examine how the same framework can be conveniently used to compare nested partitions, that is, different levels of a hierarchy. We also deal with the issue of setting the hyperparameters. Finally, we show how the Bayesian solution can naturally provide for an automatic stopping rule.
Let be a -dimensional multivariate normal variable with known mean^{2}^{2}2For the sake of simplicity, we assumed a known mean in the following theoretical development. If the mean is unknown (as will be the case in the simulation and real data sections), this development is still valid, with replaced by and by
we propose to quantify the similarity between and as the log Bayes factor, that is, the log ratio of the marginal model likelihoods of versus
(1) |
Each marginal model likelihood can be expressed as an integral over the model parameters as described below.
Let us first calculate , the marginal model likelihood under the hypothesis of dependence. Expressing this quantity as a function of the model parameters yields
(2) |
Calculation of the integral requires to know the likelihood and the prior distribution of the covariance matrix under . With multivariate normal data, given is Wishart distributed with degrees of freedom and scale matrix (Anderson, 2003, Corollary 7.2.2), leading to the following likelihood
(3) |
where is the inverse of the normalization constant
As to the prior distribution, this quantity is here set as a conjugate prior, namely an inverse-Wishart distribution with
degrees of freedom and inverse scale matrix (Gelman et al., 2004, §3.6)(4) |
Bringing Equations (3) and (4) together into Equation (2) yields for the marginal model likelihood
The integrand is proportional to an inverse-Wishart distribution with degrees of freedom and scale matrix , leading to
(5) |
We can now calculate , the marginal model likelihood under the hypothesis of independence. If holds, then is block-diagonal with two blocks and the submatrix restrictions of to and , respectively. Introduction of the model parameters therefore yields for the marginal likelihood
(6) |
To calculate this integral, we again need to know the likelihood and the prior distribution of the two blocks of the covariance matrix under . The likelihood is the same as for and has the form of Equation (3). Furthermore, since is here block diagonal, we have and , where and are the restrictions of to and , respectively. Consequently, the likelihood can be further expanded as
(7) |
As to the prior distribution, assuming no prior dependence between and yields
(8) |
For the sake of consistency, and are set equal to and , respectively, which are in turn obtained by marginalization of as given by Equation (4). For , can be found to have an inverse-Wishart distribution with degrees of freedom and inverse scale matrix the restriction of to (Press, 2005, §5.2)
(9) |
Incorporating Equations (7), (8), and (9) into Equation (6) yields
Each integrand is proportional to an inverse-Wishart distribution with degrees of freedom and scale matrix , leading to
(10) |
Let us now express the Bayesian similarity measure by incorporating Equations (5) and (10) into Equation (1), yielding
(11) |
with
Another form for is
(13) |
with
(14) |
and
(15) |
quantifies, up to a constant that cancels out in , the amount by which the data support a model of multivariate normal distributions with unrestricted covariance matrix for .
We can now provide an asymptotic expression for . Define as the standard sample covariance matrix, i.e., . When , we can use Stirling approximation (Abramowitz and Stegun, 1972, p. 257) to expand the expression of Equation (15), leading to (see Appendix B)
(16) | |||||
where
is the plug-in estimator of mutual information for a multivariate normal distribution. Alternatively, can be seen as the minimum discrimination information for the independence of and (Kullback, 1968, Chap. 12, §3.6).
A hierarchy on a set of variables is a nested set of partitions , where is a partition of into clusters (Duda et al., 2000). A hierarchical agglomerative clustering (AHC) is a generic procedure to generate such a hierarchy, outlined in pseudo-code in Algorithm 1. The main steps of the algorithm are: (1) Initialize the partition with singletons (line 2); (2) derive a matrix where each element represents the similarity between two clusters and of , based on an arbitrary function (line 4); (3) identify the two clusters that are most similar^{3}^{3}3There may be more than one pair of clusters which maximize the similarity function. Most implementations of AHC proceed by selecting arbitrarily one such pair (e.g., the first one to be detected). In the in-house implementation we used, the pair was selected randomly amongst all these pairs. This was done to properly capture the instability of the algorithm. In such a form, AHC may not be deterministic anymore, in that two runs of the same algorithm on the same dataset may result in different hierarchies. (line 5); (4) form a new partition identical to the previous one, except that the two most similar clusters are now merged (lines 6–7); (5) iterate Steps 2–4 until the partition has only one single element which covers the whole set of variables (line 3). In the case of the methods proposed here, the similarity measure is given by either Equation (11) for the exact formulation or Equation (29) for the asymptotic BIC approximation.
The last development aims at providing a way to compare nested partitions, i.e., different levels of the hierarchy. Once the hierarchical clustering has been performed, each step is associated with a partition of . Assume that, at level , the partition reads and that, at step , and are clustered. Denote by the assumption that the partition at level is the correct partition of . The global improvement brought by the clustering of and between steps and can be quantified by the log ratio between the marginal model likelihoods
where both quantities can be computed in a manner similar to what was done for the similarity measure. For instance, if is true, then is block-diagonal with blocks ’s, the submatrix restrictions of to . Introducing the model parameters then yields
(17) |
In a way similar to what was done previously, the likelihood can be expanded as
(18) |
Turning to the prior distribution and assuming no prior dependence between the ’s, we can set
(19) |
Each can be obtained by the marginalization of , which is here taken as an inverse-Wishart distribution with degrees of freedom and inverse scale matrix the -by- diagonal matrix . Note that such a prior on is compatible with the prior used earlier for if one sets and if is the restriction of to (Press, 2005, §5.2). We then have
(20) |
Incorporating Equations (18), (19), and (20) into Equation (17) and integrating leads to
(21) |
The same calculation can be performed for . The result is the same as in Equation (21), except that the product is composed of terms. Of these terms, correspond to clusters that are unchanged from to and, as a consequence, are identical to those of Equation (21). The th term corresponds to the cluster obtained by the merging of and . As a consequence, the log Bayes factor reads
But this quantity is nothing else than . In other words, we proved that
(22) |
i.e., , the local measure of similarity between and , can be used to compute the global
measure of relative probability between two successive levels of the hierarchy.
Hyperparameter selection is often a thorny issue in Bayesian analysis. We here considered two approaches. The first approach (coined BayesCov) is to set the degree of freedom to the lowest value that still corresponds to a well-defined distribution, that is , and a diagonal scale matrix that optimizes the marginal model likelihood of Equation (21) before any clustering (that is, with clusters and for all ), yielding (see Appendix C)
where and are the diagonal elements of the prior inverse scale matrix and sum-of-square matrix , respectively. An alternative approach (coined BayesCorr) is to work with the sample correlation matrix instead of the sample covariance matrix. One can then set the number of degrees of freedom to
and the scale matrix to the identity matrix. The corresponding prior distribution yields uniform marginal distributions for the correlation coefficients
(Barnard et al., 2000). Note that clustering with the asymptotic form of Equation (29) (coined Bic) does not involve hyperparameters; it is also insensitive to the fact that the input is the covariance matrix or the correlation matrix.An advantage of the Bayesian clustering scheme proposed here and its BIC approximation is that they come naturally with an automatic stopping rule. By definition of in Equation (1), the fact that for the pair that is selected for clustering also means that the marginal model likelihood for is larger than that for . As such, and are more likely to belong to the same cluster than not and, as a consequence, it indeed makes sense to cluster them. By contrast, if we have for the same pair, the marginal model likelihood for is smaller than that for . and are therefore more likely to belong to different clusters. If, at a given step of the clustering, the pair with highest similarity measure has a negative similarity measures, then all pairs do, meaning that all pairs of clusters tested more probably belong to different clusters. It therefore makes sense to stop the clustering procedure at that point. This shows that an automatic stopping rule can simply be implemented into the clustering algorithm: Stop the clustering if the pair selected for clustering at a given step has . Note that, according to Equation (22), the resulting clustering corresponds to the one in the hierarchy that has largest marginal likelihood. We will refer to BayesCovAuto, BayesCorrAuto and BicAuto when applying the clustering schemes with this automated stopping rule.
To assess the behavior of the method expounded here, we examined how it fared on synthetic data. We used the two variants of the Bayes factor mentioned above (BayesCov and BayesCorr), Bic, as well as the same methods with automatic stopping rule (BayesCovAuto, BayesCorrAuto and BicAuto). As a mean of comparison, we also used the following methods—
A random hierarchical clustering scheme, where variables were clustered uniformly at random at each step. This category contains only one algorithm: Random, which was implemented for the purpose of the present study.
Hierarchical clustering schemes with similarity measures given by either Pearson correlation or absolute Pearson correlation, and a merging rule based on either the single, average, or complete linkage, or using Ward criterion. This category contains 8 algorithms: Single, Average, Complete, Ward, SingleAbs, AverageAbs, CompleteAbs, WardAbs. We used the implementations of these methods proposed in NIAK^{4}^{4}4https://github.com/SIMEXP/niak
Hierarchical clustering with a similarity measure given by mutual information, with and without normalization. This category contains 2 algorithms: Infomut and InfomutNorm. These methods were implemented for the purpose of the present study. Note that neither algorithm can run on small samples.
An approach where the clusters are estimated as the blocks of the precision (i.e., inverse covariance) matrix estimated with the graphical lasso—essentially a maximum-likelihood with -norm penalization (Friedman et al., 2008). The penalization parameter was set in by step of 0.1, then to 5, 10, 20, 50, 100, 200 (Smith et al., 2011). A version that optimizes with a BIC criterion was also used (Lian, 2011). Since the graphical lasso is not invariant by transformation of a covariance matrix into a correlation matrix, we used either the covariance matrix or the correlation matrix as input. Note that this approach automatically determined a number of clusters. Also, for (unconstrained case), the algorithm cannot run on small samples. This category contains 34 algorithms: 16 algorithms gLassoCov-x and 16 algorithms gLassoCorr-x, where x is the value of , and 2 algorithms gLassoCovOpt and gLassoCorrOpt. For this category of algorithms, we used a package freely available^{5}^{5}5http://www.cs.ubc.ca/schmidtm/Software/L1precision.html and already used in (Smith et al., 2011).
An approach based on the spectral clustering
(von Luxburg, 2006) of either the raw value or the absolute value of either the correlation or the partial correlation matrix. This approach required the number of clusters as input. Since this clustering has a step of -means, which is stochastic by nature, we considered 2 variants: one with 1 step of -means and the other one with 10 repetitions of -means and selection of the best clustering in terms of inertia. The similarity measures were defined so that the range would be the same (namely in ) when using the signed or the absolute value of correlation: and , respectively. This category contains 8 algorithms: 2 algorithms SpecCorr-x, 2 algorithms SpecCorrAbs-x, 2 algorithms SpecCorrpar-x and 2 algorithms SpecCorrparAbs-x, where x is the number of times that -means is performed. The spectral clustering algorithms of this category were coded for the purpose of the present study, while we used the implementations of the -means algorithm proposed in NIAK.All in all, 59 variants were tested.
In order to assess the performance of the Bayesian approach, we performed the following set of simulations. For each value of in , we considered partitions with increasing number of clusters (). For a given value of , we performed 500 simulations. For each simulation, the variables were randomly partitioned into clusters, all partitions having equal probability of occurrence (Nijenhuis and Wilf, 1978, chap. 12); (Wilf, 1999). For a given partition of , we generated data according to
(23) |
where all ’s were taken either as multivariate normal distributions (parameters: mean and covariance matrix ) or multivariate Student- distributions (parameters: degres of freedom , location parameter and scale matrix ). In both case, the ’s were set to and the ’s were first sampled according to a Wishart distribution with degrees of freedom and scale matrix the identity matrix and then rescaled to a correlation matrix. The sampling scheme on generated correlation matrices with uniform marginal distributions for all correlation coefficients (Barnard et al., 2000). For the multivariate Student- distributions, was set in . Equation (23) was used to generate synthetic datasets of length varying from 10 to 300 by increment of 40. Each dataset was summarized by its sample correlation matrix and hierarchical clustering was performed using the methods mentioned above. All simulations were implemented using the Pipeline System for Octave and Matlab, PSOM^{6}^{6}6https://github.com/SIMEXP/psom (Bellec et al., 2012) under Matlab 7.2 (The MathWorks, Inc.) and run on a 24-core server.
To assess the efficiency of the various methods, we thresholded each clustering at the true number of clusters, except for BayesCovAuto, BayesCorrAuto, BicAuto and gLasso for which we used the clustering obtained by the method. We then quantified the quality of the resulting partition using the proportion of correct classifications as well as the adjusted Rand index, which computes the fraction of variable pairs that are correctly clustered corrected for chance (Rand, 1971; Hubert and Arabie, 1985). Results corresponding to a given dimension and a given method were then pooled across numbers of clusters , lengths
and distributions (multivariate normal and Student). We classified the methods from best to worst based on these global results using the following indices (in this order): median of adjusted Rand index, 25% percentile value of adjusted Rand index, 5% percentile value of adjusted Rand index, smallest value of adjusted Rand index, and proportion of correct classifications. Note that some algorithms (
Infomut, InfomutNorm and SpecCorrpar) require the sample covariance matrix to be definite positive. As a consequence, these algorithms could not run on small samples. We therefore restrained our evaluation to cases where all algorithms were operational. Finally, we performed a Bayesian ANOVA-like regression analysis
(Gelman et al., 2004, §15.6), where we explained the adjusted Rand index of nine algorithms (BayesCov, BayesCovAuto, BayesCorr, BayesCorrAuto, Bic, BicAuto, Infomut, InfomutNorm, and AverageAbs) with the following effects: clustering algorithm (9 levels), number of variables (4 levels: ), type of distribution (4 levels: multivariate Gaussian or multivariate Student- with 1, 3, or 5 degrees of freedom), number of samples (8 levels: ), and number of clusters (68 levels: from 2 to clusters for each ). In other words, the model considered was of the form(24) |
Interactions between effects, while potentially relevant, were not considered to keep the analysis tractable. The posterior distribution for the various regression parameters were estimated using Gibbs sampling.
The results corresponding to the adjusted Rand index and fraction of correct classification are summarized in Figs 1 and 2 for the 20 best methods. Globally, and as expected, all methods were adversely affected by an increase in the number of variables. In all cases, the variants proposed in the present paper performed very well compared to other methods. BayesCov and BayesCorr were always classified as the two best algorithms and Bic was never outperformed by a method already published. The methods with automatic thresholding of the hierarchy performed surprinsingly well, considered that they were compared against methods with oracle. In particular, they clearly outperformed all variants of gLasso, the only method that was able to automatically determine the number of clusters. Of note, all variants of gLasso proved too slow to be applied to our simulation data for .
The results of the regression analysis are represented in Fig 3. The 9 algorithms selected included the ones proposed in the present manuscript (BayesCov, BayesCovAuto, BayesCorr, BayesCorrAuto, Bic, and BicAuto), the best-performing classic algorithm in the previous analysis (AverageAbs) as wel as the algorithms based on mutual information (Infomut and InfomutNorm). We found that increasing the dataset size () increased the performance of the algorithm, while increasing the dimensionality of the problem () and the number of clusters () decreased it. Note that dimension was found to have a negative influence on the adjusted Rand index, even though this index was partly proposed as a modification from the raw Rand index to overcome this limitation. Finally, this analysis confirmed the superior behavior of the methods proposed here, even when the method included the automatic stopping rule.
This data set was first used in Roverato (1999) and later re-analyzed in Marrelec and Benali (2006) in the context of conditional independence graphs. It originates from a study investigating early diagnosis of HIV infection in children from HIV positive mothers. The variables are related to various measures on blood and its components: and immunoglobin G and A, respectively; the platelet count; , lymphocyte B and T4, respectively; and the T4/T8 lymphocyte ratio. The sample summary statistics are given in Table 1. Roverato (1999) found that the correlations between and any other variable had strong probability around zero and hypothesized that the model was overparametrized. Based on the simulation study, we performed clustering of the data with the following methods: BayesCov(Auto), BayesCorr(Auto), Bic(Auto), Infomut, InfomutNorm, SingleAbs, AverageAbs, CompleteAbs, WardAbs, SpecCorrAbs and SpecCorrparAbs. For spectral clustering, we used either 1 or 10 repetitions of the -means step; since the results obtained for 1 step of -means were highly variable for 3, 4, and 5 clusters, we discarded these results.
8.8374 | 0.479 | 0.356 | ||||
0.483 | 0.1919 | 0.068 | ||||
0.220 | 0.057 | 8924231.9 | 0.085 | 0.552 | ||
0.149 | 20392.4 | 0.091 | 0.013 | |||
0.253 | 0.523 | 0.179 | 1952795.2 | 0.384 | ||
0.064 | 0.213 | 1.378 | ||||
Summary statistics for the HIV data: sample variances (main diagonal), correlations (lower triangle) and partial correlations (upper triangle) [from
Roverato (1999)].The resulting clusterings are given in Fig 4 and Table 2. All hierarchical clusterings started by clustering and (lymphocite). This was confirmed by SpecCorrAbs-10 when required to provide 5 clusters. All hierarchical clustering methods then clustered and (immunoglobin). This result was in agreement with both SpecCorrAbs-10 and SpecCorrparAbs-10 when required to provide 4 clusters. After the second step, we observed four behaviors for the AHC algorithms and classified them accordingly:
BayesCov, BayesCorr, Infomut and InfomutNorm;
SingleAbs and AverageAbs;
Bic;
CompleteAbs and WardAbs.
While not a hierarchical clustering, SpecCorrAbs-10 provided successive clusterings that were in agreement with methods in (G2). Algorithms in (G1) and (G3) clustered with , creating a cluster of variables related to lymphocite. Algorithms in (G1) and (G2) (and SpecCorrAbs-10) agreed that a partitioning of the variables in two clusters should leads to on the one hand and on the other hand. This clustering was also found by SpecCorrpar-10 when constrained to extract two clusters from the data. It was also considered the best clustering for BayesCov and BayesCorr. For Bic, the optimal partitioning was composed of three clusters, namely, , , and , which is in agreement with what would methods in (G1) yield for three clusters; furthermore, it still kept separated from the other variables.
In Fig 5, we represented the evolution of the Bayes factor during hierarchical clustering for BayesCov, BayesCorr and Bic. Note that, while the clustering steps are identical for BayesCov and BayesCorr, the log Bayes factors are similar but not identical. Likewise, while the first two clustering steps of Bic is identical to those of BayesCov and BayesCorr, one can see that, unlike BayesCov and BayesCorr, Bic considered the merging of with almost as likely as that of with . Also, while the successive clusterings of with and then as well as that of with strongly increased the Bayes factor for BayesCov and BayesCorr, the improvement brought by the clustering of with in these methods was less important.
All in all, this analysis led us to the following conclusion: it is very likely that variables and belong to the same cluster of dependent variables, and similarly for variables and . Also, there is very strong evidence in favor of the fact that is independent from the other variables. Finally, we suspect that , and could belong to the same cluster of variables.
SpecCorrAbs-10 | SpecCorrparAbs-10 | |
---|---|---|
2 clusters | ||
3 clusters | ||
4 clusters | ||
5 clusters | not reproducible |
Cluster analysis is a popular tool to study the organization of brain networks in resting-state fMRI (Yeo et al., 2011; Kelly et al., 2012), by identifying clusters of brain regions with highly correlated spontaneous activity. We applied the 13 methods that were found to have good performance on simulations (see Fig 6) to a collection of resting-state time series. The time series had 205 time samples and were recorded for 82 brain regions in 19 young healthy subjects. See Appendix D for details on data collection and preparation. The data are available online^{7}^{7}7http://figshare.com/articles/Atlanta_resting_state_fMRI_time_series
_preprocessed_using_the_AAL_template/1521155.
We first aimed at establishing which clustering algorithms yielded similar results on these datasets. We more specifically investigated a 7-cluster solution, as this level of decomposition has been examined several times in the literature (Bellec et al., 2010; Power et al., 2011; Yeo et al., 2011). Each clustering algorithm was applied to the time series of each subject independently. For a given pair of methods, the similarity between the cluster solutions generated by both methods on the same subject was evaluated with the Rand. Note that the raw, unnormalized Rand index was used here, as we did not compare cluster solutions with different numbers of clusters, which is the main motivation of the normalization. The unnormalized Rand index has a more intuitive interpretation than its adjusted counterpart. The Rand indices were averaged across subjects, hence resulting into a method-by-method matrix capturing the (average) similarity of clustering outputs for each pair of methods (see Fig 6). An AHC with Ward’s criterion was applied to this matrix in order to identify clusters of methods with similar cluster outputs. We visually identified five clusters of methods that had high () average within-cluster Rand index. The largest cluster included classical AHCs such as CompleteAbs, AverageAbs, WardAbs, as well as the Bic and BayesCov methods proposed here. It should be noted that this class of algorithms generated solutions for this problem that were very close to one another (average within-cluster Rand index ). The BayesCorr method was also close to that large group of methods, but not quite as much as the aforementioned methods (average Rand index of about 0.7), and was thus singled out as a separate cluster. The spectral methods were split into two clusters, depending on whether they were based on correlation (SpecCorrAbs-1 and SpecCorrAbs-10) or partial correlation (SpecCorrparAbs-1 and SpecCorrparAbs-10). Finally, the two variants of mutual information (Infomut and InfomutNorm) generated solutions that were highly similar to SingleAbs. It was reassuring that the variants of Bayes methods proposed here performed similarly to algorithms known to produce physiologically plausible solutions, such as Ward (Thirion et al., 2014; Orban et al., in press). While we found some analogy between BayesCorr, BayesCov, Infomut and InfomutNorm, it was intriguing that the variants of mutual information tested seemed to generate markedly different classes of solutions from the Bayes methods. We decided to examine the cluster solutions of these methods in more details.
As a reference, we examined the cluster solutions generated by WardAbs, in addition to two variants of Bayes clustering that yielded slightly different solutions (BayesCov and BayesCorr), and nomalized mutual information, InfomutNorm. To represent the typical behavior of each method across subjects, we generated a “group” consensus clustering summarizing the stable features of the ensemble of individual cluster solutions. This consensus clustering was generated by the evidence accumulation algorithm (Fred and Jain, 2005) outlined below. First, each partition of each subject was represented as a binary 81-by-81 adjacency matrix , where for each pair of brain regions, if areas and were in the same cluster, and otherwise. The adjacency matrices were then averaged across subjects and selected methods, yielding a 81-by-81 stability matrix where each element coded for the frequency at which brain areas and fell in the same cluster. Finally this stability matrix was used as a similarity matrix in a AHC with Ward’s criterion to generate one consensus partition. The brain regions were grouped into the same consensus cluster if they had a high probability of falling into the same cluster on average across subjects and methods, hence the name consensus clustering.
Fig 7 represents the stability matrices and consensus clusters, for the four methods of interest. As expected based on our first experiment on the similarity of cluster outputs, the WardAbs and BayesCov methods were associated with highly similar stability matrices and almost identical consensus clusters. Many areas of high consensus could be identified (with values close of 0 or 1), illustrating the very good agreement of the cluster solutions across subjects. The outline of the consensus clusters as well as a volumetric representation of the brain parcellation are presented in Fig 7b. Some of these clusters closely matched those reported in previous studies: cluster 7 can be recognized as being the visual network, cluster 2 the sensorimotor network, and clusters 6 and 3 the anterior and posterior parts of the default-mode network, respectively (Salvador et al., 2005; van den Heuvel et al., 2008; Bellec et al., 2010). By contrast with WardAbs and BayesCorr, the InfomutNorm tended to generate very large clusters, which was apparent both on the stability matrix and the consensus clusters. The BayesCorr method was intermediate between BayesCov and InfomutNorm in terms of cluster size. These decompositions into very large clusters do not fit current views on the organization of resting-state networks.
Overall, our analysis on real fMRI data led to the following conclusions. Three big subsets of methods emerged: spectral methods, mutual information (with SingleAbs), and finally all the other methods. Application of this last group of methods, which included the Bayes variants proposed here, resulted in a plausible decomposition into resting-state networks. In the absence of ground truth, it is not possible to further comment on the relevance of the differences in cluster solutions identified by the three groups of methods. We still noted that the methods based on mutual information led to large clusters that were difficult to interpret. Our interpretation is that the strategies implemented in Infomut and InfomutNorm did not behave well for these datasets.
As a final computational note, the time required by the different methods to cluster data is summarized in Table 3. Although the differences in execution time may reflect the quality of the implementation, the methods proposed here were the slowest of the hierarchical methods, but were still faster than spectral clustering.
method | minimum | median | maximum |
---|---|---|---|
CompleteAbs | 10.9 ms | 11.5 ms | 24.5 ms |
AverageAbs | 11.0 ms | 11.6 ms | 25.2 ms |
SingleAbs | 11.1 ms | 11.7 ms | 49.5 ms |
WardAbs | 14.1 ms | 14.8 ms | 18.9 ms |
InfomutNorm | 159 ms | 170 ms | 261 ms |
InfoMut | 299 ms | 322 ms | 416 ms |
Bic | 666 ms | 673 ms | 810 ms |
BayesCov | 1.083 s | 1.094 s | 1.263 s |
BayesCorr | 1.108 s | 1.151 s | 1.133 s |
SpecCorrparAbs-1 | 1.928 s | 2.065 s | 2.291 s |
SpecCorrAbs-1 | 1.992 s | 2.176 s | 2.417 s |
SpecCorrparAbs-10 | 13.251 s | 13.939 s | 14.301 s |
SpecCorrAbs-10 | 14.456 s | 15.381 s | 16.101 s |
We here proposed some novel similarity measures well suited for the agglomerative hierarchical clustering of dependent variables. These measures rely on a Bayesian model comparison for multivariate normal random variables. On synthetic data with a known (ground truth) partition, hierarchical clustering based on the Bayesian measures was found to outperform several standard clustering procedures in terms of adjusted Rand index and classification accuracy. On the toy example, the Bayesian approaches led to result similar to those of mutual information clustering techniques, with the advantage of an automated thresholding. On the real fMRI data, the Bayesian measures led to results consistent with standard clustering methods, in contrast to methods based on mutual information, which exhibited a highly atypical behavior.
A key feature of the Bayesian approach is its ability to take the dimension of the clusters into account. Dimensionality is an important issue in two respects (see Appendix A for an illustration). First, mutual information is an extensive measure that depends on the dimension of the variables. This has motivated the introduction of a normalization factor in the application of mutual information to hierarchical clustering (Kraskov et al., 2005; Kraskov and Grassberger, 2009)
. A second issue is the existence of a bias in the estimation of mutual information. This bias mechanically increases with the dimensionality of the variables. Because of the two issues described above, hierarchical clustering based on mutual information will tend to cluster unrelated but large variables rather than correlated variables of lower dimensions. As demonstrated on real fMRI data, the heuristic proposed by
(Kraskov et al., 2005) does not provide a general solution to the issue of dimensionality. Furthermore, it removes some interesting features of mutual information, such as the additivity of the clustering measure. By contrast, the Bayesian approach takes the dimensionality into account in a principled way, providing a quantitative version of Occam’s factor (Jaynes, 2003, Chap. 20). The Bayesian normalization is additive rather than multiplicative, thus preserving the additive properties of mutual information.We defined the similarity measure between any two pairs of variables and as a log Bayes factor. At each step, the pair that had the largest similarity measure was merged. Taking into account the unique features of as the log Bayes factor defined in Equation (1) allowed us to have access to a global measure of fit as defined in Equation (22) as well as to provide an automatic stopping rule that behaved very well on simulated data. Going from a similarity measure to a log Bayes factor has other advantages that could take the clustering proposed here even further (see below).
The new alternatives to mutual information introduced in this paper (i.e., Bayes and Bic) proved useful for the analysis of resting-state fMRI. The benefits were particularly clear when compared to InfomutNorm, which tended to create large, inhomogeneous clusters. By contrast, both Bayes and Bic had a behavior similar to standard hierarchical clustering based on Pearson’s linear correlation. The possible benefits of Bayes and Bic over those canonical methods are still substantial. The mutual information first provides a multivariate measure of interaction that is well adapted to hierarchical brain decomposition (Tononi et al., 1994; Marrelec et al., 2008) and which has a clear interpretation in information theory (Watanabe, 1960; Joe, 1989; Studený and Vejnarová, 1998). For these reasons, the mutual information for Gaussian variables is more appealing than a simple average of pairwise correlation coefficients across clusters. Because mutual information is measured between clusters, it is natural to build the clusters themselves based on this metric. A second benefit of the proposed approach is that Bic proved to be the most stable of all tested methods in the range of 5–15 clusters on real fMRI datasets. The increase in stability over Ward’s was modest, but significant. This advantage may become even more substantial if the clustering is performed in higher dimension, i.e., with smaller areas than the AAL brain parcellation or even at the voxel level.
Clustering techniques are based on either a similarity measure or a distance measure. While the description of the present manuscript mostly relied on the notion of similarity, going from one concept to the other one can generally be done with minor changes. For instance, standard hierarchical procedures which rely on the minimization of a distance to perform clustering (e.g., single, complete and average linkage) can be applied to cases where closeness is quantified by a measure of similarity, simply by using the opposite of the similarity matrix as a distance matrix. Although the resulting measure may not define a mathematically valid distance, it is not required for the procedure to work. Similarly, in a Bayesian framework, switching from a similarity measure to a distance measure tantamounts to switching from
to
that is, from the log ratio of the marginal model likelihoods in favor of dependent variables to the log ratio of the marginal model likelihoods in favor of independent variables.
Any Bayesian analysis requires the introduction of prior distributions. In the present study, we needed the prior distribution for the covariance matrix associated to any clustering of . Our choices were guided by one assumption, one rule of consistency, and one rule of simplicity. First, our assumption was to not assume a priori any sort of dependence between covariance matrices associated to different clusters. This allowed for the decomposition of any prior as a product of independent priors, see Equations (8) and (19). The rule of consistency was to consider that the prior for a given covariance matrix should not be contradictory at different levels of the hierarchy. This is why we set the prior distribution for the global covariance matrix as an inverse-Wishart distribution and then derived the prior for any covariance matrix as the prior distribution for integrated over all parameters that do not appear in ; using such an approach, the resulting prior turned out to be an inverse-Wishart distribution as well. Last, the rule of simplicity is the one that dictated the use of inverse-Wishart distributions as prior distributions for the covariance matrices. Such a family of priors had the twofold advantage of being closed by marginalization and allowing for closed-form expressions of the quantities of interest. An inverse-Wishart distribution is characterized by two parameters: the degrees of freedom and the inverse scale matrix . If is a -by- matrix, we must have for the distribution to be normalized. Also, must be positive definite. From there, we had two strategies. The first one was to set the degree of freedom to the lowest value that still corresponded to a well defined distribution () and a diagonal scale matrix that optimized the marginal model likelihood. An alternative approach was to work with the sample correlation matrix, set and equal with the -by- identity matrix , since this choice corresponds to a distribution that is associated with uniform marginal distributions of the correlation coefficients (Barnard et al., 2000). While we believe that our assumption and the rule of consistency are sensible choices, we must admit that we are not quite as content with the choice of inverse-Wishart distributions for priors. The major issue with such a family of priors is that they simultaneously constrain the structure of correlation and the variances. By contrast, it seems intuitive that clustering should depend on the correlation structure only, not on the variances. As such, the prior on the variances should ideally be set independently from the correlation structure. Priors that separate variance and correlation have already been proposed (Barnard et al., 2000). Unfortunately, the use of such priors would make it impossible to provide a closed form for the marginal model likelihood. While numeric schemes could be implemented to circumvent this issue, it would render the procedure proposed much more complex and computationally burdensome. By contrast, the algorithm detailed here is rather straightforward. Also, the influence of the prior vanishes when the sample size increases. From a practical perspective, the three methods proposed here (two, BayesCov and BayesCorr, with different priors and one, Bic, not influenced by the prior) exhibited similar behaviors and still outperformed other existing methods in the simulation study. We take it as a proof of the robustness of the method to the choice of prior.
The presence of clusters of variables that are mutually independent is equivalent to having a covariance matrix that is block diagonal, which is itself equivalent to having a precision (or concentration) matrix that is also block diagonal. One could therefore solve the problem using precision matrices instead of covariance matrices. The corresponding calculations can be found in Appendix E. The main difference between the two approaches stems from the fact that a submatrix of the inverse covariance matrix is not equal to the inverse of the corresponding submatrix of the covariance matrix, that is, , unless is block diagonal; also, Wishart and inverse-Wishart distributions do not marginalize the same way. These differences are to be related to the fact that a submatrix of a covariance matrix is better estimated than the whole covariance matrix, while the same does not hold for a precision matrix. From there, we could expect the covariance-based approach to perform better than the precision-based approach, and the difference to increase with increasing . This was confirmed on our synthetic data, where the precision-based approach behaved as well as BayesCov and BayesCorr for but had worse results than these two approaches for . Besides, performance of the automated stopping rule was much more efficient with BayesCov and BayesCorr than it was with the precision-based approach. As a final note, basing the method on concentration matrices yielded a slower algorithm, arguably because of the matrix inversions that are required.
Intuitively, the structure of dependence of a multivarite normal distribution is included in its covariance matrix. All existing algorithms that we used here do not need the full dataset but only the sample covariance (or correlation) matrix. This is why we started with a likelihood model that only considers the covariance matrix [see Equations (3) and (7)]. Rigorously, this model is only valid for ; for , one should resort to a model of the full data as being multivariate normal with a mean and covariance matrix . Nonetheless, we kept the ’intuitive’ approach as it has fewer hyperparameters, is easier to deal with, and leads to formulas that are simpler to interpret. Also, the resulting similarity measure [Equation (11)] is not restricted to , but is well defined for as well. From a practical perspective, the ’intuitive’ algorithm only requires the sample covariance matrix as input, while a full model would also require the sample mean. Finally, this simpler model exhibited good behavior on our synthetic data, even for small datasets.
Regarding the computational cost, measures derived from mutual information or a Bayesian approach are more demanding than standard methods such as Average or Ward. The derivation of the similarity matrix and the search for the most similar pairs of clusters are the two critical operations that can be optimized to decrease computation time. It is always possible to speed up these two steps by taking advantage of the fact that the similarity matrices of two successive iterations of the algorithm have many elements in common, as all but one element of a partition at a given iteration are identical to the elements of the partition of the previous iteration. Critically, in the case of Average and Ward methods, it is in addition possible to derive the similarity matrix at every iteration only based on the similarity matrix at initialization through successive updates using the Lance-Jambu-Williams recursive formula (Batagelj, 1988). By contrast, other measures, including BayesCov, BayesCorr, Bic, Infomut, and InfomutNorm, have to be re-evaluated independently at every step, which means that the determinant of a covariance matrix with increasing size has to be computed. Finding an update formula analogous to Lance-Jambu-Williams for clustering methods based on variants of the mutual information would substantially accelerate these algorithms.
Another extension would be to replace the deterministic rule of selecting the pair with largest similarity measure for merging by a probabilistic rule where the probability to cluster a given pair is given by the posterior probability of the resulting clustering.
A last extension that could be easily implemented in the present framework is the generalization of the method to account for different entities (e.g., subjects in fMRI) sharing the same structure but with potentially different covariance matrices for that structure. If each entity is characterized by a variable and corresponding sample sum-of-square matrix , one can perform AHC on each using and the corresponding similarity measure . However, with a straightforward modification of the present method, one could also perform global AHC of all covariance matrices considered simultaneously. Assuming that the covariance matrices of the different elements are independent given the common underlying structure, then the resulting similarity measure is the sum of all individual similarity measures ’s.
Altogether, the Bayesian framework that we used provides a principled way to generalize our approach to distributions other than multivariate normal ones. Such generalization would potentially account for a wide variety of situations, such as nonlinear dependencies or discrete distributions. This would widen the scope of possible applications of the technique, e.g., genetics (Butte and Kohane, 2000; Zhou et al., 2004; Dawy et al., 2006; Priness et al., 2007)
. The issues related to this generalization are twofold. First, one needs a model of dependence. In the discrete case, one could think of multinomial distributions, originating from categorical, i.e., generalized Bernoulli distributions
(Papoulis, 1965, §3.4). In the continuous case, multivariate normal distributions are a first choice model beyond which it is not clear what to use. Multivariate Student- distributions could be considered, even though our results on simulated data would tend to hint that the difference with multivariate normal distriubtions might not be that large. One could also consider using models where dependence is controled independently of the marginal distributions, such as multivariate copulas (Nelsen, 1999; Fischer, 2011). Another issue is the possibility to express the marginal posterior likelihood of the data given the model selected. For multivariate discrete variables, we expect it to be feasible, albeit computationally very challenging and sensitive to the type of prior distribution. For other distributions, obtaining a closed form might be out of reach. Nonetheless, the marginal posterior likelihood could be approximated using various criteria, such as the AIC (Akaike, 1974) or variants thereof—AICc (Burnham and Anderson, 2004); (McQuarrie and Tsai, 1998, §2.3.1) or AICu (McQuarrie and Tsai, 1998, §2.4.1)—, or the BIC (Schwarz, 1978), which naturally appeared in the present derivation. In any case, any approach beyond multivariate normal distributions would drastically increase the complexity of our approach, both in terms of inference and computation.In the present manuscript, we used a hierarchical algorithm as a way to extract an underlying structure of dependence from data. Our assumption was that there was one such structure. Such an approach provided a simple and efficient clustering algorithm with an interesting connection to mutual information. However, the method as presented here is not able to deal with data that are truly organized hierarchically. Extending it to deal with such data would improve the scope of the algorithm. One way to do would be to use Dirichlet process mixtures (Heller and Ghahramani, 2005; Heller, 2007; Savage et al., 2009; Cooke et al., 2011; Darkins et al., 2013; Sirinukunwattana et al., 2013), together with a model of dependent variables.
In this paper, we proposed a procedure based on Bayesian model comparison to decide whether or not to merge Gaussian multivariate variables in an agglomerative hierarchical clustering procedure. The resulting similarity measure was found to be closely related to the standard mutual information, with some additional corrections for the dimensionality of the datasets. These new Bayesian alternatives to mutual information turned out to be beneficial to hierarchical clustering on Gaussian simulations and real datasets alike. Because of the simplicity of its implementation, its good practical performance and the potential generalizations to other types of random variables, we believe that the approach presented here is a useful new tool in the context of hierarchical clustering.
The authors are grateful to the members of the 1000 functional connectome consortium for publicly releasing the ’Atlanta’ data sample, to the UNF (Unité de neuroimagerie fonctionnelle, Université de Montréal, Montréal, Qc, Canada) for providing them with computational resources, to Mathieu Desrosiers (UNF) for his expertise and suppport running the simulations on the grid engine. Part of this study was performed while G.M. was at the UNF.
G.M. and P.B. conceived and designed the experiments, performed the experiments, and analyzed the data. G.M., A.M., and P.B. contributed reagents/materials/analysis tools and wrote the manuscript. G.M. conceived the theoretical model.
P.B. offers consulting in brain image analysis for two Contract Research Organizations, NeuroRX and Biospective, Montreal, Canada.
A.M. is supported by Deutsche Forschungsgemeinschaft (DFG) grant SFB 936/Z3. P.B. is supported by the Fonds de Recherche du Québec – Santé (FRQS). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript
In the case of a -dimensional variable with a multivariate normal distribution, the mutual information between two subvectors and is given by