Clustering, hands on

Clustering groups records that resemble each other when nobody has labelled anything. It is genuinely useful and easy to misuse, because the algorithm will hand you groups whether or not any exist. This section covers running k-means, why scaling the features first is not optional, choosing the number of groups with evidence rather than by eye, a density-based method for data that does not come in tidy blobs, and judging the result by reading what landed in each group.

k-means#

from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score

X_scaled = StandardScaler().fit_transform(X)       # otherwise the largest-valued feature dominates
for k in range(2, 11):
    labels = KMeans(n_clusters=k, n_init="auto", random_state=0).fit_predict(X_scaled)
    print(k, round(silhouette_score(X_scaled, labels), 3))

k-means places k centres and assigns each point to the nearest, then moves each centre to the middle of its points, repeating until stable. It assumes roughly round clusters of similar size. The silhouette score compares how close points are to their own cluster against the next nearest; higher is better, and it helps choose k.

The scaling line in that snippet is doing more work than it looks. Cluster 300 people on income in euros and age in years without scaling and the distance between two people is, for practical purposes, the difference in their incomes: a gap of 20,000 euros swamps a gap of 30 years. The result is three income bands whose average ages are 39.8, 39.9 and 39.9 — the age column had no effect whatsoever on a result presented as a segmentation by income and age.

import numpy as np
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler

rng = np.random.default_rng(0)
income = rng.normal(50_000, 15_000, 300)
age = rng.normal(40, 12, 300)
people = np.c_[income, age]

for name, data in (("unscaled", people), ("scaled", StandardScaler().fit_transform(people))):
    labels = KMeans(n_clusters=3, n_init=10, random_state=0).fit_predict(data)
    summary = [(round(income[labels == c].mean()), round(float(age[labels == c].mean()), 1))
               for c in range(3)]
    print(name, summary)

# unscaled [(67168, 39.8), (47117, 39.9), (29254, 39.9)]
# scaled   [(66167, 41.3), (40840, 49.1), (42458, 28.2)]

The n_init argument in that snippet is not decoration either. k-means starts from random centres and only finds a local optimum, so a single unlucky start can leave two centres inside one group and none in another. Run it ten times with different starts and keep the best, which is what n_init does. On six well-separated blobs, ten single-start runs gave inertias of 272, 300, 300, 300, 300, 272, 300, 2649, 3356 and 2350: three of the ten landed somewhere catastrophic, and four more settled ten per cent worse than the best.

from sklearn.datasets import make_blobs

blobs, _ = make_blobs(n_samples=600, centers=6, cluster_std=0.5, random_state=3)

singles = [round(KMeans(n_clusters=6, n_init=1, init="random", random_state=s)
                 .fit(blobs).inertia_, 1) for s in range(10)]
print(singles)          # [299.8, 300.2, 272.3, 2649.1, 299.8, 300.0, 271.9, 3355.7, 299.6, 2350.3]
print(round(KMeans(n_clusters=6, n_init=10, random_state=0).fit(blobs).inertia_, 1))
# 271.9   the best of the ten, which is what n_init gives you

When clusters are not round#

from sklearn.cluster import HDBSCAN

labels = HDBSCAN(min_cluster_size=20).fit_predict(X_scaled)
# -1 means 'noise': points that belong to no dense region

Density-based methods find clusters of any shape and leave outliers unassigned instead of forcing them into a group. You choose a minimum cluster size rather than the number of clusters, which is often easier to reason about.

The two interleaved crescents of the two-moons dataset are the standard case where this matters, and running both methods on it also shows why you cannot pick between them with a score. k-means gets it badly wrong, with an adjusted Rand index against the true grouping of 0.25, and HDBSCAN gets it exactly right, at 1.00. The silhouette scores go the other way: 0.489 for the wrong answer and 0.330 for the right one.

from sklearn.datasets import make_moons
from sklearn.metrics import adjusted_rand_score, silhouette_score

moons, truth = make_moons(n_samples=500, noise=0.06, random_state=0)

km = KMeans(n_clusters=2, n_init=10, random_state=0).fit_predict(moons)
hd = HDBSCAN(min_cluster_size=20).fit_predict(moons)

print(round(adjusted_rand_score(truth, km), 3), round(silhouette_score(moons, km), 3))
# 0.253 0.489      wrong grouping, better silhouette
print(round(adjusted_rand_score(truth, hd), 3), round(silhouette_score(moons, hd), 3))
# 1.0 0.33         correct grouping, worse silhouette

The practical rule that falls out of this: use an internal score to compare runs of the same method at different settings, never to choose between methods that make different shape assumptions. Choosing between methods is done by looking at the result.

what your data looks like                   reach for            the setting you choose
round groups, similar sizes, you know k     k-means              k
round groups, you do not know k             k-means + gap stat   k, by comparison with noise
elliptical or overlapping groups            Gaussian mixture     number of components
odd shapes, varying density, real outliers  HDBSCAN              minimum cluster size
you want a hierarchy to cut at any level    agglomerative        linkage, then where to cut
millions of rows                            MiniBatchKMeans      k, and the batch size

What to skip here: the dozens of other clustering algorithms, the internal indices beyond silhouette and inertia, and any attempt to prove a clustering correct. Two methods with different assumptions, an honest comparison against noise, and the discipline of naming every cluster is the whole of what this subject gives you in practice.

Judge by reading#

Scores say clusters are separated, not that they mean anything. For each cluster, look at typical members and the features that distinguish it from the rest, and try to name it. Clusters you cannot name are often artefacts of scaling or of the method.

Make that a procedure rather than an intention, because it is the step everyone skips. For each cluster print its size, the five features whose mean differs most from the overall mean in standard deviations, and the three members closest to its centre. Then write a name for it in a sentence, in a document. A cluster you cannot describe in a sentence is not a finding, and the number of clusters you can name is a better answer to the question of how many clusters there are than any score on this page.

def describe(data, labels, feature_names, X_original):
    overall = data.mean(axis=0)
    spread = data.std(axis=0)
    for c in sorted(set(labels) - {-1}):
        members = labels == c
        gap = (data[members].mean(axis=0) - overall) / spread     # in standard deviations
        top = np.argsort(-np.abs(gap))[:5]
        print(f"cluster {c}: {members.sum()} members")
        for i in top:
            print(f"   {feature_names[i]:20} {gap[i]:+.2f} sd")
        centre = data[members].mean(axis=0)
        closest = np.argsort(((data[members] - centre) ** 2).sum(axis=1))[:3]
        print("   typical members:", X_original[members][closest])