SOCR ≫ DSPA ≫ DSPA3 Topics ≫ 

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, adjustedRandIndex

How this chapter uses graphics

Every two-dimensional figure is drawn with ggplot2 and rendered statically. Immediately after each one, the equivalent plot_ly() code appears in a chunk marked eval=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:

  1. 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.
  2. Interpret an average silhouette against the standard thresholds, and recognize when a partition provides no evidence of structure.
  3. Measure agreement with known labels using the adjusted Rand index, and explain what the adjustment corrects for.
  4. State the \(k\)-means objective, explain why the problem is NP-hard, and prove that Lloyd’s algorithm converges.
  5. Explain why nstart matters and state the \(O(\log k)\) guarantee that \(k\)-means++ provides.
  6. Select \(k\) using the elbow, silhouette, and gap statistic, and explain what null each compares against.
  7. Predict when \(k\)-means will fail from the shape of the data, and choose a method that will not.
  8. Apply DBSCAN, explain the role of \(\varepsilon\) and minPts, and read a \(k\)-distance plot.
  9. Fit a Gaussian mixture, select its covariance parameterization by BIC, and distinguish soft from hard assignment.
  10. 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

tab_hd <- table(cluster = km_hd$cluster, type = hotdog$Type)
tab_hd
#>        type
#> cluster Beef Meat Poultry
#>       1    0    0       9
#>       2   12    9       3
#>       3    8    8       5
c(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.4699

The 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.000

Compare 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.0014

There 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
round(fit_l$wss_trace, 2)
#>  [1] 943.86 863.06 789.91 750.45 721.78 684.81 665.88 664.29 663.88 663.88

The 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.00
ggplot(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 is nstart = 1. Set nstart = 25 (or more) so the function runs the algorithm repeatedly and returns the best result it found. This costs a factor of nstart in 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")))