Skip to content

Clustering Methods

Clustering is an unsupervised method that groups similar observations together based on their feature values, without using predefined labels. The goal is to find structure in data where none has been imposed externally.

Two broad families of methods are covered here: hierarchical clustering, which builds a tree of nested groups, and k-means, which partitions data into a fixed number of groups. The two methods use similarity differently, but in both cases the definition of what counts as "similar" strongly affects the result.

Distance Metrics

Before clustering, you must define how dissimilarity between observations is measured. The right choice depends on your data type: Euclidean for scaled continuous variables, Bray-Curtis for species counts and OTU tables, correlation distance for expression profiles. A full reference with formulas and R code is on the Distance Metrics page.

Always scale continuous data before clustering

Variables on different scales will dominate distance calculations. Use scale() unless you are working with a distance metric designed for count data such as Bray-Curtis.

Hierarchical Clustering

Hierarchical clustering builds a tree (dendrogram) by iteratively merging the two most similar clusters. Starting with each observation as its own cluster, it continues until all observations are joined into one.

Linkage methods

The linkage method defines how distance between clusters is calculated as they grow:

Method Distance used Tends to produce
Complete Maximum between any two points Compact, spherical clusters
Single Minimum between any two points Elongated chains
Average Mean of all pairwise distances Compromise
Ward.D2 Minimises within-cluster variance Compact, relatively balanced clusters

Ward's method is often a good choice for scaled continuous biological data. It merges the two clusters whose fusion produces the smallest increase in total within-cluster variance (error sum of squares), rather than defining inter-cluster distance directly from pairwise distances as single, complete, or average linkage do. It tends to produce compact, relatively balanced clusters and often yields dendrograms with pronounced height differences, making it easier to choose a cut level.

Because both Ward's method and k-means minimise within-cluster variance, they often produce similar solutions when applied to Euclidean data; disagreement between them can indicate weak or ambiguous cluster structure.

A key limitation is that Ward's method is derived for Euclidean distances. When working with non-Euclidean dissimilarities such as Bray-Curtis, average linkage is typically preferred.

Implementation

library(factoextra)

# Scale and compute distances
df_scaled <- scale(iris[, 1:4])
dist_matrix <- dist(df_scaled, method = "euclidean")

# Cluster with Ward's method
hc <- hclust(dist_matrix, method = "ward.D2")

# Plot dendrogram
plot(hc, main = "Hierarchical Clustering", xlab = "", sub = "")

# Cut into k clusters
rect.hclust(hc, k = 3, border = 2:4)

# You can also cut by height instead of by number of clusters.
# On this dendrogram, h = 15 merges versicolor and virginica,
# leaving just the very well-separated setosa split from the rest.
rect.hclust(hc, h = 15, border = 1:2)

Reading a dendrogram

  • Height on the y-axis is the distance at which two clusters merge. Greater height means greater dissimilarity.
  • Large vertical gaps between merges suggest groups that are relatively well separated.
  • Where to cut the tree determines the number of clusters. Look for large jumps in height.

K-means Clustering

K-means partitions data into exactly k groups by minimising within-cluster variance. Unlike hierarchical clustering, which works from a precomputed distance matrix, standard k-means works directly in Euclidean feature space: it assigns each observation to the nearest cluster centroid, recalculates centroids, and repeats until convergence.

K-means is faster than hierarchical clustering and scales well to large datasets, but requires you to specify k in advance and works best when clusters are roughly compact and similarly shaped.

# Always set a seed for reproducibility
set.seed(123)

km <- kmeans(df_scaled, centers = 3, nstart = 25)
# nstart = 25 runs the algorithm from 25 random starts
# and returns the best solution. Use nstart > 1

# Visualise
fviz_cluster(km, data = df_scaled,
             palette = c("#E7B800", "#00AFBB", "#FC4E07"),
             ggtheme = theme_minimal())

Limitations

K-means works poorly when clusters are elongated or irregularly shaped, when cluster sizes differ greatly, or when outliers are present. For these cases, hierarchical clustering with complete or Ward linkage is more robust, or consider PAM (see Going Further).

Choosing the Number of Clusters

No single method reliably identifies the true number of clusters. Use at least two approaches and consider biological plausibility alongside statistical criteria.

Elbow plot shows within-cluster sum of squares against k. Look for a bend where adding more clusters gives diminishing returns:

fviz_nbclust(df_scaled, kmeans, method = "wss", k.max = 10)

Silhouette analysis compares average silhouette width across candidate values of k. Higher values indicate better separation; see Validation below for how silhouette width is calculated.

fviz_nbclust(df_scaled, kmeans, method = "silhouette", k.max = 10)

Gap statistic compares observed within-cluster dispersion to a null reference distribution. Choose the smallest k where the gap statistic is within one standard error of its maximum:

library(cluster)

gap_stat <- clusGap(df_scaled, FUN = kmeans, K.max = 10, B = 50)
fviz_gap_stat(gap_stat)

Practical advice

The three methods will not always agree. A biologically interpretable solution with k = 3 is often more useful than a statistically optimal solution with k = 7 that cannot be explained.

Validation

Clustering will always find groups, even in random data. This is especially true for methods that require fixing the number of clusters in advance, such as k-means or a forced cutree() cut: they will partition data into exactly k groups whether or not real structure exists. Density-based methods such as DBSCAN are a genuine exception, since they can legitimately report that no real clusters exist rather than being forced to produce some; see the DBSCAN section of Going Further for a worked example. Validation checks whether the clusters are meaningful, and it is worth distinguishing several separate questions:

  • Internal validation: are the clusters geometrically well separated? (silhouette, gap statistic)
  • Stability: does the clustering persist if the data, algorithm, or starting points change?
  • External validation: does the clustering correspond to independently known groups?
  • Biological interpretation: do the clusters make scientific sense?

Agreement on one of these does not guarantee the others, and the sections below work through each in turn.

Silhouette width

Silhouette width measures how well each sample fits its assigned cluster compared to neighbouring clusters. For each sample it asks two questions: how similar is this sample to the other samples in its own cluster, and how similar is it to samples in the nearest other cluster. The silhouette width is the difference between those two numbers, scaled to a range from -1 to 1.

library(cluster)

sil <- silhouette(km$cluster, dist_matrix)
fviz_silhouette(sil)
#   cluster size ave.sil.width
# 1       1   50          0.64
# 2       2   53          0.39
# 3       3   47          0.35
summary(sil)
plot(sil)

Interpreting the value:

Close to  1   sample fits well in its cluster
Close to  0   sample sits on the boundary between two clusters
Close to -1   sample probably belongs in the neighbouring cluster

As a rough guide, average silhouette widths above 0.5 suggest reasonably well-separated clusters, whereas values below 0.25 indicate weak structure. These are rules of thumb rather than universal cutoffs. On the iris example above, this is a good illustration of the point: cluster 1 (0.64) is clearly well separated, while clusters 2 and 3 (0.39 and 0.35) fall into the more ambiguous middle range rather than cleanly above or below a single threshold, consistent with two of the three species overlapping in feature space.

Silhouette width answers the question you always want to ask after clustering: did the algorithm actually find meaningful groups or did it just divide the data arbitrarily? It is the most widely used internal validation metric precisely because it does not require any external labels.

Stability across methods

The Adjusted Rand Index (mclust::adjustedRandIndex) measures the similarity between two cluster assignments. ARI = 1 indicates perfect agreement, ARI ≈ 0 indicates agreement no better than expected by chance, and negative values indicate less agreement than expected by chance.

library(mclust)

clusters_hc <- cutree(hc, k = 3)
clusters_km <- km$cluster

adjustedRandIndex(clusters_km, clusters_hc)
# [1] 0.63

If hierarchical and k-means clustering broadly agree, the result is less dependent on the choice of algorithm. This supports the stability of the solution, but does not by itself demonstrate biological reality: both algorithms are applied to the same data and can consistently pick up on the same geometric features, real or not.

Why not just use PERMANOVA?

It is tempting to test whether the clusters differ significantly using PERMANOVA (vegan::adonis2(dist_matrix ~ km$cluster)). This is circular, however, when the clusters being tested were generated from the same data. The clustering algorithm has already defined groups that maximise some aspect of separation, so testing those same groups on the same data does not provide independent evidence that the groups are meaningful.

PERMANOVA is genuinely useful for testing group differences, but only when the groups were defined independently of the multivariate data being tested, for example by treatment, habitat, or experimental condition. See the PERMANOVA page.

External validation

If true labels are available, compare cluster assignments directly against known groups:

table(iris$Species, km$cluster)
#              1  2  3
#   setosa    50  0  0
#   versicolor 0 39 11
#   virginica  0 14 36
adjustedRandIndex(iris$Species, km$cluster)
# [1] 0.62

The table() shows which species land in which cluster: setosa is perfectly separated, while versicolor and virginica overlap partially, with 11 versicolor and 14 virginica observations landing in the "wrong" cluster. The Adjusted Rand Index (0.62) reflects that partial overlap rather than perfect agreement.

External validation is only possible when true labels are known. In most real clustering applications no ground truth exists, which is why internal validation metrics such as silhouette width remain essential.

Common Pitfalls

Not scaling data. Variables with large ranges dominate Euclidean distance.

Using Euclidean distance for OTU counts. Use Bray-Curtis instead.

Setting nstart = 1 in k-means. The algorithm may converge to a poor local solution. Use multiple random starts, for example nstart = 25 or higher, to reduce this risk.

Over-interpreting clusters. Clustering is exploratory. Clusters found in the data are not necessarily biologically real. Validate them.

Exercise

Using the USArrests dataset:

  1. Scale the data and compute Euclidean distances
  2. Cluster with Ward's method and plot the dendrogram
  3. Cut the tree into 4 clusters and identify which states group together
  4. Confirm with k-means using the same k
  5. Compare the two solutions with the Adjusted Rand Index (mclust::adjustedRandIndex)
Solution 
library(factoextra)
library(mclust)

## 1. Scale and compute distances
df_scaled <- scale(USArrests)
dist_matrix <- dist(df_scaled)

## 2. Hierarchical clustering
hc <- hclust(dist_matrix, method = "ward.D2")
plot(hc, main = "US Arrests")

## 3. Cut and inspect
clusters_hc <- cutree(hc, k = 4)
rect.hclust(hc, k = 4, border = 2:5)
lapply(1:4, function(i) names(clusters_hc[clusters_hc == i]))

## 4. K-means
set.seed(123)
km <- kmeans(df_scaled, centers = 4, nstart = 25)

## 5. Compare
adjustedRandIndex(clusters_hc, km$cluster)
# This metric compares two cluster assignments: how much the two methods
# agree on which data points belong together, corrected for chance
# (random ~ 0, perfect = 1).
# A score of ~0.82 indicates strong agreement between the two methods.
# This suggests that despite using different algorithms, both Hierarchical
# Clustering and K-Means identified similar structure in the data, which is
# a good sign of stability, though not on its own proof of biological reality.