Cluster analysis using the posterior distribution of a Bayesian correlation matrix Background and Problem
I recently ran a Bayesian multivariate epidemiological meta-analysis on prevalence estimates for several disorders. This analysis included a probit-based model to deal with the individual patient data (N ~ 1600) I was able to dig up. In addition to providing an estimate of the individual prevalences for each of the 8 disorders included in my model, I also end up with an 8 x 8 tetrachoric correlation matrix between the latent variables representing the disorders in probit-space. One of my collaborators pointed out that it would be interesting to use this correlation matrix to conduct a cluster analysis on the disorders themselves, to see which disorders clustered together. I am (vaguely) familiar with approaches like the hclust(), factanal() or principal() functions in the R programming language (or similar approaches in Python, etc.) that can take correlation matrices (or derivations thereof) and use them to estimate clusters, factors or components. However, a unique challenge posed by my data is that – being a Bayesian correlation matrix – I do not have a single estimate for the matrix, but rather 2000 correlation matrices representing the posterior distribution of that matrix. As such, I am unsure how best to proceed.
Question
I was hoping that someone could suggest a clustering approach that:


*

*is easily fit using a correlation matrix (or variation thereof);

*produces output that can somehow be combined across fits (in this case combining estimates across posterior samples);


If such a thing exists. I must admit, I am not even sure how to approach this one, so any advice concerning the situation would be most welcome. I am leaning towards saying it is simply not tenable.
 A: It is possible to do Bayesian cluster analysis using a distance matrix, see this paper. It's fairly common to use correlation as a distance measure for clustering (obviously you must use $1-\rho$ rather than $\rho$ as high correlation implies similarity).
As for a way to account for the uncertainty in posterior samples, I'm less sure. I would lean towards something like what Pedro suggested; to cluster each posterior sample and then do some form of averaging. However, given 2 levels of uncertainty, and the inherent variability of clustering algorithms, that could get very messy. I don't know that you would get sensible or interpretable results from the latter approach.
A: A few suggestions:  
The straight-forward way:  Take the posterior means/MAP of each of the pairwise correlations.  If you think there should exist an hierarchical relationship among the data-points then something like hierarchical clustering would make sense. If you believe a more flat structure should be present, then you could do a factor analysis and/or run PCA.   
the Bayesian way:    rather than make a single point estimate for your correlation matrix, run an analysis on each one of the posterior samples then take the average.  So, for instance, you could run PCA, see how many eigenvalues are important, run k-means on the important eigenvectors, and assign it to a cluster. Repeat for each sample (you'll have to do a bit of work to match the clusters labels across samples, since the ordering of the labels will be random each time.)  
Finally, after you have  2000 instances of class labels for each disorder, you can produce statistics such as probability that disorder A falls into cluster 1, etc.  
In a nutshell, for each sample you can assign a hard-assignment for each disorder. By averaging over all your 2000 samples, you end up with soft-labels/ mixed membership cluster analysis.  
