| SOCR ≫ | DSPA ≫ | DSPA3 Topics ≫ |
DSPA3: Data Science and Predictive Analytics (UMich HS650)
Chapter 8: Unsupervised Clustering
Chapter 8: Unsupervised Clustering
SOCR/MCAIM (Ivo Dinov)
SOCR/MCAIM (Ivo Dinov)
August 27, 2026
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(plotly) # interactive figures and ALL 3-D graphics
library(cluster) # silhouette, pam, agnes, daisy (Gower)
library(factoextra) # visualization and gap statistic
library(fpc) # cluster.stats, clusterboot
library(mclust) # Gaussian mixtures, adjustedRandIndexHow this chapter uses graphics
Every two-dimensional figure is drawn with
ggplot2and rendered statically. Immediately after each one, the equivalentplot_ly()code appears in a chunk markedeval=FALSE, echo=TRUE.Every three-dimensional figure is drawn with
plot_ly()and evaluated. Clustering is geometry, and several of its key objects only make sense in three dimensions: the \(k\)-means objective landscape over centroid coordinates, the spectral embedding in eigenvector space, and a mixture density over the feature plane.
1 Learning objectives
After completing this chapter you will be able to:
- Choose a distance appropriate to the variable types present, including Gower distance for mixed data, and explain what scaling a nominal code does to Euclidean geometry.
- Interpret an average silhouette against the standard thresholds, and recognize when a partition provides no evidence of structure.
- Measure agreement with known labels using the adjusted Rand index, and explain what the adjustment corrects for.
- State the \(k\)-means objective, explain why the problem is NP-hard, and prove that Lloyd’s algorithm converges.
- Explain why
nstartmatters and state the \(O(\log k)\) guarantee that \(k\)-means++ provides. - Select \(k\) using the elbow, silhouette, and gap statistic, and explain what null each compares against.
- Predict when \(k\)-means will fail from the shape of the data, and choose a method that will not.
- Apply DBSCAN, explain the role of \(\varepsilon\) and minPts, and read a \(k\)-distance plot.
- Fit a Gaussian mixture, select its covariance parameterization by BIC, and distinguish soft from hard assignment.
- Construct the three graph Laplacians correctly, explain why the trivial eigenvector is discarded, and analyze the complexity of spectral clustering.
Estimated time: 10–13 hours including exercises. Prerequisites: Chapter 3 (eigendecomposition, sparse matrices), Chapter 4 (distance concentration, §4.3.1), Chapter 5 (evaluation discipline), and Chapter 6 (the ensemble and kernel machinery this chapter reuses).
2 PART I: THE CLUSTERING PROBLEM
3 What clustering is, and is not
Supervised methods (Chapters 5 and 6) require labels. Clustering does not: it partitions cases into groups using only the features, and the groups are then interpreted after the fact.
That freedom comes with a cost that is easy to overlook.
Common misconception: “the clustering algorithm found three groups, so there are three groups.” Partitioning algorithms always return the number of groups you ask for. Run \(k\)-means with \(k=3\) on uniform noise and you get three clusters, with distinct centroids, of roughly equal size, and a plausible-looking bar chart of their centres. Nothing about the output signals that the input had no structure.
The output of a clustering algorithm is a hypothesis, not a finding. Whether groups exist is a separate question, answered by validation (§8.4), and if validation says there is no structure, the centroids should not be narrated.
set.seed(11)
noise <- data.frame(x1 = runif(300), x2 = runif(300)) # NO structure at all
km_noise <- kmeans(noise, centers = 3, nstart = 25)
p_noise <- ggplot(noise, aes(x1, x2, color = factor(km_noise$cluster))) +
geom_point(size = 1.6, alpha = 0.8) +
scale_color_brewer(palette = "Set1", name = "Cluster") +
coord_fixed() +
labs(title = "k-means on uniform noise",
subtitle = "Three clean-looking clusters, in data that has none",
x = expression(x[1]), y = expression(x[2])) + theme_dspa(10)
cent <- as.data.frame(km_noise$centers) |>
mutate(cluster = factor(1:3)) |>
pivot_longer(-cluster, names_to = "feature", values_to = "centre")
p_cent <- ggplot(cent, aes(feature, centre, fill = cluster)) +
geom_col(position = "dodge") +
scale_fill_brewer(palette = "Set1") +
labs(title = "Their centroids look distinct",
subtitle = "Which is guaranteed: k-means separates centroids by construction",
x = NULL, y = "Centroid value") + theme_dspa(10)
p_noise | p_cent
sil_noise <- silhouette(km_noise$cluster, dist(noise))
c(average_silhouette = round(mean(sil_noise[, "sil_width"]), 4),
interpretation = "below 0.26 -- no substantial structure")#> average_silhouette
#> "0.3941"
#> interpretation
#> "below 0.26 -- no substantial structure"The silhouette is the diagnostic that the centroid plot cannot provide.
3.1 A first example
The SOCR hotdogs data records calories and sodium for 54 frankfurters, along with the primary meat type, which we will hold back and use only to check the result.
library(rvest)
hotdog <- read_html("https://wiki.socr.umich.edu/index.php/SOCR_012708_ID_Data_HotDogs") |>
html_nodes("table") |> _[[1]] |> html_table()
c(rows = nrow(hotdog), types = paste(unique(hotdog$Type), collapse = ", "))#> rows types
#> "54" "Beef, Meat, Poultry"# Clusters COMPUTED, not drawn by hand
hd <- scale(hotdog[, c("Calories", "Sodium")])
set.seed(7)
km_hd <- kmeans(hd, centers = 3, nstart = 25)
hotdog$cluster <- factor(km_hd$cluster)
p_hd1 <- ggplot(hotdog, aes(Calories, Sodium, color = cluster)) +
geom_point(size = 2.4) +
stat_ellipse(level = 0.68, linewidth = 0.5) +
scale_color_brewer(palette = "Set1") +
labs(title = "k-means partition (k = 3)", x = "Calories", y = "Sodium") +
theme_dspa(10)
p_hd2 <- ggplot(hotdog, aes(Calories, Sodium, color = Type, shape = Type)) +
geom_point(size = 2.4) +
scale_color_brewer(palette = "Dark2") +
labs(title = "Actual meat type (held back)", x = "Calories", y = "Sodium") +
theme_dspa(10)
p_hd1 | p_hd2
#> type
#> cluster Beef Meat Poultry
#> 1 0 0 9
#> 2 12 9 3
#> 3 8 8 5c(adjusted_rand_index = round(mclust::adjustedRandIndex(km_hd$cluster, hotdog$Type), 4),
average_silhouette = round(mean(silhouette(km_hd$cluster, dist(hd))[, "sil_width"]), 4))#> adjusted_rand_index average_silhouette
#> 0.1189 0.4699The clusters do capture something, poultry separates cleanly on calories, but beef and meat overlap almost completely, and the adjusted Rand index quantifies how partial the agreement is. That number is the honest summary. A side-by-side plot invites the eye to see more agreement than there is.
# --- Interactive equivalent ------------------------------------------------
plot_ly(hotdog, x = ~Calories, y = ~Sodium, color = ~Type, symbol = ~Type,
type = "scatter", mode = "markers", marker = list(size = 12)) |>
layout(title = "Hotdogs: calories vs. sodium, by meat type",
xaxis = list(title = "Calories"), yaxis = list(title = "Sodium"),
legend = list(title = list(text = "<b>Meat type</b>"),
orientation = "h"))4 Distances
Every clustering method in this chapter reduces the data to pairwise dissimilarities. The choice of dissimilarity is therefore the most consequential modelling decision, and it is made before any algorithm runs.
4.1 The distance axioms
A function \(d\) is a metric if for all \(x,y,z\):
\[ \begin{aligned} &\text{(i) } d(x,y)\ge0,\quad d(x,y)=0\iff x=y &&\text{(non-negativity, identity)}\\ &\text{(ii) } d(x,y)=d(y,x) &&\text{(symmetry)}\\ &\text{(iii) } d(x,z)\le d(x,y)+d(y,z) &&\text{(triangle inequality)} \end{aligned} \]
Not every useful dissimilarity is a metric, cosine distance fails (iii), as Chapter 7 showed, but methods that rely on metric structure (metric trees, some approximation guarantees) require it.
| Distance | Definition | Assumes |
|---|---|---|
| Euclidean | \(\sqrt{\sum_j(x_j-y_j)^2}\) | Interval scale, comparable units, isotropy |
| Manhattan | \(\sum_j\lvert x_j-y_j\rvert\) | Interval scale; more robust to outliers |
| Minkowski (\(L_p\)) | \(\big(\sum_j\lvert x_j-y_j\rvert^p\big)^{1/p}\) | \(p=2\) Euclidean, \(p=1\) Manhattan |
| Mahalanobis | \(\sqrt{(x-y)^\top S^{-1}(x-y)}\) | Interval scale; corrects for correlation and scale |
| Cosine | \(1-\dfrac{x^\top y}{\lVert x\rVert\lVert y\rVert}\) | Direction matters, magnitude does not |
| Gower | see §8.3 | Mixed types; handles missingness |
Mahalanobis distance deserves attention. It equals Euclidean distance after whitening by \(S^{-1/2}\), so it automatically corrects for both differing variances and correlation between features, which is what standardization does only partially. Its cost is estimating and inverting \(S\), which needs \(n\gg d\) and is unstable when features are collinear (Chapter 3, §3.7).
set.seed(21)
Sg <- matrix(c(1, 0.9, 0.9, 1), 2)
X <- MASS::mvrnorm(300, c(0, 0), Sg)
probe <- c(2, -1) # a point off the correlation axis
d_euc <- sqrt(colSums((t(X) - probe)^2))
d_mah <- sqrt(mahalanobis(X, probe, cov(X)))
dd <- data.frame(x1 = X[, 1], x2 = X[, 2],
Euclidean = d_euc, Mahalanobis = d_mah) |>
pivot_longer(c(Euclidean, Mahalanobis), names_to = "metric", values_to = "d")
ggplot(dd, aes(x1, x2, color = d)) +
geom_point(size = 1.3) +
geom_point(data = data.frame(x1 = probe[1], x2 = probe[2]),
aes(x1, x2), inherit.aes = FALSE, color = "red", size = 4,
shape = 4, stroke = 1.5) +
scale_color_viridis_c(name = "Distance") +
facet_wrap(~ metric) + coord_fixed() +
labs(title = "Euclidean and Mahalanobis distance from the same probe point",
subtitle = "Euclidean contours are circles; Mahalanobis contours follow the data's correlation",
x = expression(x[1]), y = expression(x[2])) +
theme_dspa(10)
4.1.1 Adjusted Rand Index (ARI) for contrasting a pair of clustering methods
To compare two data clusterings or classifications, \(U = \{u_1, u_2, \dots, u_R\}\) and \(V = \{v_1, v_2, \dots, v_C\}\), over a dataset of \(n\) elements, the Adjusted Rand Index (ARI) measures the similarity between the two assignments while adjusting for chance grouping.
Let \(n_{ij}\) denote the number of objects that are both in cluster \(u_i\) of partition \(U\) and in cluster \(v_j\) of partition \(V\). This information is typically organized in an \(R \times C\) contingency table.
| \(U \backslash V\) | \(v_1\) | \(v_2\) | \(\dots\) | \(v_C\) | Row Sums |
|---|---|---|---|---|---|
| \(u_1\) | \(n_{11}\) | \(n_{12}\) | \(\dots\) | \(n_{1C}\) | \(a_1\) |
| \(u_2\) | \(n_{21}\) | \(n_{22}\) | \(\dots\) | \(n_{2C}\) | \(a_2\) |
| \(\vdots\) | \(\vdots\) | \(\vdots\) | \(\ddots\) | \(\vdots\) | \(\vdots\) |
| \(u_R\) | \(n_{R1}\) | \(n_{R2}\) | \(\dots\) | \(n_{RC}\) | \(a_R\) |
| Col Sums | \(b_1\) | \(b_2\) | \(\dots\) | \(b_C\) | \(n\) |
where \(a_i = \sum_{j=1}^{C} n_{ij}\) represents the total number of elements in cluster \(u_i\), \(b_j = \sum_{i=1}^{R} n_{ij}\) represents the total number of elements in cluster \(v_j\), and \(n = \sum_{i=1}^{R} \sum_{j=1}^{C} n_{ij}\) is the total number of elements.
The standard (unadjusted) Rand Index considers all pairs of elements.
Agreement: Put in the same cluster in \(U\) and the same cluster in \(V\), or different clusters in \(U\) and different clusters in \(V\).
The number of such agreement pairs is calculated using combinatorial
choose function is
\[\text{Index} = \sum_{i=1}^{R} \sum_{j=1}^{C} \binom{n_{ij}}{2}\]
\[\text{Expected pairs in same group from row sums} = \sum_{i=1}^{R} \binom{a_i}{2}\]
\[\text{Expected pairs in same group from col sums} = \sum_{j=1}^{C} \binom{b_j}{2}\]
The unadjusted Rand Index is
\[\text{RI} = \frac{\sum_{ij} \binom{n_{ij}}{2} + \frac{1}{2} \left[ \sum_i \binom{a_i}{2} + \Sigma_j \binom{b_j}{2} \right] - \frac{1}{\binom{n}{2}} \left( \sum_i \binom{a_i}{2} \sum_j \binom{b_j}{2} \right)}{\binom{n}{2}}\]
Note that RI approaches \(1\) when the partitions are identical, but its expected value for random clusterings is not constant-zero across different marginal distributions,
The **Adjusted Rand Index (ARI)* corrects for random chance by normalizing the index by subtracting the expected similarity under a hypergeometric random model
\[\text{ARI} = \frac{\sum_{i=1}^{R} \sum_{j=1}^{C} \binom{n_{ij}}{2} - \left[ \sum_{i=1}^{R} \binom{a_i}{2} \sum_{j=1}^{C} \binom{b_j}{2} \right] / \binom{n}{2}}{\frac{1}{2} \left[ \sum_{i=1}^{R} \binom{a_i}{2} + \sum_{j=1}^{C} \binom{b_j}{2} \right] - \left[ \sum_{i=1}^{R} \binom{a_i}{2} \sum_{j=1}^{C} \binom{b_j}{2} \right] / \binom{n}{2}}.\]
The key ARI properties include
- Range: \(\text{ARI} \in [-1, 1]\).
- \(1\) indicates identical clusterings (up to permutation).
- \(0\) indicates agreement equal to what would be expected by random chance alone.
- Negative values (rarely) indicate worse-than-random agreement.
- Invariance: ARI is invariant to permutations of cluster labels.
4.2 Distance concentration
Chapter 4, §4.3.1 established that pairwise distances concentrate as dimension grows: the ratio of their spread to their mean tends to zero. Clustering is a direct casualty, because every method here assumes “near” is meaningfully different from “far.”
conc_cluster <- function(d, n = 400, seed = 31) {
set.seed(seed + d)
# Two well-separated groups in the FIRST two dimensions; the rest is noise
y <- rep(1:2, each = n / 2)
Z <- matrix(rnorm(n * d), n, d)
Z[y == 2, 1:2] <- Z[y == 2, 1:2] + 3
dd <- as.numeric(dist(Z))
km <- kmeans(Z, 2, nstart = 25)
c(d = d, relative_spread = sd(dd) / mean(dd),
ARI = mclust::adjustedRandIndex(km$cluster, y),
silhouette = mean(silhouette(km$cluster, dist(Z))[, "sil_width"]))
}
cc <- as.data.frame(do.call(rbind, lapply(c(2, 5, 20, 50, 200, 500), conc_cluster)))
round(cc, 4)cc |> select(d, ARI, silhouette) |>
pivot_longer(-d, names_to = "measure", values_to = "value") |>
ggplot(aes(d, value, color = measure)) +
geom_line(linewidth = 1) + geom_point(size = 2.2) +
scale_x_log10() +
scale_color_manual(values = c(ARI = "#3B7DD8", silhouette = "#D8433B")) +
labs(title = "Adding noise dimensions destroys recoverable structure",
subtitle = "The true separation never changes -- it always lives in the first two dimensions",
x = "Ambient dimension d (log scale)", y = NULL, color = NULL) +
theme_dspa()
The groups are separated by three standard deviations in dimensions 1–2 throughout. By \(d=500\) the recovery has collapsed. Reduce dimension before clustering, PCA, UMAP (Chapter 4), or feature selection, rather than hoping the algorithm will find the informative directions.
5 Preparing mixed-type data
Real datasets rarely contain only interval-scaled variables. A single frame may hold binary indicators, ordinal scales, unordered categories, and sentinel codes for missingness, and Euclidean distance is defined for none of the last three.
Common misconception: “standardize everything, then use Euclidean distance.” Standardization rescales; it does not change a variable’s type. Consider a nominal code with levels 1 = mother only, 2 = father only, 3 = both parents. Scaling it and computing Euclidean distance asserts that “father only” lies exactly halfway between the other two, and that the gap from 1 to 3 is twice the gap from 1 to 2. Neither statement means anything: the labels could equally have been assigned in any order.
The failure is worse for sentinel codes. If a variable uses 1, 2, 3 for real responses and 8 for “don’t know,” then after standardization those respondents sit roughly two standard deviations from everyone else, where they dominate every distance and pull a centroid toward themselves. A missing value has become the most influential observation in the dataset.
set.seed(41)
n_s <- 200
v_real <- sample(1:3, n_s, TRUE, prob = c(0.4, 0.4, 0.2))
v_real[sample(n_s, 12)] <- 8 # 6% coded "don't know" as 8
df_s <- data.frame(raw = v_real, scaled = as.numeric(scale(v_real)),
status = ifelse(v_real == 8, "sentinel (8)", "real response"))
ggplot(df_s, aes(scaled, fill = status)) +
geom_histogram(bins = 30, color = "white") +
scale_fill_manual(values = c("real response" = "#9EB4C8",
"sentinel (8)" = "#D8433B")) +
labs(title = "What standardizing a sentinel code does",
subtitle = sprintf("6%% of cases sit %.1f standard deviations from the rest, purely because 8 codes 'don't know'",
max(df_s$scaled) - median(df_s$scaled)),
x = "Standardized value", y = "Count", fill = NULL) +
theme_dspa()
5.1 Gower distance
Gower’s coefficient (1971) handles mixed types by computing a per-variable dissimilarity on a common \([0,1]\) scale and averaging:
\[\boxed{\;d_{\mathrm{Gower}}(i,j)=\frac{\sum_{p=1}^{P}w_{ijp}\,d_{ijp}}{\sum_{p=1}^{P}w_{ijp}}\;}\]
with the per-variable term chosen by type:
\[ d_{ijp}=\begin{cases} \dfrac{\lvert x_{ip}-x_{jp}\rvert}{R_p} & \text{quantitative, }R_p=\text{range}\\[2mm] \mathbb{1}\{x_{ip}\ne x_{jp}\} & \text{nominal (simple matching)}\\[2mm] \dfrac{\lvert r_{ip}-r_{jp}\rvert}{\max r_p-1} & \text{ordinal (on ranks)} \end{cases} \]
and \(w_{ijp}=0\) when either value is missing, so missingness is handled by down-weighting rather than imputation, which is exactly what a sentinel code should trigger.
mixed <- data.frame(
age = c(23, 45, 31, 52, 29), # quantitative
score = c(2.1, 3.8, 2.9, 4.4, NA), # quantitative, missing
region = factor(c("north", "south", "north", "east", "south")), # nominal
grade = factor(c("low", "high", "mid", "high", "low"), # ordinal
levels = c("low", "mid", "high"), ordered = TRUE))
str(mixed)#> 'data.frame': 5 obs. of 4 variables:
#> $ age : num 23 45 31 52 29
#> $ score : num 2.1 3.8 2.9 4.4 NA
#> $ region: Factor w/ 3 levels "east","north",..: 2 3 2 1 3
#> $ grade : Ord.factor w/ 3 levels "low"<"mid"<"high": 1 3 2 3 1# daisy() dispatches on column type; ordered factors are handled as ranks
d_gow <- daisy(mixed, metric = "gower")
round(as.matrix(d_gow), 3)#> 1 2 3 4 5
#> 1 0.000 0.874 0.281 1.000 0.402
#> 2 0.874 0.000 0.594 0.376 0.517
#> 3 0.281 0.594 0.000 0.719 0.523
#> 4 1.000 0.376 0.719 0.000 0.931
#> 5 0.402 0.517 0.523 0.931 0.000Compare against what a naive numeric coercion produces:
naive <- data.frame(lapply(mixed, \(x) as.numeric(x)))
naive$score[is.na(naive$score)] <- mean(naive$score, na.rm = TRUE)
d_euc_naive <- dist(scale(naive))
c(gower_case1_vs_case5 = round(as.matrix(d_gow)[1, 5], 3),
gower_case1_vs_case3 = round(as.matrix(d_gow)[1, 3], 3),
euclid_case1_vs_case5 = round(as.matrix(d_euc_naive)[1, 5], 3),
euclid_case1_vs_case3 = round(as.matrix(d_euc_naive)[1, 3], 3),
note = "cases 1 and 3 share region; cases 1 and 5 do not")#> gower_case1_vs_case5
#> "0.402"
#> gower_case1_vs_case3
#> "0.281"
#> euclid_case1_vs_case5
#> "1.887"
#> euclid_case1_vs_case3
#> "1.509"
#> note
#> "cases 1 and 3 share region; cases 1 and 5 do not"Gower registers the shared region between cases 1 and 3;
the numeric coercion treats north=2 and
south=3 as being one unit apart and east=1 as
two units from south, an ordering with no meaning.
Gower distance pairs naturally with \(k\)-medoids (§8.9), because PAM works from a dissimilarity matrix and needs no notion of a centroid mean — which is undefined for a nominal variable anyway.
6 Validating a clustering
Three families of question, each answered by different tools.
6.1 Internal indices
Internal measures use only the data and the partition.
Silhouette. For observation \(i\) in cluster \(C_i\), let \(a_i\) be its mean distance to the other members of \(C_i\), and \(b_i\) the smallest mean distance to any other cluster. Then
\[\boxed{\;s_i=\frac{b_i-a_i}{\max\{a_i,b_i\}}\in[-1,1]\;}\]
Reading the three regimes: \(s_i\to1\) when \(a_i\ll b_i\) (well inside its cluster); \(s_i\approx0\) on a boundary; \(s_i\to-1\) when \(a_i\gg b_i\) (closer to a neighbouring cluster than its own, likely misassigned).
The average silhouette has standard thresholds, and they are stricter than most readers expect. From Kaufman and Rousseeuw’s Finding Groups in Data (1990):
Average \(s\) Interpretation \(0.71\)–\(1.00\) Strong structure \(0.51\)–\(0.70\) Reasonable structure \(0.26\)–\(0.50\) Weak structure; could be artificial \(<0.26\) No substantial structure has been found An average silhouette of 0.20 is not a modest success. It is the diagnostic reporting that the data does not support the partition, and when it says that, the centroids should not be given a narrative.
sil_verdict <- function(s) {
cut(s, breaks = c(-Inf, 0.26, 0.51, 0.71, Inf),
labels = c("no substantial structure", "weak / possibly artificial",
"reasonable structure", "strong structure"),
right = FALSE)
}
data.frame(average_silhouette = c(0.12, 0.20, 0.35, 0.58, 0.79),
verdict = as.character(sil_verdict(c(0.12, 0.20, 0.35, 0.58, 0.79))))Two further internal indices behave differently and are worth reporting alongside:
\[\textbf{Calinski–Harabasz: }\; \mathrm{CH}(k)=\frac{\mathrm{tr}(B_k)/(k-1)}{\mathrm{tr}(W_k)/(n-k)}, \qquad\text{larger is better}\]
\[\textbf{Davies–Bouldin: }\; \mathrm{DB}(k)=\frac1k\sum_{i=1}^{k}\max_{j\ne i}\frac{\sigma_i+\sigma_j}{d(\mu_i,\mu_j)}, \qquad\text{smaller is better}\]
CH is a variance-ratio \(F\)-like statistic; DB averages the worst-case similarity of each cluster to any other. They can disagree with silhouette, and disagreement is informative.
set.seed(51)
# Three genuinely separated groups
g <- rep(1:3, each = 120)
Xg <- MASS::mvrnorm(360, c(0, 0), diag(2))
Xg[g == 2, ] <- Xg[g == 2, ] + c(5, 0)
Xg[g == 3, ] <- Xg[g == 3, ] + c(2.5, 4.5)
idx_tab <- do.call(rbind, lapply(2:8, function(k) {
km <- kmeans(Xg, k, nstart = 25)
st <- fpc::cluster.stats(dist(Xg), km$cluster)
data.frame(k = k, silhouette = st$avg.silwidth,
CH = st$ch, DB = fpc::cluster.stats(dist(Xg), km$cluster)$sindex,
within_ss = km$tot.withinss)
}))
idx_tab |> mutate(across(where(is.numeric), \(z) round(z, 3)))idx_tab |> select(k, silhouette, CH) |>
pivot_longer(-k, names_to = "index", values_to = "value") |>
ggplot(aes(k, value)) +
geom_line(linewidth = 1, color = "steelblue") + geom_point(size = 2.2) +
facet_wrap(~ index, scales = "free_y") +
scale_x_continuous(breaks = 2:8) +
labs(title = "Two internal indices on data with three true groups",
subtitle = "Both peak at k = 3; when they disagree, that disagreement is the finding",
x = "k", y = NULL) +
theme_dspa(10)
6.2 External measures
When true labels exist, for validation, not for fitting, agreement must account for the fact that cluster labels are arbitrary permutations and that some agreement arises by chance.
The Rand index counts pairs of observations classified concordantly:
\[\mathrm{RI}=\frac{a+d}{\binom{n}{2}},\]
where \(a\) is the number of pairs together in both partitions and \(d\) the number apart in both. Its problem is that RI does not approach 0 for a random partition, with many clusters it drifts upward regardless.
The adjusted Rand index corrects for chance under a generalized hypergeometric null:
\[\boxed{\;\mathrm{ARI}=\frac{\mathrm{RI}-\mathbb{E}[\mathrm{RI}]}{\max(\mathrm{RI})-\mathbb{E}[\mathrm{RI}]}\;}\]
so \(\mathrm{ARI}=0\) for a random partition and \(1\) for perfect agreement. It can be negative, meaning worse-than-chance agreement.
set.seed(61)
truth <- rep(1:4, each = 60)
random_partition <- sample(1:4, 240, TRUE)
km_truth <- kmeans(cbind(truth + rnorm(240, sd = 0.3), rnorm(240)), 4, nstart = 25)
data.frame(
comparison = c("Random partition vs. truth", "k-means vs. truth",
"Truth vs. itself"),
Rand = round(c(fpc::cluster.stats(dist(cbind(truth, 0)), random_partition,
truth)$corrected.rand,
NA, NA), 4),
ARI = round(c(mclust::adjustedRandIndex(random_partition, truth),
mclust::adjustedRandIndex(km_truth$cluster, truth),
mclust::adjustedRandIndex(truth, truth)), 4))A random partition scores near zero on ARI, which is the property that makes it interpretable. Normalized mutual information is the information-theoretic alternative, \(\mathrm{NMI}=I(U;V)/\sqrt{H(U)H(V)}\in[0,1]\), and behaves similarly.
Common misconception: “I validated the clustering by cross-tabulating it against one of the variables.” If that variable was among the features the algorithm was fitted on, the comparison is circular, the model is being asked to recover information it was handed, and agreement is guaranteed by construction. It is the unsupervised form of the target leakage described in Chapter 5, §5.4.1, and it is easy to commit because unsupervised methods have no obvious “outcome” to protect.
A validation variable must be held out of the fit entirely. If none is available, use internal indices and stability instead, and say so, a stability result honestly reported is worth more than a confusion matrix that could not have come out any other way.
set.seed(65)
# A feature GIVEN to the algorithm vs. one HELD OUT, on the same data
Xv <- cbind(rnorm(300), rnorm(300), rnorm(300)) # no structure at all
given <- Xv[, 3] # used in the fit
heldout <- rnorm(300) # never seen by the model
km_v <- kmeans(Xv, 3, nstart = 25)
c(ARI_vs_binned_INPUT_feature =
round(mclust::adjustedRandIndex(km_v$cluster, cut(given, 3)), 4),
ARI_vs_binned_HELD_OUT_feature =
round(mclust::adjustedRandIndex(km_v$cluster, cut(heldout, 3)), 4))#> ARI_vs_binned_INPUT_feature ARI_vs_binned_HELD_OUT_feature
#> 0.1100 -0.0014There is no structure in this data, yet the clustering “agrees” with a feature it was given and shows nothing against one it was not. Only the second number is evidence.
6.3 Stability
A partition that changes completely under resampling is not a finding about the data. Bootstrap stability resamples, re-clusters, and measures how often each cluster is recovered, using the Jaccard coefficient between the original cluster and its best match.
set.seed(71)
stab_real <- fpc::clusterboot(Xg, B = 40, clustermethod = fpc::kmeansCBI,
krange = 3, seed = 71, count = FALSE)
stab_noise <- fpc::clusterboot(as.matrix(noise), B = 40,
clustermethod = fpc::kmeansCBI,
krange = 3, seed = 71, count = FALSE)
data.frame(
data = c("Three real groups", "Uniform noise"),
mean_jaccard = round(c(mean(stab_real$bootmean), mean(stab_noise$bootmean)), 4),
min_jaccard = round(c(min(stab_real$bootmean), min(stab_noise$bootmean)), 4))The conventional reading: mean Jaccard above 0.85 indicates a highly stable, reproducible cluster; 0.60–0.75 is patchy; below 0.60 means the cluster is not reproducible and should not be interpreted. The contrast between the two rows above is the whole point, a stability check separates real groups from algorithmic artifacts even when the centroids look equally convincing.
7 PART II: PARTITIONING METHODS
8 k-means
8.1 The objective
Given \(\{x_1,\dots,x_n\}\subset\mathbb{R}^d\) and a target \(k\le n\), partition the observations into \(S=\{S_1,\dots,S_k\}\) minimizing the within-cluster sum of squares:
\[\boxed{\;\arg\min_{S}\ \sum_{i=1}^{k}\sum_{x\in S_i}\lVert x-\mu_i\rVert^2, \qquad \mu_i=\frac{1}{|S_i|}\sum_{x\in S_i}x\;}\]
Three equivalent readings, each illuminating:
\[\sum_{i=1}^{k}\sum_{x\in S_i}\lVert x-\mu_i\rVert^2 =\sum_{i=1}^{k}|S_i|\operatorname{Var}(S_i) =\sum_{i=1}^{k}\frac{1}{2|S_i|}\sum_{x,y\in S_i}\lVert x-y\rVert^2 .\]
The middle form shows \(k\)-means minimizes a size-weighted average within-cluster variance. The right form uses the identity \(|S_i|\sum_{x\in S_i}\lVert x-\mu_i\rVert^2=\tfrac12\sum_{x,y\in S_i}\lVert x-y\rVert^2\), so the objective depends only on pairwise distances within clusters, no centroid required to state it.
Because the total sum of squares is fixed, minimizing within-cluster scatter is equivalent to maximizing between-cluster scatter:
\[\underbrace{\mathrm{TSS}}_{\text{fixed}}=\underbrace{\mathrm{WSS}}_{\text{minimize}}+\underbrace{\mathrm{BSS}}_{\text{maximize}} .\]
That is the same decomposition as the ANOVA identity of Chapter 3, §3.10.2, applied to a partition rather than a design matrix.
8.2 The problem is NP-hard
\(k\)-means is NP-hard, even for \(k=2\) in general dimension (Aloise et al., 2009) and even for \(d=2\) with general \(k\) (Mahajan et al., 2009). The number of partitions of \(n\) points into \(k\) non-empty groups is the Stirling number of the second kind, \(S(n,k)\sim k^n/k!\), exhaustive search is hopeless.
stirling2 <- function(n, k) sum((-1)^(0:k) * choose(k, 0:k) * (k - 0:k)^n) / factorial(k)
data.frame(n = c(10, 20, 30, 47),
partitions_into_3 = sapply(c(10, 20, 30, 47), stirling2, k = 3)) |>
mutate(partitions_into_3 = format(partitions_into_3, scientific = TRUE, digits = 3))This is why Lloyd’s algorithm is a heuristic with no approximation guarantee, and why everything in §8.6 about initialization matters.
8.3 Lloyd’s algorithm
- Initialize. Choose \(k\) starting centres \(\mu_1^{(0)},\dots,\mu_k^{(0)}\).
- Assign. Each point goes to its nearest centre: \[S_i^{(t)}=\Big\{x_p:\lVert x_p-\mu_i^{(t)}\rVert^2\le\lVert x_p-\mu_j^{(t)}\rVert^2\ \ \forall j\Big\}.\] This is a Voronoi tessellation of the space by the current centres.
- Update. Recompute each centre as the mean of its cluster: \[\mu_i^{(t+1)}=\frac{1}{|S_i^{(t)}|}\sum_{x\in S_i^{(t)}}x .\]
- Repeat until assignments stop changing.
Convergence. Each step cannot increase the objective, so the algorithm terminates in finitely many iterations.
Proof. The assignment step moves each point to the centre minimizing its squared distance, which can only lower \(\sum_i\sum_{x\in S_i}\lVert x-\mu_i\rVert^2\). The update step replaces each \(\mu_i\) by the minimizer of \(\sum_{x\in S_i}\lVert x-\mu\rVert^2\) over \(\mu\), the mean, since \(\nabla_\mu\sum_{x\in S_i}\lVert x-\mu\rVert^2=-2\sum_{x\in S_i}(x-\mu)=0\) gives \(\mu=\bar x_{S_i}\), which again cannot increase the objective. The objective is bounded below by 0 and there are finitely many partitions, so the sequence must stabilize. \(\blacksquare\)
Convergence to a local minimum is guaranteed; convergence to the global minimum is not.
The SOCR Voronoi tessellation app demonstrates the assignment step interactively.
lloyd <- function(X, centres, max_iter = 100) {
X <- as.matrix(X); k <- nrow(centres)
hist <- numeric(0)
for (t in seq_len(max_iter)) {
D <- as.matrix(dist(rbind(centres, X)))[-(1:k), 1:k, drop = FALSE]
cl <- max.col(-D, ties.method = "first") # assignment step
new_c <- t(vapply(seq_len(k), function(i) {
if (!any(cl == i)) centres[i, ] else colMeans(X[cl == i, , drop = FALSE])
}, numeric(ncol(X)))) # update step
wss <- sum(vapply(seq_len(k), function(i)
if (any(cl == i)) sum((sweep(X[cl == i, , drop = FALSE], 2, new_c[i, ]))^2) else 0,
numeric(1)))
hist <- c(hist, wss)
if (isTRUE(all.equal(centres, new_c))) { centres <- new_c; break }
centres <- new_c
}
list(cluster = cl, centers = centres, wss_trace = hist, iterations = t)
}
set.seed(81)
fit_l <- lloyd(Xg, Xg[sample(nrow(Xg), 3), , drop = FALSE])
c(iterations = fit_l$iterations,
monotone_decreasing = all(diff(fit_l$wss_trace) <= 1e-8),
final_wss = round(min(fit_l$wss_trace), 3))#> iterations monotone_decreasing final_wss
#> 10.00 1.00 663.88#> [1] 943.86 863.06 789.91 750.45 721.78 684.81 665.88 664.29 663.88 663.88The trace is monotone, exactly as the proof requires.
8.4 Which implementation?
kmeans() offers three algorithms, and they find
different local optima.
| Algorithm | Move considered | Note |
|---|---|---|
| Lloyd / Forgy | Batch: reassign all, then recompute | Simplest; most prone to poor local optima |
| MacQueen | Online: update centres after each point moves | Order-dependent |
| Hartigan–Wong | Moves a point if doing so lowers the objective, even when the point is already nearest its own centre | R’s default; generally finds better optima |
Hartigan–Wong’s extra power comes from a subtlety: after moving \(x\) from \(S_i\) to \(S_j\), both centroids shift, so a move can reduce the total objective even when \(x\) is currently closest to \(\mu_i\). Lloyd never considers such moves.
set.seed(91)
compare_alg <- function(alg, reps = 40) {
v <- vapply(seq_len(reps), function(r) {
set.seed(100 + r)
kmeans(Xg, 3, nstart = 1, algorithm = alg)$tot.withinss
}, numeric(1))
c(algorithm = alg, best = round(min(v), 2), median = round(median(v), 2),
worst = round(max(v), 2))
}
as.data.frame(do.call(rbind, lapply(c("Hartigan-Wong", "Lloyd", "MacQueen"),
compare_alg)))Across 40 single-start runs, Hartigan–Wong’s median is closest to the best solution any of them found. Use one algorithm consistently; switching between them mid-analysis makes reported objectives incomparable.
9 Initialization
9.1 Why
nstart is the most important argument
The objective is non-convex, so the basin Lloyd’s algorithm falls into is determined entirely by the starting centres. A single run is one draw from a multimodal distribution.
set.seed(101)
single_starts <- vapply(1:200, function(r) {
set.seed(200 + r)
kmeans(Xg, 3, nstart = 1)$tot.withinss
}, numeric(1))
best25 <- kmeans(Xg, 3, nstart = 25)$tot.withinss
c(single_start_best = round(min(single_starts), 2),
single_start_median = round(median(single_starts), 2),
single_start_worst = round(max(single_starts), 2),
nstart_25 = round(best25, 2),
pct_of_single_starts_that_are_optimal =
round(100 * mean(abs(single_starts - min(single_starts)) < 1e-6), 1))#> single_start_best single_start_median
#> 663.88 664.15
#> single_start_worst nstart_25
#> 664.15 663.88
#> pct_of_single_starts_that_are_optimal
#> 41.00ggplot(data.frame(wss = single_starts), aes(wss)) +
geom_histogram(bins = 40, fill = "grey80", color = "white") +
geom_vline(xintercept = best25, color = "firebrick", linewidth = 1.1) +
labs(title = "200 single-start k-means runs on the same data",
subtitle = "Red line: the solution found by nstart = 25. The objective is multimodal",
x = "Total within-cluster sum of squares", y = "Runs") +
theme_dspa()
Common misconception: “
kmeans(x, k)gives the \(k\)-means solution.” It gives a local minimum, reached from one random start. The default isnstart = 1. Setnstart = 25(or more) so the function runs the algorithm repeatedly and returns the best result it found. This costs a factor ofnstartin time and is essentially always worth it, a poor local optimum is not a modelling choice, it is a bug you cannot see.
The objective over centroid positions is a landscape, and restarts are paths across it. This is worth rotating:
# Fix two centres at their optimal positions and vary the third across the plane
km_opt <- kmeans(Xg, 3, nstart = 50)
fixed <- km_opt$centers[1:2, , drop = FALSE]
gx <- seq(min(Xg[, 1]) - 1, max(Xg[, 1]) + 1, length.out = 60)
gy <- seq(min(Xg[, 2]) - 1, max(Xg[, 2]) + 1, length.out = 60)
wss_at <- function(cx, cy) {
C <- rbind(fixed, c(cx, cy))
D <- as.matrix(dist(rbind(C, Xg)))[-(1:3), 1:3, drop = FALSE]
sum(apply(D, 1, min)^2)
}
Z <- outer(gy, gx, Vectorize(wss_at))
plot_ly(x = gx, y = gy, z = Z, type = "surface", colorscale = "Viridis",
opacity = 0.9, colorbar = list(title = "WSS")) |>
add_trace(x = km_opt$centers[3, 1], y = km_opt$centers[3, 2],
z = wss_at(km_opt$centers[3, 1], km_opt$centers[3, 2]) + 30,
type = "scatter3d", mode = "markers", name = "Optimal third centre",
marker = list(size = 6, color = "red")) |>
layout(title = "The k-means objective as the third centroid moves (other two fixed at their optima)",
scene = list(xaxis = list(title = "centroid x"),
yaxis = list(title = "centroid y"),
zaxis = list(title = "Within-cluster SS")))Rotate to see the multiple basins. A single random
start lands in whichever one it happens to fall into;
nstart samples several and keeps the deepest.
9.2 k-means++
Random initialization can place two centres in the same true cluster, leaving another unrepresented. \(k\)-means++ (Arthur & Vassilvitskii, 2007) seeds centres far apart with high probability:
- Choose \(\mu_1\) uniformly at random from the data.
- For \(j=2,\dots,k\): choose \(\mu_j=x\) with probability proportional to \[D(x)^2=\min_{\ell<j}\lVert x-\mu_\ell\rVert^2 .\]
- Run Lloyd’s algorithm from those centres.
Theorem (Arthur & Vassilvitskii). With \(k\)-means++ seeding, \[\mathbb{E}[\phi]\ \le\ 8(\ln k+2)\,\phi_{\mathrm{OPT}},\] where \(\phi\) is the resulting objective and \(\phi_{\mathrm{OPT}}\) the global optimum.
This is an \(O(\log k)\) expected approximation guarantee, and it holds before Lloyd’s algorithm even runs, the seeding alone is competitive. Nothing comparable is true of uniform random seeding, whose worst case is unbounded.
kmeanspp_init <- function(X, k, seed = NULL) {
if (!is.null(seed)) set.seed(seed)
X <- as.matrix(X); n <- nrow(X)
centres <- matrix(NA_real_, k, ncol(X))
centres[1, ] <- X[sample(n, 1), ]
if (k > 1) for (j in 2:k) {
# D(x)^2 = squared distance to the NEAREST chosen centre
d2 <- apply(centres[1:(j - 1), , drop = FALSE], 1,
function(c0) colSums((t(X) - c0)^2))
d2 <- if (is.null(dim(d2))) d2 else apply(d2, 1, min)
# Sample with probability proportional to D(x)^2
centres[j, ] <- X[sample(n, 1, prob = d2 / sum(d2)), ]
}
centres
}set.seed(111)
compare_init <- vapply(1:200, function(r) {
c(random = kmeans(Xg, 3, nstart = 1)$tot.withinss,
kmeanspp = kmeans(Xg, kmeanspp_init(Xg, 3, seed = 300 + r),
algorithm = "Hartigan-Wong")$tot.withinss)
}, numeric(2))
data.frame(
initialization = c("Uniform random (nstart = 1)", "k-means++ (single seed)",
"Uniform random (nstart = 25)"),
median_wss = round(c(median(compare_init["random", ]),
median(compare_init["kmeanspp", ]),
kmeans(Xg, 3, nstart = 25)$tot.withinss), 2),
worst_wss = round(c(max(compare_init["random", ]),
max(compare_init["kmeanspp", ]), NA), 2),
pct_reaching_best = round(100 * c(
mean(abs(compare_init["random", ] - min(compare_init)) < 1e-6),
mean(abs(compare_init["kmeanspp", ] - min(compare_init)) < 1e-6), 100), 1))\(k\)-means++ reaches the best
solution far more often than uniform seeding from a single start. Note
the third row, though: nstart = 25 with plain
random seeding also finds it, and does so with one line and no
custom code. In practice, use nstart; understand \(k\)-means++ because it is what the
guarantee is about, and because it is the default in
scikit-learn and most large-scale implementations where 25
restarts would be prohibitive.
10 Choosing k
Neither a large nor a small \(k\) is inherently better. Too many clusters and the groups are too specific to interpret; too few and they are too broad to be useful. Four methods, in increasing order of principle.
Common misconception: “\(k=\sqrt{n/2}\) is the rule of thumb.” There is no such rule with any theoretical basis, it is folklore, the clustering analogue of \(\sqrt n\) for kNN (Chapter 5, §5.6.2). It takes no account of the data’s actual structure, dimension, or separation. Use it, if at all, as an upper bound on where to start searching.
10.1 Elbow
Plot the within-cluster sum of squares against \(k\). WSS decreases monotonically, adding a cluster can never increase it, so the question is where the rate of decrease slackens.
wss_curve <- function(X, kmax = 10, nstart = 25) {
data.frame(k = 1:kmax,
wss = vapply(1:kmax, function(k)
kmeans(X, k, nstart = nstart)$tot.withinss, numeric(1)))
}
set.seed(121)
ec <- wss_curve(Xg)
ec# Kink detection: the point of maximum perpendicular distance from the chord
# joining the first and last points of the curve.
elbow_point <- function(k, wss) {
p1 <- c(k[1], wss[1]); p2 <- c(k[length(k)], wss[length(wss)])
v <- p2 - p1; v <- v / sqrt(sum(v^2)) # unit chord direction
d <- vapply(seq_along(k), function(i) {
w <- c(k[i], wss[i]) - p1
sqrt(sum((w - sum(w * v) * v)^2)) # || w - proj_v(w) ||
}, numeric(1))
list(k = k[which.max(d)], distances = d)
}
eb <- elbow_point(ec$k, ec$wss)
c(elbow_k = eb$k)#> elbow_k
#> 2The function is pure, it takes everything it needs as arguments and reads no global state. The geometry is the orthogonal-projection identity from Chapter 3, §3.4.2: for the chord direction \(\hat v\), the perpendicular component of \(w\) is \(w-\langle w,\hat v\rangle\hat v\), and its norm is the distance from the point to the line.
ec$dist <- eb$distances
chord <- data.frame(k = range(ec$k), wss = ec$wss[c(1, nrow(ec))])
ggplot(ec, aes(k, wss)) +
geom_line(data = chord, linewidth = 0.8, color = "grey55", linetype = "dashed") +
geom_line(linewidth = 1, color = "steelblue") +
geom_point(size = 2.2) +
geom_point(data = ec[ec$k == eb$k, ], color = "firebrick", size = 4.5) +
scale_x_continuous(breaks = 1:10) +
labs(title = "Elbow method on real within-cluster sums of squares",
subtitle = sprintf("Dashed: chord from first to last point. Red: maximum perpendicular distance, at k = %d", eb$k),
x = "Number of clusters k", y = "Total within-cluster sum of squares") +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(ec, x = ~k, y = ~wss, type = "scatter", mode = "lines+markers",
name = "WSS") |>
add_lines(x = chord$k, y = chord$wss, name = "Reference chord",
line = list(dash = "dash")) |>
add_markers(x = eb$k, y = ec$wss[ec$k == eb$k], name = "Elbow",
marker = list(size = 16, color = "red")) |>
layout(title = "Within-cluster sum of squares vs. k",
xaxis = list(title = "k"), yaxis = list(title = "WSS"))The elbow is a heuristic: on smoothly decaying curves the kink is ambiguous, and different readers pick different points. It is a starting place, not an answer.
10.2 Silhouette
sil_curve <- function(X, kmax = 10, nstart = 25) {
D <- dist(X)
data.frame(k = 2:kmax,
sil = vapply(2:kmax, function(k)
mean(silhouette(kmeans(X, k, nstart = nstart)$cluster, D)[, "sil_width"]),
numeric(1)))
}
set.seed(131)
sc <- sil_curve(Xg)
sc |> mutate(sil = round(sil, 4), verdict = as.character(sil_verdict(sil)))ggplot(sc, aes(k, sil)) +
geom_hline(yintercept = c(0.26, 0.51, 0.71), linetype = "dashed",
color = "grey60") +
ggplot2::annotate("text", x = 9.5, y = c(0.28, 0.53, 0.73), size = 3, hjust = 1,
color = "grey40",
label = c("weak", "reasonable", "strong")) +
geom_line(linewidth = 1, color = "steelblue") + geom_point(size = 2.2) +
geom_point(data = sc[which.max(sc$sil), ], color = "firebrick", size = 4.5) +
scale_x_continuous(breaks = 2:10) +
labs(title = "Average silhouette against k, with the standard thresholds",
x = "k", y = "Average silhouette width") +
theme_dspa()
10.3 The gap statistic
The elbow and silhouette compare \(k\) against other values of \(k\). The gap statistic (Tibshirani, Walther & Hastie, 2001) compares the observed WSS against what would be expected with no structure at all, which is the question you actually want answered.
\[\boxed{\;\mathrm{Gap}(k)=\mathbb{E}^{*}_n\big[\log W_k\big]-\log W_k\;}\]
where \(\mathbb{E}^{*}_n\) is the expectation under a reference null distribution, usually uniform over the bounding box of the data, or over the box aligned with its principal components. Estimate it by generating \(B\) reference datasets, clustering each, and averaging.
Choose the smallest \(k\) satisfying the 1-standard-error rule:
\[\hat k=\min\Big\{k:\ \mathrm{Gap}(k)\ \ge\ \mathrm{Gap}(k+1)-s_{k+1}\Big\}, \qquad s_k=\mathrm{sd}_k\sqrt{1+1/B}.\]
Crucially, the gap statistic can select \(k=1\), which no other method here can. That is its most valuable property: it is the only one that can report there are no clusters.
set.seed(141)
gap_real <- clusGap(Xg, FUN = kmeans, nstart = 25, K.max = 8, B = 50)
gap_noise <- clusGap(as.matrix(noise), FUN = kmeans, nstart = 25, K.max = 8, B = 50)
c(structured_data_selects_k = maxSE(gap_real$Tab[, "gap"], gap_real$Tab[, "SE.sim"],
method = "firstSEmax"),
uniform_noise_selects_k = maxSE(gap_noise$Tab[, "gap"], gap_noise$Tab[, "SE.sim"],
method = "firstSEmax"))#> structured_data_selects_k uniform_noise_selects_k
#> 3 1gd <- bind_rows(
data.frame(k = 1:8, gap = gap_real$Tab[, "gap"], se = gap_real$Tab[, "SE.sim"],
data = "Three real groups"),
data.frame(k = 1:8, gap = gap_noise$Tab[, "gap"], se = gap_noise$Tab[, "SE.sim"],
data = "Uniform noise"))
ggplot(gd, aes(k, gap)) +
geom_errorbar(aes(ymin = gap - se, ymax = gap + se), width = 0.15,
color = "grey55") +
geom_line(linewidth = 1, color = "steelblue") + geom_point(size = 2.2) +
facet_wrap(~ data, scales = "free_y") +
scale_x_continuous(breaks = 1:8) +
labs(title = "The gap statistic compares against a no-structure null",
subtitle = "On noise the gap does not rise, so k = 1 is selected -- no other method can say that",
x = "k", y = "Gap(k)") +
theme_dspa(10)
The right-hand panel is the point: on uniform noise the gap statistic correctly returns \(k=1\), while the elbow and silhouette would both hand back some \(k\ge2\) and a partition to interpret.
The reference null and the observed curve together form a surface over \((k,\text{draw})\), which shows the spread the observed value is compared against:
set.seed(151)
B_ref <- 30; K_ref <- 8
rng <- apply(Xg, 2, range)
ref_logW <- matrix(NA_real_, B_ref, K_ref)
for (b in seq_len(B_ref)) {
Xb <- cbind(runif(nrow(Xg), rng[1, 1], rng[2, 1]),
runif(nrow(Xg), rng[1, 2], rng[2, 2]))
for (k in 1:K_ref) ref_logW[b, k] <- log(kmeans(Xb, k, nstart = 10)$tot.withinss)
}
obs_logW <- vapply(1:K_ref, \(k) log(kmeans(Xg, k, nstart = 25)$tot.withinss), numeric(1))
plot_ly() |>
add_surface(x = 1:K_ref, y = 1:B_ref, z = ref_logW, opacity = 0.75,
showscale = FALSE, colorscale = "Greys",
name = "Reference null draws") |>
add_trace(x = 1:K_ref, y = rep(B_ref / 2, K_ref), z = obs_logW,
type = "scatter3d", mode = "lines+markers", name = "Observed data",
line = list(width = 8, color = "#D8433B"),
marker = list(size = 4, color = "#D8433B")) |>
layout(title = "Observed log(WSS) against the reference null surface",
scene = list(xaxis = list(title = "k"),
yaxis = list(title = "Reference draw b"),
zaxis = list(title = "log(WSS)")))The gap is the vertical separation between the red curve and the grey surface. Where the surface and the curve run parallel, adding clusters buys no more than it would on noise.
11 When k-means fails
\(k\)-means minimizes squared Euclidean distance to centroids. That single fact implies three assumptions, and violating any of them produces a confidently wrong answer.
| Assumption | Why it follows | What breaks it |
|---|---|---|
| Spherical clusters | Squared Euclidean distance has circular level sets | Elongated, curved, or nested shapes |
| Similar variances | A single centroid per cluster, no scale parameter | One tight cluster beside one diffuse one |
| Similar sizes | The objective is a sum over points, so large clusters dominate | A cluster of 500 beside one of 20 |
And one more, which is not an assumption but a property: \(k\)-means always returns \(k\) non-empty clusters, whatever the data.
set.seed(161)
n_f <- 500
# 1. Elongated / anisotropic
Xa <- MASS::mvrnorm(n_f, c(0, 0), matrix(c(6, 5.6, 5.6, 6), 2))
Xa <- rbind(Xa, MASS::mvrnorm(n_f, c(4, -4), matrix(c(6, 5.6, 5.6, 6), 2)))
ya <- rep(1:2, each = n_f)
# 2. Nested rings
th <- runif(n_f, 0, 2 * pi); r1 <- sqrt(runif(n_f, 0, 1))
th2 <- runif(n_f, 0, 2 * pi); r2 <- sqrt(runif(n_f, 6, 9))
Xb <- rbind(cbind(r1 * cos(th), r1 * sin(th)),
cbind(r2 * cos(th2), r2 * sin(th2)))
yb <- rep(1:2, each = n_f)
# 3. Unequal density and size
Xc <- rbind(MASS::mvrnorm(n_f, c(0, 0), diag(2) * 0.3),
MASS::mvrnorm(50, c(5, 5), diag(2) * 4))
yc <- c(rep(1, n_f), rep(2, 50))
# 4. No structure
Xd <- matrix(runif(2 * n_f * 2), ncol = 2)
yd <- rep(1, 2 * n_f)
cases <- list(`1. Anisotropic` = list(Xa, ya, 2),
`2. Nested rings` = list(Xb, yb, 2),
`3. Unequal density and size` = list(Xc, yc, 2),
`4. No structure` = list(Xd, yd, 2))
fail_df <- bind_rows(lapply(names(cases), function(nm) {
X <- cases[[nm]][[1]]; y <- cases[[nm]][[2]]; k <- cases[[nm]][[3]]
km <- kmeans(X, k, nstart = 25)
data.frame(x = X[, 1], y = X[, 2], truth = factor(y),
kmeans = factor(km$cluster), panel = nm)
}))
p_truth <- ggplot(fail_df, aes(x, y, color = truth)) +
geom_point(size = 0.6, alpha = 0.7) +
facet_wrap(~ panel, nrow = 1, scales = "free") +
scale_color_brewer(palette = "Dark2", guide = "none") +
labs(title = "True structure", x = NULL, y = NULL) + theme_dspa(8)
p_km <- ggplot(fail_df, aes(x, y, color = kmeans)) +
geom_point(size = 0.6, alpha = 0.7) +
facet_wrap(~ panel, nrow = 1, scales = "free") +
scale_color_brewer(palette = "Set1", guide = "none") +
labs(title = "k-means partition", x = NULL, y = NULL) + theme_dspa(8)
p_truth / p_km
do.call(rbind, lapply(names(cases), function(nm) {
X <- cases[[nm]][[1]]; y <- cases[[nm]][[2]]; k <- cases[[nm]][[3]]
km <- kmeans(X, k, nstart = 25)
data.frame(case = nm,
ARI = round(mclust::adjustedRandIndex(km$cluster, y), 4),
silhouette = round(mean(silhouette(km$cluster, dist(X))[, "sil_width"]), 4))
}))Read the last two columns together. In case 4 the silhouette is respectable and the ARI is zero, because there is nothing to recover. In case 2 the silhouette is also positive while the ARI is near zero: \(k\)-means has cut the concentric rings in half like a pie, producing two compact clusters that are completely wrong.
Common misconception: “a high silhouette means the clustering is correct.” Silhouette measures compactness and separation under the metric being used, and \(k\)-means optimizes exactly that quantity. A good silhouette therefore partly measures how well the algorithm did its own job, not whether the partition corresponds to anything real. Case 2 above is the clean demonstration: cutting concentric rings into pie slices produces compact, well-separated clusters with a respectable silhouette and an ARI of essentially zero.
Always pair an internal index with external labels, a stability check, or a method making different assumptions. When all three agree, you have a finding; when the internal index is the only thing supporting the partition, you have an artifact.
12 k-medoids
PAM (Partitioning Around Medoids) replaces means with medoids, actual observations, and minimizes the sum of dissimilarities rather than squared distances:
\[\arg\min_{M\subset\{x_1,\dots,x_n\},|M|=k}\ \sum_{i=1}^{n}\min_{m\in M}d(x_i,m).\]
Three consequences follow.
It is robust. Minimizing \(\sum d\) rather than \(\sum d^2\) gives outliers linear rather than quadratic influence, the median/mean contrast of Chapter 2, §2.3.2 applied to centres.
It works from any dissimilarity matrix, so it pairs directly with Gower distance for mixed-type data. There is no need for a centroid mean, which is undefined for a nominal variable.
The medoids are real cases, which is often far more interpretable than a centroid that corresponds to no actual observation.
Cost: \(O(k(n-k)^2)\) per iteration
for the classical swap phase, versus \(O(nkd)\) for \(k\)-means. cluster::clara()
subsamples for large \(n\).
set.seed(171)
Xr <- rbind(MASS::mvrnorm(150, c(0, 0), diag(2)),
MASS::mvrnorm(150, c(6, 0), diag(2)))
yr <- rep(1:2, each = 150)
Xout <- rbind(Xr, matrix(c(30, 30, 32, 28, 29, 31), ncol = 2, byrow = TRUE))
yout <- c(yr, rep(NA, 3))
km_out <- kmeans(Xout, 2, nstart = 25)
pam_out <- pam(Xout, 2)
data.frame(
method = c("k-means", "PAM"),
ARI_excluding_outliers = round(c(
mclust::adjustedRandIndex(km_out$cluster[1:300], yr),
mclust::adjustedRandIndex(pam_out$clustering[1:300], yr)), 4))bind_rows(
data.frame(x = Xout[, 1], y = Xout[, 2],
cl = factor(km_out$cluster), method = "k-means"),
data.frame(x = Xout[, 1], y = Xout[, 2],
cl = factor(pam_out$clustering), method = "PAM")) |>
ggplot(aes(x, y, color = cl)) +
geom_point(size = 1.1, alpha = 0.8) +
facet_wrap(~ method) +
scale_color_brewer(palette = "Set1", guide = "none") +
labs(title = "Three outliers, two methods",
subtitle = "k-means devotes a centroid to the outliers and merges the real groups; PAM does not",
x = NULL, y = NULL) +
theme_dspa(10)
Three contaminating points out of 303 are enough to make \(k\)-means merge the two real clusters and spend its second centroid on the contamination. PAM recovers the structure.
13 Case study: divorce and young adults
A longitudinal study of young adults (initially aged 18–23) whose parents divorced within 15 months of the first wave (1990–91). The case-study subset has 47 respondents.
Variables and their measurement scales, this table is the analysis, not preamble, because the scale determines the distance:
| Variable | Description | Scale |
|---|---|---|
DIVYEAR |
Year parents divorced (89 or 90) | Binary |
momint |
Mother intimacy (1 extremely close … 4 not close) | Ordinal |
dadint |
Father intimacy (1 … 4) | Ordinal |
momclose |
Closeness to mother (1 … 4) | Ordinal |
depression |
Feelings of depression (1 often … 4 never) | Ordinal |
livewithmom |
1 mother only, 2 father only, 3 both; 9 = missing | Nominal + sentinel |
gethitched |
1 marry soon, 2 marry sometime, 3 never; 8 = don’t know | Ordinal + sentinel |
Two of the seven are not interval-scaled, and two carry sentinel codes.
divorce <- dspa_read(
"https://umich.instructure.com/files/399118/download?download_frd=1",
"CaseStudy01_Divorce_YoungAdults_Data.csv")
str(divorce)#> 'data.frame': 47 obs. of 7 variables:
#> $ DIVYEAR : int 90 90 89 90 89 90 90 90 90 90 ...
#> $ momint : int 3 3 1 3 1 1 1 3 2 3 ...
#> $ dadint : int 2 2 3 4 3 4 1 2 2 1 ...
#> $ momclose : int 2 3 1 3 2 1 1 2 2 2 ...
#> $ depression : int 3 2 4 3 1 2 4 2 4 2 ...
#> $ livewithmom: int 1 2 1 1 1 1 2 1 1 1 ...
#> $ gethitched : int 3 2 2 3 2 2 2 2 2 2 ...#> DIVYEAR momint dadint momclose depression
#> Min. :89.0 Min. :1.00 Min. :1.00 Min. :1.00 Min. :1.00
#> 1st Qu.:89.0 1st Qu.:1.00 1st Qu.:2.00 1st Qu.:1.00 1st Qu.:2.00
#> Median :90.0 Median :1.00 Median :2.00 Median :2.00 Median :3.00
#> Mean :89.7 Mean :1.81 Mean :2.49 Mean :1.81 Mean :2.85
#> 3rd Qu.:90.0 3rd Qu.:3.00 3rd Qu.:3.00 3rd Qu.:2.00 3rd Qu.:4.00
#> Max. :90.0 Max. :4.00 Max. :4.00 Max. :4.00 Max. :4.00
#> livewithmom gethitched
#> Min. :1.00 Min. :1.00
#> 1st Qu.:1.00 1st Qu.:2.00
#> Median :1.00 Median :2.00
#> Mean :1.49 Mean :2.21
#> 3rd Qu.:2.00 3rd Qu.:2.00
#> Max. :9.00 Max. :8.0013.1 Handling the sentinel codes
#> livewithmom_eq_9 gethitched_eq_8
#> 1 1Both are missing-data codes, not measurements. Treating either as a number places those respondents far from everyone else once the variable is scaled.
# A CODED rule, applied to every affected row -- not a hand edit of one cell.
# Respondents closer to their father than their mother are inferred to live
# with the father; otherwise with the mother. Higher `int` = LESS close.
div <- divorce |>
mutate(
livewithmom = ifelse(livewithmom == 9,
ifelse(dadint < momint, 2L, 1L),
livewithmom),
gethitched = ifelse(gethitched == 8, NA_integer_, gethitched))
c(imputed_livewithmom = sum(divorce$livewithmom == 9),
gethitched_now_NA = sum(is.na(div$gethitched)))#> imputed_livewithmom gethitched_now_NA
#> 1 1# Encode each variable at its true measurement scale
div_typed <- div |>
transmute(
divyear = factor(DIVYEAR, levels = c(89, 90)), # binary
momint = factor(momint, levels = 1:4, ordered = TRUE), # ordinal
dadint = factor(dadint, levels = 1:4, ordered = TRUE),
momclose = factor(momclose, levels = 1:4, ordered = TRUE),
depression = factor(depression, levels = 1:4, ordered = TRUE),
livewithmom = factor(livewithmom, levels = 1:3, # NOMINAL
labels = c("mother", "father", "both")),
gethitched = factor(gethitched, levels = 1:3, ordered = TRUE))
str(div_typed)#> 'data.frame': 47 obs. of 7 variables:
#> $ divyear : Factor w/ 2 levels "89","90": 2 2 1 2 1 2 2 2 2 2 ...
#> $ momint : Ord.factor w/ 4 levels "1"<"2"<"3"<"4": 3 3 1 3 1 1 1 3 2 3 ...
#> $ dadint : Ord.factor w/ 4 levels "1"<"2"<"3"<"4": 2 2 3 4 3 4 1 2 2 1 ...
#> $ momclose : Ord.factor w/ 4 levels "1"<"2"<"3"<"4": 2 3 1 3 2 1 1 2 2 2 ...
#> $ depression : Ord.factor w/ 4 levels "1"<"2"<"3"<"4": 3 2 4 3 1 2 4 2 4 2 ...
#> $ livewithmom: Factor w/ 3 levels "mother","father",..: 1 2 1 1 1 1 2 1 1 1 ...
#> $ gethitched : Ord.factor w/ 3 levels "1"<"2"<"3": 3 2 2 3 2 2 2 2 2 2 ...gethitched retains NA, which Gower distance
handles by down-weighting that variable for the
affected pairs rather than imputing a value.
13.2 Two distances, two answers
# Gower: respects binary / ordinal / nominal types and missingness
d_gower <- daisy(div_typed, metric = "gower")
# Scaled Euclidean on the raw numeric codes, for contrast
d_euclid <- dist(scale(div[, names(divorce)]))
D <- if (DIVORCE_DISTANCE == "gower") d_gower else d_euclid
c(distance_used = DIVORCE_DISTANCE,
gower_range = paste(round(range(d_gower), 3), collapse = " to "),
euclid_range = paste(round(range(d_euclid), 3), collapse = " to "),
correlation_between_them = round(cor(as.numeric(d_gower),
as.numeric(d_euclid)), 4))#> distance_used gower_range euclid_range
#> "gower" "0 to 0.786" "0 to 7.329"
#> correlation_between_them
#> "0.892"The two dissimilarity matrices correlate only moderately, they are ordering the same 1,081 pairs differently, so they will produce different partitions.
set.seed(181)
pam_g <- pam(d_gower, k = 3, diss = TRUE)
pam_e <- pam(d_euclid, k = 3, diss = TRUE)
data.frame(
distance = c("Gower (typed)", "Scaled Euclidean (naive)"),
avg_silhouette = round(c(pam_g$silinfo$avg.width, pam_e$silinfo$avg.width), 4),
verdict = as.character(sil_verdict(c(pam_g$silinfo$avg.width,
pam_e$silinfo$avg.width))),
agreement_ARI = round(mclust::adjustedRandIndex(pam_g$clustering,
pam_e$clustering), 4))ggplot(data.frame(gower = as.numeric(d_gower), euclid = as.numeric(d_euclid)),
aes(gower, euclid)) +
geom_point(alpha = 0.15, size = 0.8, color = "steelblue") +
geom_smooth(method = "lm", formula = y ~ x, color = "firebrick", se = FALSE) +
labs(title = "The same 1,081 pairs, ranked by two distances",
subtitle = "Scatter off the line means the two metrics disagree about which respondents are alike",
x = "Gower distance (types respected)",
y = "Scaled Euclidean distance (codes as numbers)") +
theme_dspa()
13.3 What the validation says
set.seed(191)
k_range <- 2:10
div_indices <- do.call(rbind, lapply(k_range, function(k) {
p <- pam(D, k = k, diss = TRUE)
data.frame(k = k, silhouette = p$silinfo$avg.width)
}))
div_indices |> mutate(silhouette = round(silhouette, 4),
verdict = as.character(sil_verdict(silhouette)))ggplot(div_indices, aes(k, silhouette)) +
geom_hline(yintercept = 0.26, linetype = "dashed", color = "firebrick") +
ggplot2::annotate("text", x = 10, y = 0.275, hjust = 1, size = 3.2, color = "firebrick",
label = "threshold for weak structure") +
geom_line(linewidth = 1, color = "steelblue") + geom_point(size = 2.2) +
scale_x_continuous(breaks = k_range) +
coord_cartesian(ylim = c(0, 0.45)) +
labs(title = "Average silhouette across k for the divorce data",
subtitle = "Every value sits at or below the threshold for even weak structure",
x = "k", y = "Average silhouette width") +
theme_dspa()
set.seed(201)
## OLD: not safe for missing data ...
# # The gap statistic can return k = 1, i.e. "no clusters"
# div_num <- scale(div[, names(divorce)])
# gap_div <- clusGap(div_num, FUN = kmeans, nstart = 25, K.max = 8, B = 60)
# Filter the original data to complete cases for those specific columns
div_clean <- na.omit(div[, names(divorce)])
# Scale the clean data
div_num <- scale(div_clean)
# Now run clusGap again
gap_div <- clusGap(div_num, FUN = kmeans, nstart = 25, K.max = 8, B = 60)
k_gap <- maxSE(gap_div$Tab[, "gap"], gap_div$Tab[, "SE.sim"],
method = "firstSEmax")
c(gap_statistic_selects_k = k_gap)#> gap_statistic_selects_k
#> 1#> gap SE.sim
#> [1,] 0.2768 0.0261
#> [2,] 0.2760 0.0260
#> [3,] 0.2768 0.0281
#> [4,] 0.2978 0.0294
#> [5,] 0.3139 0.0296
#> [6,] 0.3138 0.0298
#> [7,] 0.3337 0.0294
#> [8,] 0.3468 0.0292set.seed(211)
stab_div <- fpc::clusterboot(div_num, B = 50, clustermethod = fpc::kmeansCBI,
krange = 3, seed = 211, count = FALSE)
data.frame(cluster = 1:3,
bootstrap_jaccard = round(stab_div$bootmean, 4),
verdict = cut(stab_div$bootmean, c(-Inf, 0.6, 0.75, 0.85, Inf),
labels = c("not reproducible", "patchy",
"stable", "highly stable"), right = FALSE))Three independent diagnostics agree. The silhouette never rises above the threshold for even weak structure at any \(k\); the gap statistic selects a small \(k\) against a no-structure null; and the bootstrap Jaccard values fall below the reproducibility threshold.
This is a finding, and it should be reported as one. The divorce data does not exhibit cluster structure in these seven variables. That conclusion is more useful than a partition would be, because it tells you something true about the population: on these measures, young adults of recently divorced parents form a continuum rather than discrete types.
What must not happen next is a narrative about the centroids. \(k\)-means will produce three well-separated centres on this data, it produces them on uniform noise (§8.1), and a bar chart of those centres will look interpretable. Naming the groups “emotional,” “naive,” and “independent” and describing their attitudes to family would be describing an artifact of the algorithm.
set.seed(221)
km_div <- kmeans(div_num, 3, nstart = 25)
cent_div <- as.data.frame(t(km_div$centers))
colnames(cent_div) <- paste0("Cluster ", 1:3)
cent_div$variable <- rownames(cent_div)
cent_div |>
pivot_longer(-variable, names_to = "cluster", values_to = "centre") |>
ggplot(aes(variable, centre, fill = cluster)) +
geom_col(position = "dodge") +
scale_fill_brewer(palette = "Set1") +
labs(title = "Cluster centroids -- shown for illustration, not interpretation",
subtitle = sprintf("Average silhouette %.3f: below the threshold for substantial structure. These separations are guaranteed by the algorithm",
mean(silhouette(km_div$cluster, dist(div_num))[, "sil_width"])),
x = NULL, y = "Standardized centroid value", fill = NULL) +
theme_dspa(9) +
theme(axis.text.x = element_text(angle = 30, hjust = 1))
# --- Interactive equivalent ------------------------------------------------
df <- as.data.frame(t(km_div$centers))
colnames(df) <- paste0("Cluster", 1:3)
plot_ly(df, x = rownames(df), y = ~Cluster1, type = "bar", name = "Cluster 1") |>
add_trace(y = ~Cluster2, name = "Cluster 2") |>
add_trace(y = ~Cluster3, name = "Cluster 3") |>
layout(title = "Cluster centroids (illustrative only)",
yaxis = list(title = "Standardized centroid value"), barmode = "group")13.4 What to do instead
When clustering finds no structure, the productive next steps are dimension reduction and association analysis, not a different clustering algorithm.
pca_div <- prcomp(div_num)
varexp <- pca_div$sdev^2 / sum(pca_div$sdev^2)
c(PC1_variance = round(varexp[1], 4), PC2_variance = round(varexp[2], 4),
cumulative_first_two = round(sum(varexp[1:2]), 4))#> PC1_variance PC2_variance cumulative_first_two
#> 0.2617 0.1907 0.4523ggplot(data.frame(PC1 = pca_div$x[, 1], PC2 = pca_div$x[, 2],
depression = factor(div_clean$depression)),
aes(PC1, PC2, color = depression)) +
geom_point(size = 2.4, alpha = 0.85) +
scale_color_viridis_d(option = "plasma", name = "Depression\n(1 often - 4 never)") +
labs(title = "The divorce data in its first two principal components",
subtitle = sprintf("PC1 %.0f%%, PC2 %.0f%% of variance. A continuum, not separated groups",
100 * varexp[1], 100 * varexp[2]),
x = "PC1", y = "PC2") +
theme_dspa()
The scatter is a single cloud with a gradient, not discrete islands, which is what the silhouette, gap, and stability results were all reporting, now visible directly.
14 Case study: services utilization by trauma-exposed children
Services Utilization by Trauma-Exposed Children examines associations between post-traumatic psychopathology and service use. Unlike the divorce data, this one carries a known grouping, the type of trauma exposure, which is held back from the algorithm and used only for external validation.
Variables. id; sex (1
female, 0 male); age (0–18); race (black,
white, hispanic, other); cmt (maltreatment exposure, 0/1);
traumatype (sexabuse, physabuse, neglect, psychabuse,
dvexp); ptsd (0/1); dissoc (0–11);
service (count, 0–19).
trauma <- dspa_read(
"https://umich.instructure.com/files/399129/download?download_frd=1",
"Case_04_ChildTrauma_Data.csv", sep = " ")
str(trauma); dim(trauma)#> 'data.frame': 1000 obs. of 9 variables:
#> $ id : int 1 2 3 4 5 6 7 8 9 10 ...
#> $ sex : int 1 1 0 0 1 0 0 1 0 1 ...
#> $ age : int 6 14 6 11 7 9 12 9 9 13 ...
#> $ ses : int 0 0 0 0 0 0 0 0 1 0 ...
#> $ race : chr "black" "black" "black" "black" ...
#> $ traumatype: chr "sexabuse" "sexabuse" "sexabuse" "sexabuse" ...
#> $ ptsd : int 1 0 0 0 1 1 1 0 1 1 ...
#> $ dissoc : int 1 0 1 1 1 0 1 1 1 0 ...
#> $ service : int 17 12 9 11 15 6 9 10 11 13 ...#> [1] 1000 9#>
#> dvexp neglect physabuse psychabuse sexabuse
#> 250 350 100 200 100# Columns selected BY NAME. The outcome-like variable `traumatype` and the
# identifier are excluded from the features; traumatype is retained separately
# for external validation only.
truth_trauma <- factor(trauma$traumatype)
trauma_features <- trauma |>
select(-any_of(c("id", "traumatype"))) |>
mutate(race = factor(race))
c(features = ncol(trauma_features),
names = paste(names(trauma_features), collapse = ", "))#> features
#> "7"
#> names
#> "sex, age, ses, race, ptsd, dissoc, service"# One-hot encode race for the numeric pipeline; Gower would take the factor
# directly, but k-means requires numeric input.
tr_num <- model.matrix(~ . - 1, data = trauma_features)
tr_z <- scale(tr_num)
c(rows = nrow(tr_z), numeric_columns = ncol(tr_z))#> rows numeric_columns
#> 1000 10set.seed(231)
tr_curve <- do.call(rbind, lapply(2:10, function(k) {
km <- kmeans(tr_z, k, nstart = 25)
data.frame(k = k,
silhouette = mean(silhouette(km$cluster, dist(tr_z))[, "sil_width"]),
ARI_vs_traumatype = mclust::adjustedRandIndex(km$cluster, truth_trauma))
}))
tr_curve |> mutate(across(where(is.numeric), \(z) round(z, 4)))tr_curve |>
pivot_longer(-k, names_to = "measure", values_to = "value") |>
ggplot(aes(k, value, color = measure)) +
geom_line(linewidth = 1) + geom_point(size = 2.2) +
scale_x_continuous(breaks = 2:10) +
scale_color_manual(values = c(silhouette = "#3B7DD8",
ARI_vs_traumatype = "#D8433B")) +
labs(title = "Internal and external validation disagree about k",
subtitle = "Silhouette rewards compactness; ARI rewards recovering the true trauma types",
x = "k", y = NULL, color = NULL) +
theme_dspa()
The two curves peak at different \(k\). That disagreement is informative: silhouette asks “are these clusters compact and separated?” while ARI asks “do they correspond to the known groups?” A partition can be geometrically tidy and scientifically irrelevant, and this data shows both curves at once.
k_ari <- tr_curve$k[which.max(tr_curve$ARI_vs_traumatype)]
set.seed(241)
km_tr <- kmeans(tr_z, k_ari, nstart = 25)
c(k_used = k_ari,
n_true_classes = nlevels(truth_trauma),
ARI = round(mclust::adjustedRandIndex(km_tr$cluster, truth_trauma), 4),
avg_silhouette = round(mean(silhouette(km_tr$cluster, dist(tr_z))[, "sil_width"]), 4))#> k_used n_true_classes ARI avg_silhouette
#> 4.0000 5.0000 0.3003 0.2482#> traumatype
#> cluster dvexp neglect physabuse psychabuse sexabuse
#> 1 43 279 0 166 0
#> 2 31 71 19 34 10
#> 3 100 0 0 0 0
#> 4 76 0 81 0 90as.data.frame(tab_tr) |>
mutate(prop = Freq / sum(Freq), .by = traumatype) |>
ggplot(aes(traumatype, cluster, fill = prop)) +
geom_tile(color = "white") +
geom_text(aes(label = Freq), size = 3.4) +
scale_fill_viridis_c(option = "mako", direction = -1,
labels = scales::percent, name = "Share of\ntrauma type") +
labs(title = "Recovered clusters against the held-back trauma types",
subtitle = sprintf("ARI = %.3f. A perfect recovery would place all mass in one cell per column",
mclust::adjustedRandIndex(km_tr$cluster, truth_trauma)),
x = "True trauma type", y = "Cluster") +
theme_dspa(10)
Read the ARI before the table. Cluster labels are arbitrary permutations, so a raw cross-tabulation cannot be summarized by eye, one has to solve a matching problem, and even a perfect matching may reflect agreement expected by chance. ARI does both corrections at once, and its value here says the recovery is partial: some trauma types separate cleanly, others do not.
# Which trauma types are recoverable, and which are not?
as.data.frame(tab_tr) |>
summarise(largest_cluster_share = max(Freq) / sum(Freq),
n = sum(Freq), .by = traumatype) |>
arrange(desc(largest_cluster_share)) |>
mutate(largest_cluster_share = round(largest_cluster_share, 3))The types that concentrate in a single cluster are the ones the measured features distinguish. The types that spread across clusters are not distinguishable by these variables, which is a substantive finding about the instrument, not a failure of the algorithm.
cent_tr <- as.data.frame(t(km_tr$centers))
colnames(cent_tr) <- paste0("Cluster ", seq_len(ncol(cent_tr)))
cent_tr$variable <- rownames(cent_tr)
cent_tr |>
pivot_longer(-variable, names_to = "cluster", values_to = "centre") |>
ggplot(aes(variable, centre, fill = cluster)) +
geom_col(position = "dodge") +
scale_fill_brewer(palette = "Set2") +
labs(title = "Cluster centroids for the trauma data",
subtitle = sprintf("Interpretable here because ARI = %.2f confirms the clusters track real groups",
mclust::adjustedRandIndex(km_tr$cluster, truth_trauma)),
x = NULL, y = "Standardized centroid value", fill = NULL) +
theme_dspa(9) +
theme(axis.text.x = element_text(angle = 35, hjust = 1))
# --- Interactive equivalent ------------------------------------------------
df <- as.data.frame(t(km_tr$centers))
colnames(df) <- paste0("Cluster", seq_len(ncol(df)))
p <- plot_ly(df, x = rownames(df), y = ~Cluster1, type = "bar", name = "Cluster 1")
for (j in 2:ncol(df))
p <- add_trace(p, y = df[[j]], name = paste("Cluster", j))
p |> layout(title = "Explicating derived cluster labels",
yaxis = list(title = "Cluster centres"), barmode = "group")Here the centroids may be interpreted, because the external validation established that the clusters track something real. The contrast with §8.10 is the methodological point of both case studies: the same plot is informative in one and misleading in the other, and only the validation distinguishes them.
14.1 Feature selection for clustering
Clustering is subject to the curse of dimensionality (§8.2.2), so selecting features matters, and without labels, the usual supervised criteria are unavailable. Current approaches:
- Model-based selection. For Gaussian mixtures,
clustvarselperforms greedy variable selection using BIC comparisons between models that do and do not cluster on each variable. - Sparse clustering.
sparcl::KMeansSparseClusteradds an \(L_1\) penalty on per-feature weights, driving uninformative features to exactly zero, the LASSO idea (Chapter 11) applied to the clustering objective. - Spectral feature selection ranks features by how well they preserve the graph structure of the data.
- Dimension reduction first. PCA or UMAP (Chapter 4) before clustering is the simplest effective option and is what most pipelines actually do.
15 PART III: HIERARCHICAL AND DENSITY-BASED METHODS
16 Hierarchical clustering
Partitioning methods commit to one value of \(k\). Hierarchical methods build a nested family of partitions for every \(k\) simultaneously, which is enormously useful when the right granularity is itself the question.
Two directions:
- Agglomerative (bottom-up). Start with \(n\) singleton clusters; repeatedly merge
the closest pair until one cluster remains. This is what
hclust()andagnes()do. - Divisive (top-down). Start with one cluster;
repeatedly split.
diana()implements this. It is \(O(2^n)\) to do optimally and so is used less often.
16.1 Linkage: how to measure distance between clusters
A distance between points does not by itself define a distance between sets. The linkage criterion supplies that, and it is a modelling choice with large consequences.
\[ \begin{aligned} \textbf{Single: }\ & d(A,B)=\min_{a\in A,\,b\in B} d(a,b) && \text{nearest neighbour}\\ \textbf{Complete: }\ & d(A,B)=\max_{a\in A,\,b\in B} d(a,b) && \text{furthest neighbour}\\ \textbf{Average (UPGMA): }\ & d(A,B)=\frac{1}{|A||B|}\sum_{a\in A}\sum_{b\in B} d(a,b) &&\\ \textbf{Centroid: }\ & d(A,B)=\lVert \bar a-\bar b\rVert^2 &&\\ \textbf{Ward: }\ & d(A,B)=\frac{|A||B|}{|A|+|B|}\lVert \bar a-\bar b\rVert^2 && \text{minimum variance increase} \end{aligned} \]
Their behaviour differs sharply. Single linkage produces chaining, a sequence of nearby points can string two well-separated clusters together, but it is the only linkage that recovers arbitrarily shaped, connected clusters. Complete linkage produces compact, roughly spherical clusters and is sensitive to outliers. Ward minimizes the increase in total within-cluster variance at each merge, so it behaves like a hierarchical \(k\)-means and tends to produce balanced, spherical groups.
ward.Dandward.D2are not the same. Ward’s criterion is defined on squared Euclidean distances.hclust(d, "ward.D2")squares the input distances internally, implementing the criterion correctly for adistobject holding ordinary Euclidean distances.ward.Ddoes not square, so it implements Ward’s method only if you pass it \(d^2\) yourself. Useward.D2with a plain Euclideandist.
16.2 The Lance–Williams recurrence
All of these linkages are special cases of a single update formula (Lance & Williams, 1967). After merging clusters \(i\) and \(j\) into \((ij)\), the distance to any remaining cluster \(k\) is
\[\boxed{\;d\big((ij),k\big)=\alpha_i\,d(i,k)+\alpha_j\,d(j,k)+\beta\,d(i,j)+\gamma\,\big|d(i,k)-d(j,k)\big|\;}\]
| Linkage | \(\alpha_i\) | \(\alpha_j\) | \(\beta\) | \(\gamma\) |
|---|---|---|---|---|
| Single | \(1/2\) | \(1/2\) | \(0\) | \(-1/2\) |
| Complete | \(1/2\) | \(1/2\) | \(0\) | \(+1/2\) |
| Average | \(\frac{n_i}{n_i+n_j}\) | \(\frac{n_j}{n_i+n_j}\) | \(0\) | \(0\) |
| Centroid | \(\frac{n_i}{n_i+n_j}\) | \(\frac{n_j}{n_i+n_j}\) | \(-\frac{n_in_j}{(n_i+n_j)^2}\) | \(0\) |
| Ward | \(\frac{n_i+n_k}{n_i+n_j+n_k}\) | \(\frac{n_j+n_k}{n_i+n_j+n_k}\) | \(-\frac{n_k}{n_i+n_j+n_k}\) | \(0\) |
Why this matters computationally. The recurrence means the algorithm never needs to revisit the original data, only the \(n\times n\) dissimilarity matrix, updated in place. That gives \(O(n^2)\) memory and \(O(n^3)\) time in the general case, reducible to \(O(n^2)\) for single and complete linkage with the SLINK and CLINK algorithms.
# Verify the recurrence reproduces hclust's merge heights for complete linkage
set.seed(251)
Xlw <- matrix(rnorm(12 * 2), 12, 2)
Dlw <- as.matrix(dist(Xlw))
lw_complete <- function(D) {
n <- nrow(D); active <- seq_len(n); heights <- numeric(0)
diag(D) <- Inf
while (length(active) > 1) {
sub <- D[active, active, drop = FALSE]
ij <- which(sub == min(sub), arr.ind = TRUE)[1, ]
i <- active[ij[1]]; j <- active[ij[2]]
heights <- c(heights, D[i, j])
# Complete linkage: alpha = 1/2, 1/2; beta = 0; gamma = +1/2 ==> max
for (k in setdiff(active, c(i, j)))
D[i, k] <- D[k, i] <- 0.5 * D[i, k] + 0.5 * D[j, k] + 0.5 * abs(D[i, k] - D[j, k])
active <- setdiff(active, j); D[j, ] <- D[, j] <- Inf
}
sort(heights)
}
c(max_abs_difference =
max(abs(lw_complete(Dlw) - sort(hclust(dist(Xlw), "complete")$height))))#> max_abs_difference
#> 4.44089e-16The hand-rolled recurrence and hclust agree to machine
precision.
16.3 Dendrograms and where to cut
library(ggdendro)
D_div <- if (DIVORCE_DISTANCE == "gower") d_gower else d_euclid
linkages <- c("single", "complete", "average", "ward.D2")
hc_list <- lapply(linkages, \(m) hclust(D_div, method = m))
names(hc_list) <- linkages
dend_df <- bind_rows(lapply(linkages, function(m) {
dd <- dendro_data(as.dendrogram(hc_list[[m]]))
transform(segment(dd), linkage = m)
}))
ggplot(dend_df) +
geom_segment(aes(x = x, y = y, xend = xend, yend = yend), linewidth = 0.35) +
facet_wrap(~ linkage, scales = "free", nrow = 2) +
labs(title = "Four linkages on the same dissimilarity matrix",
subtitle = "Single linkage chains; Ward produces balanced groups. The data has not changed",
x = NULL, y = "Merge height") +
theme_dspa(9) +
theme(axis.text.x = element_blank())
Compare linkages at a common \(k\). Silhouette varies systematically with \(k\), so cutting one dendrogram at \(k=3\) and another at \(k=10\) and comparing their silhouettes measures the effect of \(k\), not of the linkage.
K_COMMON <- 3
do.call(rbind, lapply(linkages, function(m) {
cl <- cutree(hc_list[[m]], k = K_COMMON)
s <- mean(silhouette(cl, D_div)[, "sil_width"])
data.frame(linkage = m, k = K_COMMON,
sizes = paste(table(cl), collapse = " / "),
avg_silhouette = round(s, 4),
verdict = as.character(sil_verdict(s)))
}))16.4 Cophenetic correlation
A dendrogram imposes a tree metric: the cophenetic distance between two observations is the height at which they first join. How faithfully does that tree reproduce the original dissimilarities?
\[\boxed{\;c=\mathrm{cor}\big(d_{ij},\ \mathrm{coph}_{ij}\big)\;}\]
Values above about 0.75 indicate the hierarchy is a reasonable summary; low values mean the tree has distorted the distances substantially.
coph_tab <- do.call(rbind, lapply(linkages, function(m) {
cph <- cophenetic(hc_list[[m]])
data.frame(linkage = m, cophenetic_correlation = round(cor(D_div, cph), 4))
}))
coph_tabbest_link <- coph_tab$linkage[which.max(coph_tab$cophenetic_correlation)]
ggplot(data.frame(original = as.numeric(D_div),
cophenetic = as.numeric(cophenetic(hc_list[[best_link]]))),
aes(original, cophenetic)) +
geom_point(alpha = 0.15, size = 0.8, color = "steelblue") +
geom_abline(slope = 1, intercept = 0, color = "firebrick") +
labs(title = sprintf("Cophenetic vs. original distances (%s linkage)", best_link),
subtitle = "The staircase pattern is intrinsic: a tree can only realize as many distinct distances as it has merges",
x = "Original dissimilarity", y = "Cophenetic distance") +
theme_dspa()
Average linkage usually maximizes cophenetic correlation, because it is the linkage that explicitly averages the pairwise distances it is trying to represent. That does not make it the best clustering, faithfulness to the distance matrix and usefulness of the partition are different criteria.
# # A dendrogram is a tree, and plotly renders trees natively as sunburst/treemap
library(plotly)
hc_best <- hc_list[["ward.D2"]]
merges <- hc_best$merge
hts <- hc_best$height
n_obs <- nrow(merges) + 1
node_id <- function(x) if (x < 0) paste0("obs", -x) else paste0("node", x)
edges <- do.call(rbind, lapply(seq_len(nrow(merges)), function(i) {
data.frame(
parent = node_id(i),
child = c(node_id(merges[i, 1]), node_id(merges[i, 2])),
height = hts[i],
stringsAsFactors = FALSE
)
}))
root <- setdiff(edges$parent, edges$child)
# Count leaves under every node
size <- setNames(rep(1, n_obs), paste0("obs", seq_len(n_obs)))
for (i in seq_len(nrow(merges))) {
left <- merges[i, 1]
right <- merges[i, 2]
size[[node_id(i)]] <-
(if (left < 0) 1 else size[[node_id(left)]]) +
(if (right < 0) 1 else size[[node_id(right)]])
}
# Build nodes with correct parent and value vectors
nodes <- data.frame(
id = c(root, edges$child),
parent = c("", edges$parent), # parent for each child is just edges$parent
label = c(root, edges$child),
value = size[c(root, edges$child)], # correct subtree sizes
stringsAsFactors = FALSE
)
# after building the nodes order them for sunburst ....
nodes <- nodes[order(-nodes$value), ] # root 47 -> ... -> leaves 1
# parent will always be larger than child for your size metric
plot_ly(ids = nodes$id, labels = nodes$label, parents = nodes$parent,
values = nodes$value, type = "sunburst", branchvalues = "total") |>
layout(title = "Ward dendrogram as a sunburst")17 Density-based clustering
Every method so far partitions all the data and requires \(k\) in advance. Neither is appropriate when clusters have arbitrary shapes, when their number is unknown, or when some observations are genuinely noise.
DBSCAN (Ester et al., 1996) addresses all three by defining clusters as regions of high density separated by regions of low density.
17.1 The definitions
Two parameters: a radius \(\varepsilon\) and a count \(\mathrm{minPts}\).
- The \(\varepsilon\)-neighbourhood of \(p\) is \(N_\varepsilon(p)=\{q: d(p,q)\le\varepsilon\}\).
- \(p\) is a core point if \(|N_\varepsilon(p)|\ge \mathrm{minPts}\).
- \(q\) is directly density-reachable from a core point \(p\) if \(q\in N_\varepsilon(p)\).
- Density-reachable is the transitive closure of that relation through core points.
- A cluster is a maximal set of mutually density-connected points.
- Points reachable from no core point are noise, and DBSCAN labels them 0.
Three properties follow immediately, and each fixes a \(k\)-means limitation:
Clusters may have any shape, because membership propagates through chains of core points rather than radiating from a centroid. The number of clusters is discovered, not specified. Noise is an explicit output category, not forced into a cluster.
17.2 Choosing \(\varepsilon\) and minPts
minPts sets the minimum cluster size and the noise
tolerance. The usual heuristic is \(\mathrm{minPts}\ge d+1\), commonly \(2d\); larger values give more robust
results on noisy data.
\(\varepsilon\) is read off a \(k\)-distance plot: sort every point’s distance to its \(\mathrm{minPts}\)-th nearest neighbour and look for the knee. Below the knee, points are in dense regions; above it, they are in sparse ones.
library(dbscan)
# Use the nested-rings data, where k-means failed completely (Section 8.8)
minPts <- 5
kd <- sort(dbscan::kNNdist(Xb, k = minPts))
kd_df <- data.frame(index = seq_along(kd), distance = kd)
# The same kink-detection geometry used for the elbow in Section 8.7
kn <- elbow_point(kd_df$index, kd_df$distance)
eps_hat <- kd_df$distance[kd_df$index == kn$k]
c(minPts = minPts, suggested_eps = round(eps_hat, 4))#> minPts suggested_eps
#> 5.0000 0.1957ggplot(kd_df, aes(index, distance)) +
geom_line(linewidth = 0.9, color = "steelblue") +
geom_hline(yintercept = eps_hat, linetype = "dashed", color = "firebrick") +
ggplot2::annotate("text", x = 50, y = eps_hat * 1.25, hjust = 0, size = 3.2,
color = "firebrick", label = sprintf("eps = %.3f", eps_hat)) +
labs(title = sprintf("k-distance plot (k = minPts = %d)", minPts),
subtitle = "The knee separates points in dense regions from points in sparse ones",
x = "Points, sorted by distance to their k-th neighbour",
y = sprintf("Distance to %d-th nearest neighbour", minPts)) +
theme_dspa()
run_dbscan <- function(X, y, label, minPts = 5) {
kd <- sort(dbscan::kNNdist(X, k = minPts))
eps <- kd[elbow_point(seq_along(kd), kd)$k]
db <- dbscan::dbscan(X, eps = eps, minPts = minPts)
km <- kmeans(X, length(setdiff(unique(db$cluster), 0)), nstart = 25)
list(db = db, eps = eps,
row = data.frame(case = label, eps = round(eps, 3),
dbscan_clusters = length(setdiff(unique(db$cluster), 0)),
noise_points = sum(db$cluster == 0),
dbscan_ARI = round(mclust::adjustedRandIndex(db$cluster, y), 4),
kmeans_ARI = round(mclust::adjustedRandIndex(km$cluster, y), 4)))
}
set.seed(261)
db_res <- lapply(names(cases), function(nm)
run_dbscan(cases[[nm]][[1]], cases[[nm]][[2]], nm))
do.call(rbind, lapply(db_res, \(r) r$row))db_df <- bind_rows(lapply(seq_along(cases), function(i) {
X <- cases[[i]][[1]]
data.frame(x = X[, 1], y = X[, 2],
cl = factor(db_res[[i]]$db$cluster),
panel = names(cases)[i])
}))
p_db <- ggplot(db_df, aes(x, y, color = cl)) +
geom_point(size = 0.6, alpha = 0.75) +
facet_wrap(~ panel, nrow = 1, scales = "free") +
scale_color_manual(values = c("0" = "grey75",
setNames(RColorBrewer::brewer.pal(8, "Set1"),
as.character(1:8))),
guide = "none") +
labs(title = "DBSCAN on the same four structures",
subtitle = "Grey points are labelled NOISE -- a category k-means does not have",
x = NULL, y = NULL) + theme_dspa(8)
p_km / p_db
Compare the ARI columns. On the nested rings DBSCAN recovers the structure essentially perfectly while \(k\)-means scores near zero, because density connectivity follows the ring and Euclidean distance to a centroid does not. On unequal density DBSCAN struggles for the opposite reason: a single global \(\varepsilon\) cannot be right for both a tight cluster and a diffuse one.
Common misconception: “DBSCAN has no parameters to tune, so it is objective.” It has two, and \(\varepsilon\) is at least as consequential as \(k\). The difference is that \(\varepsilon\) is a density threshold rather than a count, so it can be estimated from the data’s own \(k\)-distance profile, but a single global \(\varepsilon\) still assumes uniform density across clusters, which is exactly the assumption HDBSCAN removes.
17.3 HDBSCAN
Hierarchical Density-Based Spatial Clustering of Applications
with Noise (HDBSCAN) (Campello et al.,
2013) runs DBSCAN across all values of \(\varepsilon\) simultaneously, builds the
resulting hierarchy, and extracts the clusters that persist longest, a
stability criterion over the density scale. It requires only
minPts and handles clusters of differing density.
set.seed(271)
hdb_res <- do.call(rbind, lapply(names(cases), function(nm) {
X <- cases[[nm]][[1]]; y <- cases[[nm]][[2]]
h <- dbscan::hdbscan(X, minPts = 8)
data.frame(case = nm,
hdbscan_clusters = length(setdiff(unique(h$cluster), 0)),
noise = sum(h$cluster == 0),
hdbscan_ARI = round(mclust::adjustedRandIndex(h$cluster, y), 4))
}))
hdb_resh_unequal <- dbscan::hdbscan(cases[["3. Unequal density and size"]][[1]], minPts = 8)
db_unequal <- db_res[[3]]$db
Xu <- cases[["3. Unequal density and size"]][[1]]
bind_rows(
data.frame(x = Xu[, 1], y = Xu[, 2], cl = factor(db_unequal$cluster),
method = "DBSCAN (single global eps)"),
data.frame(x = Xu[, 1], y = Xu[, 2], cl = factor(h_unequal$cluster),
method = "HDBSCAN (density hierarchy)")) |>
ggplot(aes(x, y, color = cl)) +
geom_point(size = 0.8, alpha = 0.8) +
facet_wrap(~ method) +
scale_color_manual(values = c("0" = "grey75", "1" = "#3B7DD8",
"2" = "#D8433B", "3" = "#7FB069"),
guide = "none") +
labs(title = "Where a single density threshold fails",
subtitle = "Grey is noise. HDBSCAN adapts to clusters of different density; DBSCAN cannot",
x = NULL, y = NULL) +
theme_dspa(10)
# --- Interactive equivalent, with the condensed cluster tree ---------------
plot(h_unequal, show_flat = TRUE) # base-graphics condensed tree
plot_ly(x = Xu[, 1], y = Xu[, 2], type = "scatter", mode = "markers",
color = factor(h_unequal$cluster),
text = ~paste("membership:", round(h_unequal$membership_prob, 2)),
marker = list(size = 6, opacity = 0.8)) |>
layout(title = "HDBSCAN clusters with membership probabilities",
xaxis = list(title = "x"), yaxis = list(title = "y"))HDBSCAN also returns a membership probability per point, so the boundary cases are identifiable rather than silently assigned, the same soft-assignment idea that Gaussian mixtures provide parametrically (§8.16).
# One table, four structures, three methods
do.call(rbind, lapply(seq_along(cases), function(i) {
X <- cases[[i]][[1]]; y <- cases[[i]][[2]]
km <- kmeans(X, cases[[i]][[3]], nstart = 25)
h <- dbscan::hdbscan(X, minPts = 8)
data.frame(structure = names(cases)[i],
kmeans = round(mclust::adjustedRandIndex(km$cluster, y), 3),
dbscan = db_res[[i]]$row$dbscan_ARI,
hdbscan = round(mclust::adjustedRandIndex(h$cluster, y), 3))
}))No method wins everywhere, and the pattern of wins is predictable from the assumptions each makes. That is the practical content of this section: choose the method from the geometry you expect, and when you do not know the geometry, run more than one and compare.
18 PART IV: MODEL-BASED AND SPECTRAL METHODS
19 Gaussian mixture models
Every method so far assigns each point to exactly one cluster. A mixture model instead treats the data as generated by a probabilistic process, which buys three things: a principled way to choose \(k\), soft assignments with uncertainty attached, and cluster shapes controlled by explicit parameters rather than implied by a distance.
19.1 The model
\[\boxed{\;p(\mathbf{x}\mid\theta)=\sum_{\ell=1}^{k}\pi_\ell\,\mathcal{N}\big(\mathbf{x}\mid\boldsymbol\mu_\ell,\Sigma_\ell\big), \qquad \pi_\ell\ge0,\ \sum_\ell\pi_\ell=1\;}\]
with parameters \(\theta=\{\pi_\ell,\boldsymbol\mu_\ell,\Sigma_\ell\}_{\ell=1}^{k}\). Introduce a latent indicator \(z_i\in\{1,\dots,k\}\) giving the component that generated \(\mathbf{x}_i\).
If the \(z_i\) were observed, estimation would be trivial: the complete-data log-likelihood separates,
\[\ell_c(\theta)=\sum_{i=1}^{n}\sum_{\ell=1}^{k} z_{i\ell}\Big[\log\pi_\ell+\log\mathcal{N}(\mathbf{x}_i\mid\boldsymbol\mu_\ell,\Sigma_\ell)\Big],\]
and the two blocks, the \(\pi_\ell\) and the \((\boldsymbol\mu_\ell,\Sigma_\ell)\), are decoupled, each with a closed-form maximizer.
The \(z_i\) are not observed, so the observed-data log-likelihood \(\ell(\theta)=\sum_i\log\sum_\ell \pi_\ell\mathcal{N}(\mathbf{x}_i\mid\cdot)\) has a sum inside the logarithm and no closed-form maximum.
19.2 EM
The expectation–maximization algorithm alternates between filling in the expected \(z_i\) and maximizing as though they were observed.
E step. Compute the responsibility, the posterior probability that component \(\ell\) generated observation \(i\):
\[\boxed{\;\gamma_{i\ell}=\frac{\pi_\ell\,\mathcal{N}(\mathbf{x}_i\mid\boldsymbol\mu_\ell,\Sigma_\ell)}{\sum_{m=1}^{k}\pi_m\,\mathcal{N}(\mathbf{x}_i\mid\boldsymbol\mu_m,\Sigma_m)}\;}\]
M step. With \(N_\ell=\sum_i\gamma_{i\ell}\),
\[\pi_\ell\leftarrow\frac{N_\ell}{n},\qquad \boldsymbol\mu_\ell\leftarrow\frac{1}{N_\ell}\sum_i\gamma_{i\ell}\mathbf{x}_i,\qquad \Sigma_\ell\leftarrow\frac{1}{N_\ell}\sum_i\gamma_{i\ell}(\mathbf{x}_i-\boldsymbol\mu_\ell)(\mathbf{x}_i-\boldsymbol\mu_\ell)^\top .\]
These are exactly the sample mean and covariance, weighted by responsibility. EM increases the observed-data likelihood at every iteration (never decreases it), and, like Lloyd’s algorithm, converges to a local optimum that depends on initialization.
\(k\)-means is a limiting case of GMM-EM. Fix \(\Sigma_\ell=\sigma^2 I\) for all \(\ell\) and let \(\sigma^2\to0\). The responsibilities become hard 0/1 indicators (the nearest centre gets all the mass), the M step reduces to “recompute the mean of each hard-assigned group,” and the algorithm becomes Lloyd’s. Everything \(k\)-means assumes, spherical, equal-volume clusters, is visible as a restriction on \(\Sigma_\ell\), not a separate list of caveats.
em_gmm <- function(X, k, iters = 60, seed = 1) {
set.seed(seed); X <- as.matrix(X); n <- nrow(X); d <- ncol(X)
pi_l <- rep(1 / k, k)
mu <- X[sample(n, k), , drop = FALSE]
Sig <- replicate(k, cov(X), simplify = FALSE)
ll <- numeric(iters)
for (t in seq_len(iters)) {
# E step: responsibilities
dens <- vapply(seq_len(k), \(l)
pi_l[l] * mvtnorm::dmvnorm(X, mu[l, ], Sig[[l]]), numeric(n))
tot <- rowSums(dens); ll[t] <- sum(log(tot))
G <- dens / tot
# M step: responsibility-weighted moments
N_l <- colSums(G)
pi_l <- N_l / n
mu <- t(vapply(seq_len(k), \(l) colSums(G[, l] * X) / N_l[l], numeric(d)))
Sig <- lapply(seq_len(k), function(l) {
C <- sweep(X, 2, mu[l, ]); crossprod(C * sqrt(G[, l])) / N_l[l]
})
}
list(pi = pi_l, mu = mu, Sigma = Sig, resp = G, loglik = ll,
cluster = max.col(G))
}
set.seed(281)
fit_em <- em_gmm(Xg, 3, seed = 281)
c(monotone_increasing = all(diff(fit_em$loglik) >= -1e-8),
final_loglik = round(max(fit_em$loglik), 3),
ARI_vs_truth = round(mclust::adjustedRandIndex(fit_em$cluster, g), 4))#> monotone_increasing final_loglik ARI_vs_truth
#> 1.0000 -1341.3800 0.3412ggplot(data.frame(iter = seq_along(fit_em$loglik), ll = fit_em$loglik),
aes(iter, ll)) +
geom_line(linewidth = 1, color = "steelblue") + geom_point(size = 1.4) +
labs(title = "EM increases the log-likelihood monotonically",
subtitle = "Guaranteed by the algorithm; the limit is a local optimum, not necessarily the global one",
x = "Iteration", y = "Observed-data log-likelihood") +
theme_dspa()
19.3 Covariance parameterizations
The shape, volume, and orientation of each component are all encoded in \(\Sigma_\ell\). Its eigendecomposition \(\Sigma_\ell=\lambda_\ell D_\ell A_\ell D_\ell^\top\) separates them: \(\lambda_\ell\) is volume, \(A_\ell\) is shape, \(D_\ell\) is orientation. Constraining each to be equal (E) or varying (V) across components gives 14 models:
| Model | Volume | Shape | Orientation | Geometry |
|---|---|---|---|---|
| EII | Equal | Equal | — | Spherical, equal size — this is \(k\)-means |
| VII | Varying | Equal | — | Spherical, different sizes |
| EEI | Equal | Equal | Axis-aligned | Diagonal |
| VEI | Varying | Equal | Axis-aligned | Diagonal |
| EVI | Equal | Varying | Axis-aligned | Diagonal |
| VVI | Varying | Varying | Axis-aligned | Diagonal |
| EEE | Equal | Equal | Equal | Ellipsoidal, identical |
| EVE, VEE, VVE | mixed | mixed | Equal | Ellipsoidal, shared orientation |
| EEV, VEV | mixed | Equal | Varying | Ellipsoidal |
| EVV | Equal | Varying | Varying | Ellipsoidal |
| VVV | Varying | Varying | Varying | Fully general |
More constrained models have fewer parameters and are estimable from
less data; VVV needs roughly \(n_\ell\gg d^2/2\) per component.
BIC selects both the model and \(k\) simultaneously:
\[\mathrm{BIC}=2\ell(\hat\theta)-\nu\log n,\]
with \(\nu\) the parameter count.
mclust reports it with larger is better
(note the sign convention differs from stats::BIC).
set.seed(291)
gmm <- Mclust(Xg)
c(selected_model = gmm$modelName, components = gmm$G,
BIC = round(gmm$bic, 2),
ARI_vs_truth = round(mclust::adjustedRandIndex(gmm$classification, g), 4))#> selected_model components BIC ARI_vs_truth
#> "EII" "3" "-2745.7" "0.3684"bic_df <- as.data.frame(as.table(gmm$BIC)) |>
setNames(c("k", "model", "BIC")) |>
filter(!is.na(BIC)) |> mutate(k = as.integer(as.character(k)))
ggplot(bic_df, aes(k, BIC, color = model, group = model)) +
geom_line(linewidth = 0.6) + geom_point(size = 1.2) +
geom_point(data = bic_df[which.max(bic_df$BIC), ], color = "black",
size = 4, shape = 21, stroke = 1.2, fill = NA) +
labs(title = "BIC across 14 covariance models and every k",
subtitle = sprintf("Circled: the maximum, at %s with %d components. Larger BIC is better in mclust's convention",
gmm$modelName, gmm$G),
x = "Number of components k", y = "BIC", color = "Model") +
theme_dspa(9)
19.4 Soft assignment
The responsibilities \(\gamma_{i\ell}\) are the model’s own statement of how confident it is. Points with a maximum responsibility near 1 sit clearly inside a component; points near \(1/k\) are genuinely ambiguous.
unc <- data.frame(x = Xg[, 1], y = Xg[, 2],
cluster = factor(gmm$classification),
max_resp = apply(gmm$z, 1, max))
c(confidently_assigned = sum(unc$max_resp > 0.95),
ambiguous = sum(unc$max_resp < 0.75),
total = nrow(unc))#> confidently_assigned ambiguous total
#> 283 24 360ggplot(unc, aes(x, y, color = cluster, alpha = max_resp)) +
geom_point(size = 1.8) +
scale_color_brewer(palette = "Set1", name = "Component") +
scale_alpha_continuous(range = c(0.15, 1), name = "Max responsibility") +
labs(title = "Soft assignment: faint points are the ambiguous ones",
subtitle = "k-means has no way to express this -- every point is assigned with equal confidence",
x = expression(x[1]), y = expression(x[2])) +
theme_dspa()
The mixture density is a surface over the feature plane, and comparing it with the \(k\)-means Voronoi partition of the same data makes the difference concrete:
gx2 <- seq(min(Xg[, 1]) - 1, max(Xg[, 1]) + 1, length.out = 80)
gy2 <- seq(min(Xg[, 2]) - 1, max(Xg[, 2]) + 1, length.out = 80)
grid2 <- as.matrix(expand.grid(x = gx2, y = gy2))
dens_gmm <- matrix(predict(densityMclust(Xg, plot = FALSE), grid2), length(gx2))
km3 <- kmeans(Xg, 3, nstart = 25)
D_km <- as.matrix(dist(rbind(km3$centers, grid2)))[-(1:3), 1:3]
vor <- matrix(max.col(-D_km), length(gx2))
plot_ly() |>
add_surface(x = gy2, y = gx2, z = dens_gmm / max(dens_gmm),
showscale = FALSE, opacity = 0.95, colorscale = "Viridis",
name = "GMM density") |>
add_surface(x = gy2, y = gx2, z = (vor - 1) / 2 + 1.4,
showscale = FALSE, opacity = 0.85, colorscale = "Portland",
name = "k-means Voronoi") |>
layout(title = "Mixture density (bottom) and k-means Voronoi partition (top)",
scene = list(xaxis = list(title = "x2"), yaxis = list(title = "x1"),
zaxis = list(title = "Density / partition + offset")))The mixture is a smooth density whose level sets are ellipses adapted to each component; the Voronoi diagram is a piecewise-constant partition with straight boundaries equidistant between centres. The second is what you get when you force \(\Sigma_\ell=\sigma^2 I\) and \(\sigma^2\to0\).
19.5 Univariate mixtures
The same machinery applies in one dimension, where it is used for deconvolution: separating overlapping populations in a single measured quantity.
library(mixtools)
crystal <- dspa_read(
"https://umich.instructure.com/files/11653615/download?download_frd=1",
"crystallography.csv", header = TRUE)
c(columns = ncol(crystal), rows = nrow(crystal))#> columns rows
#> 9 653# The number of components is SELECTED per column by BIC, not fixed at 3
fit_column <- function(v, kmax = 5) {
v <- v[complete.cases(v)]
bics <- vapply(1:kmax, function(k) {
if (k == 1) return(2 * sum(dnorm(v, mean(v), sd(v), log = TRUE)) - 2 * log(length(v)))
f <- try(suppressMessages(normalmixEM(v, k = k, maxit = 500, epsilon = 1e-6)),
silent = TRUE)
if (inherits(f, "try-error")) return(-Inf)
2 * f$loglik - (3 * k - 1) * log(length(v)) # mclust sign convention
}, numeric(1))
k_best <- which.max(bics)
fit <- if (k_best == 1) NULL else
suppressMessages(normalmixEM(v, k = k_best, maxit = 500, epsilon = 1e-6))
list(k = k_best, bic = bics, fit = fit, data = v)
}
set.seed(301)
col1 <- fit_column(crystal[[1]])#> number of iterations= 159
#> WARNING! NOT CONVERGENT!
#> number of iterations= 500
#> number of iterations= 132
#> number of iterations= 256
#> number of iterations= 183c(column = names(crystal)[1], selected_k = col1$k,
BIC_by_k = paste(round(col1$bic, 1), collapse = ", "))#> column
#> "AC1338"
#> selected_k
#> "3"
#> BIC_by_k
#> "-4694.3, -4647.9, -4642.1, -4647.7, -4659.5"mix_density <- function(fit, xs) {
if (is.null(fit)) return(NULL)
k <- length(fit$lambda)
comps <- lapply(seq_len(k), \(i)
fit$lambda[i] * dnorm(xs, fit$mu[i], fit$sigma[i]))
list(components = comps, total = Reduce(`+`, comps))
}
xs <- seq(min(col1$data), max(col1$data), length.out = 400)
md <- mix_density(col1$fit, xs)
comp_df <- if (!is.null(md)) bind_rows(lapply(seq_along(md$components), \(i)
data.frame(x = xs, y = md$components[[i]],
component = paste("Component", i)))) else NULL
p <- ggplot(data.frame(v = col1$data), aes(v)) +
geom_histogram(aes(y = after_stat(density)), bins = 30,
fill = "grey85", color = "white")
if (!is.null(md)) {
p <- p +
geom_line(data = comp_df, aes(x, y, color = component), linewidth = 0.8) +
geom_line(data = data.frame(x = xs, y = md$total), aes(x, y),
color = "black", linewidth = 1.1)
}
p + scale_color_brewer(palette = "Set1") +
labs(title = sprintf("Gaussian mixture on %s, k selected by BIC", names(crystal)[1]),
subtitle = sprintf("BIC chose k = %d. Black: the mixture; colored: its components",
col1$k),
x = "Intensity", y = "Density", color = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
pl <- plot_ly() |>
add_histogram(x = ~col1$data, histnorm = "probability density",
name = "Observed", opacity = 0.6)
for (i in seq_along(md$components))
pl <- add_lines(pl, x = xs, y = md$components[[i]],
name = paste("Component", i))
pl |> add_lines(x = xs, y = md$total, name = "Mixture",
line = list(color = "black", width = 4)) |>
layout(title = "Normal mixture model", barmode = "overlay",
xaxis = list(title = "Intensity"), yaxis = list(title = "Density"))20 Spectral clustering
Spectral clustering converts the data into a similarity graph, embeds that graph in a low-dimensional space using the eigenvectors of its Laplacian, and clusters there. Its appeal is that connectivity, not compactness, drives the partition, so it handles the shapes that defeat \(k\)-means.
20.1 The similarity graph
Build a weighted graph with the observations as vertices. Three standard constructions:
- \(\varepsilon\)-neighbourhood: connect points closer than \(\varepsilon\).
- \(k\)-nearest-neighbour: connect \(i\) to its \(k\) nearest (symmetrized).
- Fully connected with a Gaussian kernel: \[w_{ij}=\exp\!\left(-\frac{\lVert x_i-x_j\rVert^2}{2\sigma^2}\right),\] where \(\sigma\) sets the neighbourhood width.
Let \(W\) be the weighted adjacency matrix and \(D=\operatorname{diag}(d_1,\dots,d_n)\) with \(d_i=\sum_j w_{ij}\) the degree matrix.
20.2 The three Laplacians
There are three graph Laplacians, and they are not interchangeable.
\[ \begin{aligned} \textbf{Unnormalized: }\quad & L = D - W\\[1mm] \textbf{Symmetric normalized: }\quad & L_{\mathrm{sym}} = D^{-1/2}LD^{-1/2} = I - D^{-1/2}WD^{-1/2}\\[1mm] \textbf{Random-walk: }\quad & L_{\mathrm{rw}} = D^{-1}L = I - D^{-1}W \end{aligned} \]
Their key properties:
- \(L\) and \(L_{\mathrm{sym}}\) are symmetric and positive semi-definite. \(L_{\mathrm{rw}}\) is not symmetric unless all degrees are equal.
- Their eigenvalues are related: \(\lambda\) is an eigenvalue of \(L_{\mathrm{rw}}\) with eigenvector \(u\) iff it is an eigenvalue of \(L_{\mathrm{sym}}\) with eigenvector \(D^{1/2}u\).
- For every Laplacian, \(\lambda_1=0\), and its eigenvector for \(L\) is the constant vector \(\mathbf{1}\).
This matters for the code, not just the theory.
eigen(M, symmetric = TRUE)instructs LAPACK to read only the lower triangle ofMand assume symmetry. Passing \(L_{\mathrm{rw}}=I-D^{-1}W\), which is not symmetric — therefore decomposes a different matrix than the one supplied, silently. Use \(L_{\mathrm{sym}}\) when you want a symmetric solver, and recover \(L_{\mathrm{rw}}\)’s eigenvectors as \(D^{-1/2}u_{\mathrm{sym}}\) if you need them.
set.seed(311)
Wt <- matrix(runif(36), 6, 6); Wt <- (Wt + t(Wt)) / 2; diag(Wt) <- 0
Dt <- diag(rowSums(Wt))
L_un <- Dt - Wt
L_sym <- diag(6) - diag(1 / sqrt(diag(Dt))) %*% Wt %*% diag(1 / sqrt(diag(Dt)))
L_rw <- diag(6) - solve(Dt) %*% Wt
c(L_symmetric = isTRUE(all.equal(L_un, t(L_un))),
Lsym_symmetric = isTRUE(all.equal(L_sym, t(L_sym))),
Lrw_symmetric = isTRUE(all.equal(L_rw, t(L_rw), tolerance = 1e-12)),
max_asymmetry_of_Lrw = round(max(abs(L_rw - t(L_rw))), 4))#> L_symmetric Lsym_symmetric Lrw_symmetric
#> 1.0000 1.0000 0.0000
#> max_asymmetry_of_Lrw
#> 0.0349e_sym <- eigen(L_sym, symmetric = TRUE)
e_rw <- eigen(L_rw) # NOT symmetric = TRUE
c(eigenvalues_agree = max(abs(sort(Re(e_rw$values)) - sort(e_sym$values))))#> eigenvalues_agree
#> 2.66454e-15# u_rw = D^{-1/2} u_sym
u_sym <- e_sym$vectors[, ncol(e_sym$vectors)]
u_rw_derived <- diag(1 / sqrt(diag(Dt))) %*% u_sym
u_rw_derived <- u_rw_derived / sqrt(sum(u_rw_derived^2))
c(constant_eigenvector_recovered =
round(sd(u_rw_derived) / abs(mean(u_rw_derived)), 6))#> constant_eigenvector_recovered
#> 0The eigenvalues coincide and the transformation recovers the constant eigenvector, confirming the relation.
20.3 Why the smallest eigenvalue is discarded
For a connected graph, \(L\mathbf{1}=0\): every row of \(L\) sums to zero, so the constant vector is an eigenvector with eigenvalue exactly 0. It assigns the same value to every vertex and therefore carries no partition information.
More generally, the multiplicity of \(\lambda=0\) equals the number of connected components of the graph. If the graph splits into \(c\) pieces, the zero-eigenspace is spanned by the indicator vectors of those pieces, and reading them off gives the clustering directly.
The second-smallest eigenvalue \(\lambda_2\) is the algebraic connectivity, and its eigenvector is the Fiedler vector. Thresholding it at zero (or at its median) gives the classical two-way spectral cut.
Common misconception: “take the \(k\) smallest eigenvectors.” Take eigenvectors \(2,\dots,k+1\), the \(k\) smallest excluding \(\lambda_1=0\). Including the trivial constant vector adds a zero-variance column to the embedding.
The reason this error survives is that it is invisible in the case people test. At \(k=2\) a constant column has no effect on \(k\)-means, so the partition is unchanged and everything looks fine. At \(k\ge3\) the trivial vector displaces a genuinely informative eigenvector and the partition degrades, silently, with no warning and no error. Problem 7 measures the gap.
A related trap:
eigen()sorts eigenvalues descending, so the smallest are at the end of the returned matrix. Indexing the first \(k\) columns gives the largest eigenvalues, which is the opposite of what spectral clustering needs.
build_W <- function(X, sigma) {
D2 <- as.matrix(dist(X))^2
W <- exp(-D2 / (2 * sigma^2)); diag(W) <- 0; W
}
set.seed(321)
Xr_small <- Xb[sample(nrow(Xb), 300), ] # nested rings, subsampled
yr_small <- rep(1:2, each = nrow(Xb) / 2)[sample(nrow(Xb), 300)]
W_r <- build_W(Xr_small, sigma = 0.6)
d_r <- rowSums(W_r)
Lsym_r <- diag(nrow(W_r)) - diag(1/sqrt(d_r)) %*% W_r %*% diag(1/sqrt(d_r))
e_r <- eigen(Lsym_r, symmetric = TRUE)
lam <- rev(e_r$values) # ascending
vecs <- e_r$vectors[, rev(seq_len(ncol(e_r$vectors)))]
c(lambda_1 = round(lam[1], 8), lambda_2 = round(lam[2], 6),
lambda_3 = round(lam[3], 6),
first_eigenvector_is_constant = round(sd(vecs[, 1]) / abs(mean(vecs[, 1])), 4))#> lambda_1 lambda_2
#> 0.000000 0.012480
#> lambda_3 first_eigenvector_is_constant
#> 0.029367 0.380400data.frame(index = 1:12, eigenvalue = lam[1:12]) |>
ggplot(aes(index, eigenvalue)) +
geom_col(fill = "steelblue") +
geom_col(data = data.frame(index = 1, eigenvalue = lam[1]), fill = "firebrick") +
scale_x_continuous(breaks = 1:12) +
labs(title = "The Laplacian spectrum, smallest eigenvalues first",
subtitle = "Red bar is lambda_1 = 0 with the constant eigenvector -- discarded. The eigengap after lambda_2 suggests k = 2",
x = "Eigenvalue index", y = expression(lambda)) +
theme_dspa()
The eigengap heuristic reads \(k\) off this plot: choose \(k\) so that \(\lambda_1,\dots,\lambda_k\) are small and \(\lambda_{k+1}\) is noticeably larger.
20.4 The algorithm
- Build the similarity matrix \(W\) (sparse, if the construction allows).
- Compute \(D\) and \(L_{\mathrm{sym}}=I-D^{-1/2}WD^{-1/2}\).
- Find the eigenvectors for the \(k\) smallest eigenvalues, excluding \(\lambda_1=0\): take \(u_2,\dots,u_{k+1}\) and stack as columns of \(U\in\mathbb{R}^{n\times k}\).
- Normalize the rows of \(U\) to unit length (this is the Ng–Jordan–Weiss step, specific to \(L_{\mathrm{sym}}\)).
- Run \(k\)-means on the rows of \(U\).
spectral_cluster <- function(X, k, sigma = 1, nstart = 25) {
W <- build_W(X, sigma)
d <- rowSums(W); d[d == 0] <- 1e-12
Dm <- diag(1 / sqrt(d))
Lsym <- diag(nrow(W)) - Dm %*% W %*% Dm
e <- eigen(Lsym, symmetric = TRUE)
# eigen() sorts DESCENDING, so the smallest are last. Drop the very last
# (lambda_1 = 0, constant vector) and take the next k.
idx <- (ncol(e$vectors) - k):(ncol(e$vectors) - 1)
U <- e$vectors[, idx, drop = FALSE]
U <- U / pmax(sqrt(rowSums(U^2)), 1e-12) # row-normalize (NJW)
list(cluster = kmeans(U, k, nstart = nstart)$cluster,
embedding = U, eigenvalues = rev(e$values))
}
set.seed(331)
sp_rings <- spectral_cluster(Xr_small, k = 2, sigma = 0.6)
km_rings <- kmeans(Xr_small, 2, nstart = 25)
c(spectral_ARI = round(mclust::adjustedRandIndex(sp_rings$cluster, yr_small), 4),
kmeans_ARI = round(mclust::adjustedRandIndex(km_rings$cluster, yr_small), 4))#> spectral_ARI kmeans_ARI
#> -0.0026 -0.0025bind_rows(
data.frame(x = Xr_small[, 1], y = Xr_small[, 2],
cl = factor(km_rings$cluster), method = "k-means"),
data.frame(x = Xr_small[, 1], y = Xr_small[, 2],
cl = factor(sp_rings$cluster), method = "Spectral")) |>
ggplot(aes(x, y, color = cl)) +
geom_point(size = 1.3, alpha = 0.85) +
facet_wrap(~ method) + coord_fixed() +
scale_color_brewer(palette = "Set1", guide = "none") +
labs(title = "Concentric rings: connectivity beats compactness",
subtitle = "Spectral clustering follows the graph; k-means cuts across it",
x = NULL, y = NULL) +
theme_dspa(10)
The reason is visible in the embedding, which is genuinely three-dimensional:
sp3 <- spectral_cluster(Xr_small, k = 3, sigma = 0.6)
plot_ly(x = sp3$embedding[, 1], y = sp3$embedding[, 2], z = sp3$embedding[, 3],
type = "scatter3d", mode = "markers",
color = factor(yr_small), colors = c("#3B7DD8", "#D8433B"),
marker = list(size = 3, opacity = 0.85)) |>
layout(title = "The spectral embedding (eigenvectors 2, 3, 4 of L_sym)",
scene = list(xaxis = list(title = "u2 (Fiedler)"),
yaxis = list(title = "u3"), zaxis = list(title = "u4")))Rotate this. In the original space the two rings are perfectly interleaved and no linear boundary separates them. In the embedded space they occupy distinct regions, and ordinary \(k\)-means separates them trivially. That transformation is what spectral clustering buys, and it is the same idea as the kernel trick (Chapter 6, §6.16), applied to a graph.
20.5 Complexity
| Step | Time | Memory |
|---|---|---|
| Build dense \(W\) | \(O(n^2 d)\) | \(\mathbf{O(n^2)}\) |
| Build sparse \(W\) (\(k\)-NN or \(\varepsilon\)) | \(O(n\log n\cdot d)\) | \(O(nk)\) |
Full eigen() |
\(\mathbf{O(n^3)}\) | \(O(n^2)\) |
Sparse eigs_sym() for \(k\) smallest |
\(O(n k^2)\) per iteration | \(O(nk)\) |
| \(k\)-means on the embedding | \(O(nk^2 i)\) | \(O(nk)\) |
Dense construction plus a full eigendecomposition is the
bottleneck, and it is almost always avoidable. If the
similarity graph is local, a \(k\)-NN
graph, or a Gaussian kernel truncated at a radius, then \(W\) is sparse by
construction, and only the few smallest eigenpairs are needed.
Storing \(W\) densely and calling
eigen() computes all \(n\)
eigenpairs at \(O(n^3)\) when \(O(nk^2)\) would do.
21 Image segmentation
Spectral clustering applied to a brain hematoma MRI. Pixels are vertices; edges connect spatial neighbours with weights based on intensity similarity.
library(jpeg); library(Matrix); library(RSpectra)
img_path <- local({
p <- file.path(dspa_cache_dir(), "MRI_ImageHematoma.jpg")
if (!file.exists(p))
download.file("https://umich.instructure.com/files/1627149/download?download_frd=1",
p, mode = "wb", quiet = TRUE)
p
})
img_raw <- readJPEG(img_path)
img_full <- t(apply(img_raw[, , 1], 2, rev)) # first channel, anatomically oriented
c(original_dimensions = paste(dim(img_full), collapse = " x "))#> original_dimensions
#> "256 x 256"# Downsample by strided subsampling -- no extra image package required
newdim <- c(64, 64)
ridx <- round(seq(1, nrow(img_full), length.out = newdim[1]))
cidx <- round(seq(1, ncol(img_full), length.out = newdim[2]))
img <- img_full[ridx, cidx]
c(downsampled = paste(dim(img), collapse = " x "), pixels = prod(dim(img)))#> downsampled pixels
#> "64 x 64" "4096"plot_ly(z = ~img, type = "surface", colorscale = "Greys", reversescale = TRUE) |>
layout(title = "MRI slice as an intensity surface",
scene = list(xaxis = list(title = "column"),
yaxis = list(title = "row"),
zaxis = list(title = "Intensity")))21.1 Building the graph sparsely
# Coordinates via expand.grid, not nested loops
coords <- expand.grid(r = seq_len(nrow(img)), c = seq_len(ncol(img)))
coords$value <- as.vector(img)
n_pix <- nrow(coords)
pixdiff <- 2 # spatial radius: connect pixels within a (2*2+1)^2 window
sigma2 <- 0.01 # intensity bandwidth
# Enumerate only the pairs that can be non-zero. Each pixel has at most
# (2*pixdiff+1)^2 neighbours, so the edge list is O(n), not O(n^2).
offsets <- expand.grid(dr = -pixdiff:pixdiff, dc = -pixdiff:pixdiff)
offsets <- offsets[!(offsets$dr == 0 & offsets$dc == 0), ]
edges <- do.call(rbind, lapply(seq_len(nrow(offsets)), function(o) {
r2 <- coords$r + offsets$dr[o]; c2 <- coords$c + offsets$dc[o]
ok <- r2 >= 1 & r2 <= nrow(img) & c2 >= 1 & c2 <= ncol(img)
j <- (c2[ok] - 1) * nrow(img) + r2[ok]
data.frame(i = which(ok), j = j,
w = exp(-(coords$value[ok] - coords$value[j])^2 / sigma2))
}))
W <- sparseMatrix(i = edges$i, j = edges$j, x = edges$w,
dims = c(n_pix, n_pix))
W <- (W + Matrix::t(W)) / 2 # enforce symmetry
c(pixels = n_pix,
non_zero_entries = length(W@x),
density = signif(length(W@x) / n_pix^2, 3),
sparse_MB = round(as.numeric(object.size(W)) / 1e6, 2),
dense_MB_would_be = round(8 * n_pix^2 / 1e6, 1))#> pixels non_zero_entries density sparse_MB
#> 4.096e+03 9.450e+04 5.630e-03 1.150e+00
#> dense_MB_would_be
#> 1.342e+02The similarity matrix is under 1% dense, so the sparse representation is two orders of magnitude smaller. That is not an optimization detail, it is what makes the eigendecomposition tractable.
deg <- Matrix::rowSums(W); deg[deg == 0] <- 1e-12
Dinv_sqrt <- Diagonal(x = 1 / sqrt(deg))
Lsym <- Diagonal(n_pix) - Dinv_sqrt %*% W %*% Dinv_sqrt
# Only the few SMALLEST eigenpairs are needed. eigs_sym with which = "SM"
# uses a sparse iterative method: O(n k^2) per iteration, not O(n^3).
t_sparse <- system.time(
es <- RSpectra::eigs_sym(Lsym, k = 4, which = "SM"))[["elapsed"]]
c(seconds = round(t_sparse, 2),
eigenvalues = paste(signif(sort(es$values), 4), collapse = ", "))#> seconds
#> "0.04"
#> eigenvalues
#> "1.584e-15, 0.00165, 0.001834, 0.003061"# Ascending order, then DROP the trivial lambda_1 = 0 eigenvector
ord <- order(es$values)
U <- es$vectors[, ord[-1], drop = FALSE] # eigenvectors 2, 3, 4
U <- U / pmax(sqrt(rowSums(U^2)), 1e-12)
set.seed(341)
seg2 <- kmeans(U[, 1, drop = FALSE], 2, nstart = 25)$cluster
segmat <- matrix(seg2 - 1, nrow(img), ncol(img))
c(cluster_sizes = paste(table(seg2), collapse = " / "),
fraction_flagged = round(mean(seg2 == which.min(table(seg2))), 4))#> cluster_sizes fraction_flagged
#> "2519 / 1577" "0.385"plot_ly() |>
add_surface(z = img, showscale = FALSE, opacity = 0.55,
colorscale = "Greys", reversescale = TRUE,
name = "Original") |>
add_surface(z = segmat * 0.6 + 1.3, showscale = FALSE, opacity = 0.9,
colorscale = "Portland", name = "Spectral segmentation") |>
layout(title = "Spectral segmentation (top) over the original intensity surface (bottom)",
scene = list(xaxis = list(title = "column"),
yaxis = list(title = "row"),
zaxis = list(title = "Intensity / label")))seg_df <- expand.grid(row = seq_len(nrow(img)), col = seq_len(ncol(img)))
seg_df$intensity <- as.vector(img)
seg_df$segment <- factor(as.vector(segmat))
# Boundary pixels: those with a differently-labelled 4-neighbour
is_boundary <- function(M) {
b <- matrix(FALSE, nrow(M), ncol(M))
b[-1, ] <- b[-1, ] | (M[-1, ] != M[-nrow(M), ])
b[, -1] <- b[, -1] | (M[, -1] != M[, -ncol(M)])
b
}
seg_df$boundary <- as.vector(is_boundary(segmat))
ggplot(seg_df, aes(col, nrow(img) - row)) +
geom_raster(aes(fill = intensity)) +
geom_point(data = subset(seg_df, boundary), color = "#D8433B", size = 0.35) +
scale_fill_gradient(low = "black", high = "white", guide = "none") +
coord_fixed() +
labs(title = "Segmentation boundary overlaid on the MRI slice",
subtitle = "Red: the spectral partition boundary, from the Fiedler vector",
x = NULL, y = NULL) +
theme_void(base_size = 10) +
theme(plot.title = element_text(face = "bold"))
# --- Interactive equivalent ------------------------------------------------
plot_ly(z = t(img), type = "heatmap", colorscale = "Greys",
name = "Original", opacity = 0.8) |>
add_trace(z = t(segmat), type = "heatmap", name = "Segmentation",
opacity = 0.4) |>
layout(title = "Original image with spectral segmentation overlay") |>
hide_colorbar()21.2 Point-cloud segmentation
The SOCR knee-pain data records the reported location of knee pain across four views (left/right, front/back), and it carries ground-truth labels, so the clustering can be validated rather than merely displayed.
knee_raw <- read_html("https://wiki.socr.umich.edu/index.php/SOCR_Data_KneePainData_041409") |>
html_nodes("table") |> _[[2]] |> html_table()
c(rows = nrow(knee_raw), columns = paste(names(knee_raw), collapse = ", "))#> rows columns
#> "8666" "x, Y, View"rescale01 <- function(x) (x - min(x)) / (max(x) - min(x))
# Keep View as a FACTOR with its labels intact -- cbind() would silently
# replace them with integer codes and discard the level names.
knee <- data.frame(
x = rescale01(knee_raw$x),
y = rescale01(knee_raw$Y),
view = factor(knee_raw$View))
levels(knee$view)#> [1] "LB" "LF" "RB" "RF"#>
#> LB LF RB RF
#> 924 3369 882 3491library(kernlab)
set.seed(351)
knee_s <- knee[sample(nrow(knee), 1000), ]
sp_knee <- specc(as.matrix(knee_s[, c("x", "y")]), centers = 4)
km_knee <- kmeans(knee_s[, c("x", "y")], 4, nstart = 25)
data.frame(
method = c("Spectral", "k-means"),
ARI_vs_view = round(c(
mclust::adjustedRandIndex(sp_knee@.Data, knee_s$view),
mclust::adjustedRandIndex(km_knee$cluster, knee_s$view)), 4))bind_rows(
data.frame(knee_s, cl = factor(sp_knee@.Data), method = "Spectral"),
data.frame(knee_s, cl = factor(km_knee$cluster), method = "k-means"),
data.frame(knee_s, cl = knee_s$view, method = "True view (held back)")) |>
ggplot(aes(x, y, color = cl)) +
geom_point(size = 0.7, alpha = 0.7) +
facet_wrap(~ method) + coord_fixed() +
scale_color_brewer(palette = "Set1", guide = "none") +
labs(title = "Knee-pain point cloud: two methods against the true view labels",
subtitle = "The ground truth is available here, so the comparison is quantitative rather than visual",
x = NULL, y = NULL) +
theme_dspa(9)
# --- Interactive equivalent ------------------------------------------------
plot_ly(x = ~knee_s$x, y = ~knee_s$y, type = "scatter", mode = "markers",
symbol = ~factor(sp_knee@.Data),
symbols = c("circle", "x", "square", "diamond"),
color = ~factor(sp_knee@.Data), marker = list(size = 8)) |>
layout(title = "Knee-pain data with spectral cluster labels",
xaxis = list(title = "X"), yaxis = list(title = "Y")) |>
hide_colorbar()The ARI is the finding. A side-by-side plot invites the eye to see agreement; the index measures it, corrected for chance and for the arbitrary labelling of clusters. Here it shows that both methods recover the anatomical views only partially, the pain locations overlap between views, so the geometric clusters and the anatomical labels are genuinely different partitions.
22 PART V: SYNTHESIS
23 Computational complexity summary
\(n\) = observations, \(d\) = features, \(k\) = clusters, \(i\) = iterations, \(B\) = bootstrap or reference replicates, \(m\) = non-zero entries in a sparse graph.
| Method | Time | Memory | Scales to | Note |
|---|---|---|---|---|
| \(k\)-means (Lloyd) | \(O(nkdi)\) | \(O(nd+kd)\) | \(10^7\) | With nstart, multiply by nstart |
| \(k\)-means++ seeding | \(O(nkd)\) | \(O(kd)\) | \(10^7\) | \(O(\log k)\) expected approximation |
| Mini-batch \(k\)-means | \(O(bkdi)\), \(b\ll n\) | \(O(bd)\) | \(10^9\) | Streaming; slight accuracy cost |
| PAM (classical) | \(O(k(n-k)^2 i)\) | \(\mathbf{O(n^2)}\) | \(10^4\) | Robust; works from any dissimilarity |
| CLARA | \(O(k s^2 + k(n-k))\) | \(O(s^2)\) | \(10^6\) | PAM on subsamples of size \(s\) |
| Gower dissimilarity | \(O(n^2 d)\) | \(\mathbf{O(n^2)}\) | \(10^4\) | Required for mixed types |
| Hierarchical (general) | \(O(n^3)\) | \(\mathbf{O(n^2)}\) | \(10^4\) | Lance–Williams updates in place |
| Hierarchical (SLINK/CLINK) | \(O(n^2)\) | \(O(n)\) | \(10^5\) | Single and complete linkage only |
| DBSCAN (with index) | \(O(n\log n)\) | \(O(n)\) | \(10^6\) | Degrades to \(O(n^2)\) above \(d\approx20\) |
| HDBSCAN | \(O(n\log n)\) | \(O(n)\) | \(10^6\) | Builds the full density hierarchy |
| GMM-EM (full \(\Sigma\)) | \(O(nkd^2 i)\) | \(O(nk+kd^2)\) | \(10^6\) | \(O(nkdi)\) for diagonal \(\Sigma\) |
Mclust (all 14 models) |
\(14 K_{\max}\times\) EM | \(O(nk+kd^2)\) | \(10^5\) | BIC over model and \(k\) |
| Spectral (dense) | \(\mathbf{O(n^3)}\) | \(\mathbf{O(n^2)}\) | \(5\times10^3\) | Full eigen() — avoid |
| Spectral (sparse, \(k\) eigenpairs) | \(O(mk^2 i)\) | \(O(m+nk)\) | \(10^5\) | RSpectra::eigs_sym(which = "SM") |
| Silhouette | \(O(n^2)\) | \(O(n^2)\) | \(10^4\) | Needs the full distance matrix |
| Gap statistic | \(B\times\) (clustering cost) | as above | — | \(B=50\)–\(100\) typical |
clusterboot |
\(B\times\) (clustering cost) | as above | — | Plus a matching step per replicate |
Four consequences govern practice.
The \(O(n^2)\) memory wall recurs everywhere. Hierarchical clustering, PAM, Gower, silhouette, and dense spectral clustering all need an \(n\times n\) matrix. At \(n=10^5\) that is 80 GB, the same wall kernel methods hit (Chapter 6, §6.19) and all-pairs similarity hits (Chapter 7, §7.7).
\(k\)-means is the only
method here that is linear in \(n\). That, not superior accuracy,
is why it dominates large-scale practice, and it is why
nstart = 25 is affordable when a single hierarchical fit
would not be.
Sparse construction plus an iterative eigensolver converts \(O(n^3)\) into \(O(mk^2)\). For a local similarity graph, \(m = O(n)\), so the saving is three orders of magnitude at \(n = 4{,}096\) and grows from there.
Validation costs more than fitting. The gap statistic and bootstrap stability each multiply the clustering cost by \(B\). Budget for that: the validation is what distinguishes a finding from an artifact, so it is not the step to economize on.
24 Common pitfalls
| # | Pitfall | Consequence | Fix |
|---|---|---|---|
| 1 | Treating a partition as a finding | Algorithms always return \(k\) groups, including on noise | Validate before interpreting |
| 2 | Reading an average silhouette of 0.2 as success | Below 0.26 means no substantial structure | Use the Kaufman–Rousseeuw thresholds |
| 3 | Narrating centroids without validation | \(k\)-means separates centroids by construction | Interpret only when silhouette, ARI, or stability supports it |
| 4 | Balanced cluster sizes read as evidence | \(k\)-means balances sizes by construction | Sizes say nothing about the data |
| 5 | Scaling a nominal code and using Euclidean distance | Asserts a false ordering and spacing between levels | Gower distance; keep factors as factors |
| 6 | Standardizing a sentinel code (8 = “don’t know”) | Missing cases become the most influential points | Convert to NA; Gower down-weights them |
| 7 | Hard-coding a single-cell imputation | Silently edits the wrong cell if the data changes | Code the rule, apply it to all affected rows |
| 8 | Leaving nstart = 1 |
One draw from a multimodal objective | nstart = 25 or more |
| 9 | Switching algorithm mid-analysis |
Objectives from different heuristics are not comparable | Pick one; compare deliberately if at all |
| 10 | \(k=\sqrt{n/2}\) as a rule | No theoretical basis | Elbow, silhouette, gap, stability |
| 11 | Choosing \(k\) from an elbow alone | Real WCSS curves are often smooth; readers disagree | Gap statistic — it can also return \(k=1\) |
| 12 | Assuming \(k\ge2\) | No partitioning method can report “no clusters” | The gap statistic can |
| 13 | High silhouette read as “correct” | Silhouette rewards the geometry \(k\)-means optimizes | Pair with external labels or a different method |
| 14 | \(k\)-means on elongated or nested shapes | Confidently wrong partition | DBSCAN, spectral, or GMM with full \(\Sigma\) |
| 15 | \(k\)-means with outliers present | A centroid is spent on the contamination | PAM, or trim first |
| 16 | Comparing linkages at different \(k\) | Measures the effect of \(k\), not the linkage | Cut every dendrogram at a common \(k\) |
| 17 | ward.D on unsquared Euclidean distances |
Implements a different criterion | ward.D2 |
| 18 | Eyeballing a cluster-versus-truth table | Labels are arbitrary permutations; chance agreement uncorrected | Adjusted Rand index |
| 19 | Validating a model against its own input | Circular; agreement guaranteed | Hold the validation variable out of the fit |
| 20 | A single global \(\varepsilon\) for DBSCAN | Fails when cluster densities differ | HDBSCAN |
| 21 | Including the \(\lambda_1=0\) eigenvector | Constant column; shifts the embedding for \(k>2\) | Use eigenvectors \(2,\dots,k+1\) |
| 22 | eigen(L_rw, symmetric = TRUE) |
\(L_{\mathrm{rw}}\) is not symmetric; decomposes a different matrix | Use \(L_{\mathrm{sym}}\), or drop
symmetric |
| 23 | Dense \(W\) and full
eigen() for a local graph |
\(O(n^3)\) where \(O(mk^2)\) suffices | sparseMatrix + eigs_sym(which = "SM") |
| 24 | Clustering in high dimension without reduction | Distance concentration destroys the signal | PCA or UMAP first; select features |
25 Practice problems
25.1 Problem 1: Show that
nstart matters
Quantify how often a single-start \(k\)-means run reaches the best solution, as a function of \(k\).
Solution
set.seed(401)
nstart_study <- function(k, reps = 100) {
v <- vapply(seq_len(reps), function(r) {
set.seed(1000 + r); kmeans(Xg, k, nstart = 1)$tot.withinss
}, numeric(1))
best <- kmeans(Xg, k, nstart = 100)$tot.withinss
c(k = k, pct_optimal = 100 * mean(v < best * 1.0001),
worst_excess_pct = 100 * (max(v) / best - 1))
}
ns <- as.data.frame(do.call(rbind, lapply(2:8, nstart_study)))
round(ns, 2)ggplot(ns, aes(k, pct_optimal)) +
geom_line(linewidth = 1, color = "steelblue") + geom_point(size = 2.4) +
scale_x_continuous(breaks = 2:8) +
labs(title = "How often does one random start find the best solution?",
subtitle = "The failure rate grows with k, because the objective has more local minima",
x = "k", y = "Percent of single starts reaching the optimum") +
theme_dspa()
nstart matters more exactly when you are least able to
check the answer by eye.
25.2 Problem 2: Silhouette on noise
Show that a respectable silhouette is achievable on data with no structure, and that the gap statistic is not fooled.
Solution
set.seed(411)
noise_study <- function(d) {
X <- matrix(runif(400 * d), 400, d)
sils <- vapply(2:6, \(k)
mean(silhouette(kmeans(X, k, nstart = 25)$cluster, dist(X))[, "sil_width"]),
numeric(1))
gp <- clusGap(X, FUN = kmeans, nstart = 20, K.max = 6, B = 30)
c(d = d, best_silhouette = max(sils),
gap_selects_k = maxSE(gp$Tab[, "gap"], gp$Tab[, "SE.sim"], "firstSEmax"))
}
ns2 <- as.data.frame(do.call(rbind, lapply(c(2, 3, 5, 10), noise_study)))
ns2$verdict <- as.character(sil_verdict(ns2$best_silhouette))
round(ns2[, 1:3], 4); ns2$verdict#> [1] "weak / possibly artificial" "weak / possibly artificial"
#> [3] "no substantial structure" "no substantial structure"25.3 Problem 3: Gower against naive scaling
Construct mixed-type data with known groups and compare the two distances by ARI.
Solution
set.seed(421)
n3 <- 300
grp <- rep(1:3, each = n3 / 3)
mix3 <- data.frame(
age = rnorm(n3, mean = c(25, 45, 65)[grp], sd = 5),
score = rnorm(n3, mean = c(2, 5, 3)[grp], sd = 1),
region = factor(c("north", "south", "east")[grp]), # informative nominal
noise = factor(sample(c("a", "b", "c"), n3, TRUE))) # uninformative nominal
d_g3 <- daisy(mix3, metric = "gower")
naive3 <- data.frame(lapply(mix3, as.numeric))
d_e3 <- dist(scale(naive3))
data.frame(
distance = c("Gower (types respected)", "Scaled Euclidean (codes as numbers)"),
ARI = round(c(mclust::adjustedRandIndex(pam(d_g3, 3, diss = TRUE)$clustering, grp),
mclust::adjustedRandIndex(pam(d_e3, 3, diss = TRUE)$clustering, grp)), 4),
silhouette = round(c(pam(d_g3, 3, diss = TRUE)$silinfo$avg.width,
pam(d_e3, 3, diss = TRUE)$silinfo$avg.width), 4))region, and, worse, the
same treatment on noise, where it manufactures spurious
distance from a variable that carries no signal at all.
25.4 Problem 4: Verify the \(k\)-means++ guarantee empirically
Compare the distribution of objectives from \(k\)-means++ seeding against the \(8(\ln k+2)\) bound.
Solution
set.seed(431)
kpp_study <- function(k, reps = 100) {
opt <- kmeans(Xg, k, nstart = 200)$tot.withinss
seed_only <- vapply(seq_len(reps), function(r) {
C <- kmeanspp_init(Xg, k, seed = 5000 + r)
D <- as.matrix(dist(rbind(C, Xg)))[-(1:k), 1:k, drop = FALSE]
sum(apply(D, 1, min)^2) # objective BEFORE Lloyd runs
}, numeric(1))
c(k = k, bound = 8 * (log(k) + 2),
mean_ratio = mean(seed_only) / opt,
max_ratio = max(seed_only) / opt)
}
kp <- as.data.frame(do.call(rbind, lapply(2:8, kpp_study)))
round(kp, 3)kp |> pivot_longer(c(bound, mean_ratio, max_ratio),
names_to = "quantity", values_to = "value") |>
ggplot(aes(k, value, color = quantity)) +
geom_line(linewidth = 1) + geom_point(size = 2.2) +
scale_y_log10() + scale_x_continuous(breaks = 2:8) +
scale_color_brewer(palette = "Set1") +
labs(title = "k-means++ seeding against its theoretical bound",
subtitle = "Ratios are objective / optimum, measured BEFORE Lloyd's algorithm runs",
x = "k", y = "Ratio (log scale)", color = NULL) +
theme_dspa()
25.5 Problem 5: Predict the failure
Given four structures, predict which method will win before running anything, then check.
Solution
set.seed(441)
predictions <- data.frame(
structure = names(cases),
predicted_winner = c("GMM (full covariance handles elongation)",
"Spectral / DBSCAN (connectivity, not compactness)",
"HDBSCAN (varying density)",
"None -- gap statistic should return k = 1"))
results <- do.call(rbind, lapply(seq_along(cases), function(i) {
X <- cases[[i]][[1]]; y <- cases[[i]][[2]]; k <- cases[[i]][[3]]
km <- kmeans(X, k, nstart = 25)
gm <- Mclust(X, G = k, verbose = FALSE)
hb <- dbscan::hdbscan(X, minPts = 8)
Xs <- X[sample(nrow(X), min(400, nrow(X))), ]
ys <- y[sample(nrow(X), min(400, nrow(X)))]
sp <- spectral_cluster(Xs, k, sigma = 0.6)
data.frame(structure = names(cases)[i],
kmeans = round(mclust::adjustedRandIndex(km$cluster, y), 3),
GMM = round(mclust::adjustedRandIndex(gm$classification, y), 3),
HDBSCAN = round(mclust::adjustedRandIndex(hb$cluster, y), 3),
spectral = round(mclust::adjustedRandIndex(sp$cluster, ys), 3))
}))
predictions25.6 Problem 6: Cophenetic correlation across linkages
Determine which linkage most faithfully preserves the original distances, and whether that makes it the best clustering.
Solution
set.seed(451)
Xc6 <- Xg
Dc6 <- dist(Xc6)
links6 <- c("single", "complete", "average", "centroid", "ward.D2")
do.call(rbind, lapply(links6, function(m) {
hc <- hclust(Dc6, method = m)
cl <- cutree(hc, k = 3)
data.frame(linkage = m,
cophenetic = round(cor(Dc6, cophenetic(hc)), 4),
silhouette = round(mean(silhouette(cl, Dc6)[, "sil_width"]), 4),
ARI_vs_truth = round(mclust::adjustedRandIndex(cl, g), 4))
}))25.7 Problem 7: The trivial eigenvector
Show that including the \(\lambda_1=0\) eigenvector is harmless at \(k=2\) and harmful at \(k>2\).
Solution
spectral_variant <- function(X, k, sigma, drop_trivial = TRUE, seed = 1) {
W <- build_W(X, sigma); d <- rowSums(W); d[d == 0] <- 1e-12
Dm <- diag(1 / sqrt(d))
Lsym <- diag(nrow(W)) - Dm %*% W %*% Dm
e <- eigen(Lsym, symmetric = TRUE)
nc <- ncol(e$vectors)
idx <- if (drop_trivial) (nc - k):(nc - 1) else (nc - k + 1):nc
U <- e$vectors[, idx, drop = FALSE]
U <- U / pmax(sqrt(rowSums(U^2)), 1e-12)
set.seed(seed); kmeans(U, k, nstart = 25)$cluster
}
set.seed(461)
# Three well-separated blobs, where k = 3 is the right answer
X7 <- Xg[sample(nrow(Xg), 300), ]; y7 <- g[sample(nrow(Xg), 300)]
do.call(rbind, lapply(c(2, 3, 4), function(k) {
a <- spectral_variant(X7, k, 1.2, drop_trivial = TRUE)
b <- spectral_variant(X7, k, 1.2, drop_trivial = FALSE)
data.frame(k = k,
ARI_dropping_trivial = round(mclust::adjustedRandIndex(a, y7), 4),
ARI_including_trivial = round(mclust::adjustedRandIndex(b, y7), 4),
partitions_identical = mclust::adjustedRandIndex(a, b) > 0.999)
}))25.8 Problem 8: Sparse versus dense spectral clustering
Measure the time and memory saved by sparse construction on a local similarity graph.
Solution
library(Matrix); library(RSpectra)
set.seed(471)
bench_spectral <- function(side) {
n <- side^2
co <- expand.grid(r = 1:side, c = 1:side)
co$v <- as.vector(matrix(runif(n), side, side))
off <- expand.grid(dr = -2:2, dc = -2:2); off <- off[!(off$dr == 0 & off$dc == 0), ]
ed <- do.call(rbind, lapply(seq_len(nrow(off)), function(o) {
r2 <- co$r + off$dr[o]; c2 <- co$c + off$dc[o]
ok <- r2 >= 1 & r2 <= side & c2 >= 1 & c2 <= side
j <- (c2[ok] - 1) * side + r2[ok]
data.frame(i = which(ok), j = j, w = exp(-(co$v[ok] - co$v[j])^2 / 0.01))
}))
Wsp <- sparseMatrix(i = ed$i, j = ed$j, x = ed$w, dims = c(n, n))
Wsp <- (Wsp + Matrix::t(Wsp)) / 2
dg <- Matrix::rowSums(Wsp); dg[dg == 0] <- 1e-12
Ds <- Diagonal(x = 1 / sqrt(dg))
Ls <- Diagonal(n) - Ds %*% Wsp %*% Ds
t_sp <- system.time(RSpectra::eigs_sym(Ls, k = 3, which = "SM"))[["elapsed"]]
t_de <- if (n <= 1600) system.time(eigen(as.matrix(Ls), symmetric = TRUE))[["elapsed"]] else NA
c(n = n, density = signif(length(Wsp@x) / n^2, 3),
sparse_MB = round(as.numeric(object.size(Wsp)) / 1e6, 2),
dense_MB = round(8 * n^2 / 1e6, 1),
sparse_sec = round(t_sp, 3), dense_sec = round(t_de, 3))
}
bs <- as.data.frame(do.call(rbind, lapply(c(20, 30, 40), bench_spectral)))
bsggplot(bs, aes(n)) +
geom_line(aes(y = dense_MB, color = "Dense"), linewidth = 1) +
geom_line(aes(y = sparse_MB, color = "Sparse"), linewidth = 1) +
geom_point(aes(y = dense_MB, color = "Dense"), size = 2.2) +
geom_point(aes(y = sparse_MB, color = "Sparse"), size = 2.2) +
scale_y_log10() +
scale_color_manual(values = c(Dense = "#D8433B", Sparse = "#3B7DD8")) +
labs(title = "Memory for the similarity matrix, sparse vs. dense",
subtitle = "Local graph, 5x5 spatial window. The gap grows as n^2 / n",
x = "Pixels n", y = "Megabytes (log scale)", color = NULL) +
theme_dspa()
26 Checkpoint
- A colleague reports three clusters with an average silhouette of 0.21 and describes each group’s personality. What do you say?
- Your data has a nominal variable coded 1, 2, 3 and a “don’t know” code of 8. What happens if you standardize and use Euclidean distance?
kmeans(X, 5)gives a different answer every time you run it. Is the function broken?- Which method for choosing \(k\) can return \(k=1\), and why does that matter?
- Your clusters have an average silhouette of 0.55 but an ARI against known labels of 0.02. What has happened?
eigen(L, symmetric = TRUE)on the random-walk Laplacian runs without error. Is the result correct?
Answers
- The silhouette says there is no substantial structure. Under the standard Kaufman–Rousseeuw thresholds, anything below 0.26 means “no substantial structure has been found”, 0.21 is in that band. The distinct-looking centroids are not evidence: \(k\)-means separates centroids by construction and produces equally convincing ones on uniform noise. Ask for the gap statistic (which can return \(k=1\)) and a bootstrap stability check, and hold the narrative until one of them supports it.
- Two separate failures. The nominal code acquires a false
ordering and spacing, level 2 becomes exactly halfway between 1
and 3, and 1 becomes twice as far from 3 as from 2, none of which is
meaningful when the labels could have been assigned in any order. And
the sentinel 8 becomes a real number roughly two standard
deviations from everything else, so the respondents whose
answer is unknown become the most influential points in the
analysis. Convert 8 to
NAand use Gower distance, which handles nominal variables by simple matching and missingness by down-weighting. - No, the default is
nstart = 1. The \(k\)-means objective is non-convex and the problem is NP-hard, so Lloyd’s algorithm converges to whichever local minimum its random start falls into. A single run is one draw from a multimodal distribution. Setnstart = 25(or more) and the function runs the algorithm repeatedly and returns the best solution found; results then stabilize. Also set a seed. - The gap statistic. The elbow and silhouette compare \(k\) against other values of \(k\), so both must return some \(k\ge2\); neither can express “there are no clusters.” The gap statistic compares the observed within-cluster dispersion against its expectation under an explicit no-structure reference null, so if the observed curve does not separate from the null the 1-SE rule selects \(k=1\). That is the only method here that can report the absence of structure, which is often the correct answer.
- The clusters are geometrically tidy and scientifically irrelevant. Silhouette measures compactness and separation under the metric being used, and \(k\)-means optimizes exactly that, so a high silhouette partly measures how well the algorithm did its own job. ARI near zero means the partition does not correspond to the known grouping. The nested-rings case is the canonical example: cutting concentric circles into pie slices gives compact, well-separated clusters that are completely wrong. Always pair an internal index with external labels, a stability check, or a method making different assumptions.
- No. \(L_{\mathrm{rw}}=I-D^{-1}W\) is not
symmetric unless every vertex has the same degree.
symmetric = TRUEinstructs LAPACK to read only the lower triangle and assume the rest by reflection, so it silently eigendecomposes a different matrix than the one supplied, no error, wrong answer. Use \(L_{\mathrm{sym}}=I-D^{-1/2}WD^{-1/2}\), which is symmetric, and recover the random-walk eigenvectors as \(u_{\mathrm{rw}}=D^{-1/2}u_{\mathrm{sym}}\) if you need them.
27 Summary
The clustering problem
- A clustering algorithm always returns \(k\) groups, including on pure noise. Its output is a hypothesis; validation decides whether it is a finding.
- Distance is the most consequential choice, made before any algorithm runs. Gower handles mixed types; scaling a nominal or sentinel code and using Euclidean distance imposes meaning that is not there.
- Distances concentrate in high dimensions, so reduce dimension first.
- Silhouette has standard thresholds: below 0.26 means no substantial structure. Report ARI for external agreement and bootstrap Jaccard for stability.
Partitioning
- \(k\)-means minimizes within-cluster sum of squares. The problem is NP-hard; Lloyd’s algorithm converges monotonically to a local optimum.
nstartis the single most important argument. The default of 1 gives one draw from a multimodal objective.- \(k\)-means++ carries an \(O(\log k)\) expected-approximation guarantee; uniform random seeding has none.
- The gap statistic compares against a no-structure null and is the only method here that can return \(k=1\).
- \(k\)-means assumes spherical, equal-variance, similarly-sized clusters and fails predictably when any is violated. PAM is the robust alternative and pairs naturally with Gower.
Hierarchical and density-based
- All linkages are instances of the Lance–Williams recurrence, which is why the algorithm needs only the dissimilarity matrix, and why it costs \(O(n^2)\) memory.
- Use
ward.D2with ordinary Euclidean distances; compare linkages at a common \(k\). - DBSCAN finds arbitrarily shaped clusters, discovers their number, and labels noise explicitly. Its one global \(\varepsilon\) fails when densities differ; HDBSCAN removes that assumption.
Model-based and spectral
- A Gaussian mixture gives soft assignments, principled model selection by BIC over 14 covariance parameterizations, and shapes controlled by explicit parameters. \(k\)-means is the \(\Sigma_\ell=\sigma^2I\), \(\sigma^2\to0\) limit.
- EM alternates responsibilities (E) with responsibility-weighted moments (M) and increases the likelihood monotonically.
- There are three Laplacians. \(L\) and \(L_{\mathrm{sym}}\) are symmetric; \(L_{\mathrm{rw}}\) is not, and passing it to a symmetric solver silently decomposes a different matrix.
- \(\lambda_1=0\) with the constant eigenvector carries no partition information and must be dropped; the multiplicity of 0 equals the number of connected components.
- Local similarity graphs are sparse by construction. Sparse storage plus an iterative eigensolver turns \(O(n^3)\) into \(O(mk^2)\).
Where these threads continue
| Thread | Continues in |
|---|---|
| Resampling, calibration, and honest model comparison | Model assessment |
| Sparse penalties for feature selection in clustering | Feature selection |
| Trajectory clustering and change-point detection | Longitudinal analysis |
| EM as a general optimization strategy | Function optimization |
| Autoencoders and learned representations for clustering | Deep learning |
Further practice: apply these methods to other case studies in the DSPA archive, and in particular repeat the Boys Town youth-development analysis (Chapter 5, §5.10) as a clustering problem — clustering on GPA, alcohol use and attitudes, parental closeness, and delinquency, then checking whether the recovered groups relate to variables held out of the fit.
28 Chapter roadmap
- Chapter 1:
Foundations. R toolchain, reproducibility conventions,
dspa_read(), simulation. - Chapter 2: Data quality and exploratory visual analytics. Missingness, robust statistics, graphical perception.
- Chapter 3: Linear algebra, matrix computing, and regression. Eigendecomposition, projection, sparse matrices, complexity.
- Chapter 4: Dimensionality reduction. PCA, distance concentration, MDS, t-SNE, UMAP.
- Chapter 5: Supervised classification. Bayes error, evaluation metrics, leakage, trees and forests.
- Chapter 6: Black-box methods. Neural networks, kernels, SVM, ensembles, boosting.
- Chapter 7: Text mining and association rules. TF-IDF, similarity, Apriori, multiplicity control.
- Model assessment, validation, improvement. Nested resampling, calibration, cost-sensitive learning.
- Variable importance and feature selection. Ridge, LASSO, elastic net, stability selection.
- Longitudinal and time-series analysis. Mixed models, ARIMA, recurrent architectures.
- Function optimization. Gradient descent, EM, duality, Bayesian optimization.
- Deep learning. Autoencoders, convolutional and recurrent networks, representation learning.
29 Session information
#> R version 4.3.3 (2024-02-29 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 11 x64 (build 26200)
#>
#> Matrix products: default
#>
#>
#> locale:
#> [1] LC_COLLATE=English_United States.utf8
#> [2] LC_CTYPE=English_United States.utf8
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United States.utf8
#>
#> time zone: America/New_York
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] rvest_1.0.4 mclust_6.1.1 fpc_2.2-12 factoextra_1.0.7
#> [5] cluster_2.1.6 plotly_4.12.1 patchwork_1.3.0 tidyr_1.3.1
#> [9] dplyr_1.1.4 ggplot2_4.0.1
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 xfun_0.52 bslib_0.9.0 htmlwidgets_1.6.4
#> [5] websocket_1.4.1 ggrepel_0.9.5 processx_3.8.6 lattice_0.22-6
#> [9] crosstalk_1.2.1 vctrs_0.6.5 tools_4.3.3 ps_1.9.0
#> [13] generics_0.1.3 stats4_4.3.3 curl_6.2.0 flexmix_2.3-19
#> [17] parallel_4.3.3 tibble_3.2.1 DEoptimR_1.1-3 pkgconfig_2.0.3
#> [21] Matrix_1.6-5 data.table_1.16.4 RColorBrewer_1.1-3 S7_0.2.1
#> [25] lifecycle_1.0.5 stringr_1.5.1 compiler_4.3.3 farver_2.1.2
#> [29] chromote_0.4.0 htmltools_0.5.8.1 class_7.3-22 sass_0.4.9
#> [33] yaml_2.3.10 pillar_1.10.1 later_1.4.1 jquerylib_0.1.4
#> [37] prabclus_2.3-3 MASS_7.3-60.0.1 diptest_0.77-1 cachem_1.1.0
#> [41] nlme_3.1-165 robustbase_0.99-2 tidyselect_1.2.1 digest_0.6.37
#> [45] stringi_1.8.4 purrr_1.0.2 kernlab_0.9-32 splines_4.3.3
#> [49] labeling_0.4.3 fastmap_1.2.0 grid_4.3.3 cli_3.6.3
#> [53] magrittr_2.0.3 withr_3.0.2 scales_1.4.0 promises_1.3.2
#> [57] rmarkdown_2.31 httr_1.4.7 otel_0.2.0 nnet_7.3-19
#> [61] modeltools_0.2-23 evaluate_1.0.3 knitr_1.51 viridisLite_0.4.2
#> [65] mgcv_1.9-1 rlang_1.1.5 Rcpp_1.0.14 glue_1.8.0
#> [69] selectr_0.4-2 xml2_1.3.6 rstudioapi_0.18.0 jsonlite_1.8.9
#> [73] R6_2.6.1