Hierarchical Clustering and Dendrograms

The algorithm starts with fine units. It repeatedly combines the closest eligible groups. Every stage below follows from two decisions made before the first merge: the distance metric and the linkage rule. The tree records the joint history of those two choices. It does not record a structure guaranteed to exist in the world.

RNA-seq profiles from 15 tissues, human and mouse, grouped by species rather than by tissue. The Mouse ENCODE Consortium published that comparison in PNAS in December 2014, and the consortium's own record of the paper reproduces the finding: “the expression for many sets of genes was found to be more similar in different tissues within the same species than between species”. Read as biology, that is a striking claim about how far the two species have diverged.

The following May, an independent group at the University of Chicago reanalysed the same data. Gilad and Mizrahi-Man found that the assignment of samples to sequencing flowcells and lanes was nearly completely confounded with species. The machine run and the organism were, in practice, the same variable. Their abstract states the result plainly: “When we account for the batch effect, the corrected comparative gene expression data from human and mouse tend to cluster by tissue, not by species.” Same samples, once corrected for the sequencing batch: opposite conclusion. Nothing in the shape of the first dendrogram announced which of the two it was drawing.

A hierarchy makes uncertainty look orderly. Its branches record the chosen linkage process applied to whatever the measurements actually vary with. That is not a guaranteed taxonomy of the world.

Cutting the tree discards some of the multiscale information that made hierarchical clustering attractive.

The dendrogram's dependence on linkage and metric is not a matter of taste either, to be settled by whichever pair looks best. Several linkages are defined only for one metric. SciPy 1.18.0's documentation states it as a restriction: “Methods 'centroid', 'median', and 'ward' are correctly defined only if Euclidean pairwise metric is used.” MATLAB's method table for linkage says the same in its own words. It marks 'centroid' as “Centroid distance (UPGMC), appropriate for Euclidean distances only”, 'median' as “Weighted center of mass distance (WPGMC), appropriate for Euclidean distances only” and 'ward' as “Inner squared distance (minimum variance algorithm), appropriate for Euclidean distances only”.

Both packages will nonetheless run Ward over a cosine distance matrix and return a tree. That tree is outside the criterion's definition. No amount of downstream cutting, naming or stability checking repairs it. The metric you chose decides which column of this comparison you are entitled to build at all.

Merge five invented points on a line, at 1, 2, 4, 8 and 9, step by step under two linkage rules.

These checks focus attention on evidence rather than decorative branch labels. Each one has a documented limit, and the limits are as useful as the checks.

Preserve more than one level until the intended use requires a flat decision. Two of these steps have published tooling and published benchmarks behind them, so neither has to rest on judgement alone.

Step 3, branch stability, can be computed rather than eyeballed. The R package pvclust, published by Suzuki and Shimodaira in Bioinformatics in 2006, attaches p-values to every cluster in a dendrogram: an approximately unbiased (AU) p-value from multiscale bootstrap resampling, and a bootstrap probability (BP) from the ordinary bootstrap. The AU value is superior in bias to BP. The current CRAN release, version 2.2-0 from November 2019, reports a third. Its package description: “An implementation of multiscale bootstrap resampling for assessing the uncertainty in hierarchical cluster analysis. It provides SI (selective inference) p-value, AU (approximately unbiased) p-value and BP (bootstrap probability) value for each cluster in a dendrogram.”

Step 4, the cut, has been benchmarked head-to-head on data with known answers. Milligan and Cooper ran a Monte Carlo evaluation of 30 stopping rules in 1985, on artificial data containing 2, 3, 4 or 5 distinct nonoverlapping clusters, analysed by four hierarchical methods. Their verdict: “The simulation results for the stopping rules revealed a wide range in their ability to determine the correct number of clusters in the data. Several procedures worked fairly well, whereas others performed rather poorly.” SAS/STAT's 14.1 user guide names the three that performed best with a high degree of error in the data: the pseudo F statistic of Caliński and Harabasz (1974), the Je(2)/Je(1) statistic of Duda and Hart (1973) that can be transformed into a pseudo t², and the cubic clustering criterion. Thirty candidate rules, a handful defensible. The number of clusters is a claim to be argued, not a value to be read off.

Many datasets support coarse regimes and finer local variants rather than one privileged partition. Hierarchical clustering can preserve that structure when its branches remain stable and interpretable. The tree should remain an analytical object, not an automatic ontology. Names and actions still require evidence beyond the merge procedure.

Even the linkage name can be ambiguous. Two algorithms in the literature both go by Ward's method, and Murtagh and Legendre set them side by side in 2014: “Two algorithms are found in the literature and software, both announcing that they implement the Ward clustering method. When applied to the same distance matrix, they produce different results. One algorithm preserves Ward’s criterion, the other does not.” R ships both, and its manual is explicit. The “ward.D” option “does not implement the clustering criterion of Ward Jr. (1963)”, whereas “ward.D2” does. The authors' own preprint tells readers what follows from that: check what their favorite software package is doing.

A tree is reproducible only when the software option is recorded too — alongside the metric, the tie-breaking rule, whether leaf ordering was optimized, and the stability statistic that supported the level you acted on.