7.2 Clustering

Clustering looks for groups of similar cases without being told what the groups are. Every method needs a definition of “similar” and returns groups even when there are none, so the results need checking and judgement.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from matplotlib.axes import Axes
from matplotlib.colors import ListedColormap
from scipy.cluster.hierarchy import dendrogram, fcluster, linkage
from sklearn import datasets
from sklearn.cluster import DBSCAN, AgglomerativeClustering, KMeans
from sklearn.decomposition import PCA
from sklearn.metrics import adjusted_rand_score, silhouette_samples, silhouette_score
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler

from statly import style

style.use()
rng = np.random.default_rng(1)

wine = datasets.load_wine(as_frame=True)  # 178 Italian wines, 13 chemical measurements
cultivar = wine.target.to_numpy()  # three grape cultivars, used only to check results
CULTIVAR_COLORS = np.array([style.BLUE, style.GREEN, style.CORAL])
GROUP_COLORS = np.array(
    [style.GOLD, style.PURPLE, style.GREY, style.BLUE, style.GREEN, style.CORAL]
)
scaled = StandardScaler().fit_transform(wine.data)
plane = PCA(2).fit_transform(
    scaled
)  # the first two principal components, for drawing (Chapter 7.1)

7.2.1 K-means

K-means needs the number of groups, k, and looks for the split of the cases into k groups that makes the groups as tight as possible. Tightness is the within-cluster variation: for each group, the squared distances of its cases to the group’s centre, added up over all groups:

$$W = \sum_{k=1}^{K} \sum_{i \in C_k} \sum_{j=1}^{p} (x_{ij} - \bar x_{kj})^2$$

where $\bar x_{kj}$ is the average of measurement j in group k. Trying every split is impossible (there are about $3^{178}$ ways to split 178 wines into three groups), so the algorithm alternates two steps: give every case to its nearest centre, then move every centre to the middle of its cases. Each step can only lower W, so it stops when nothing changes.

start = plane[rng.choice(len(plane), 3, replace=False)] + np.array(
    [[0.5, 2], [-1, 1], [1, -1.5]]
)


def nearest(centres: np.ndarray) -> np.ndarray:
    """For every wine, the number of its nearest centre."""
    return np.argmin(((plane[:, None, :] - centres[None]) ** 2).sum(axis=2), axis=1)


def kmeans_rounds(centres: np.ndarray, rounds: int) -> tuple[np.ndarray, np.ndarray]:
    """Run k-means by hand for a number of rounds; return the centres and groups."""
    for _ in range(rounds):
        groups = nearest(centres)
        centres = np.array([plane[groups == g].mean(axis=0) for g in range(len(centres))])
    return centres, nearest(centres)


def plot_round(ax: Axes, rounds: int) -> None:
    """The wines coloured by group after some rounds, with the centres."""
    centres, groups = kmeans_rounds(start, rounds)
    ax.scatter(plane[:, 0], plane[:, 1], c=GROUP_COLORS[groups], s=12, alpha=0.7)
    ax.scatter(
        centres[:, 0],
        centres[:, 1],
        marker="X",
        s=160,
        c=GROUP_COLORS[:3],
        edgecolors=style.INK,
        linewidths=1,
    )
    ax.set(
        title="Start" if rounds == 0 else f"After {rounds} round{'s' if rounds > 1 else ''}",
        xlabel="component 1",
    )


fig, axes = style.subplots(3, height=3.1, width=10.5, sharey=True)
for ax, rounds in zip(axes, (0, 1, 10)):
    plot_round(ax, rounds)
axes[0].set_ylabel("component 2")
fig.tight_layout()
plt.show()

Figure 7.2.1: K-means with k = 3 on the wines, drawn on the first two principal components. From a random start (left), the centres (crosses) move to the middle of their wines, the wines join the nearest centre, and after a few rounds nothing changes (right).

7.2.2 Local optima

The algorithm finds a split where neither step can improve things, but not necessarily the best split. Different starts can end in different places:

centres_true = np.array([[0, 0], [4, 0], [0, 4], [4, 4]])
blobs = np.concatenate([c + rng.normal(0, 0.7, (60, 2)) for c in centres_true])


def plot_start(ax: Axes, seed: int) -> None:
    """One k-means run from one random start."""
    model = KMeans(4, n_init=1, init="random", random_state=seed).fit(blobs)
    ax.scatter(blobs[:, 0], blobs[:, 1], c=GROUP_COLORS[model.labels_], s=8, alpha=0.7)
    ax.scatter(
        *model.cluster_centers_.T, marker="X", s=120, c=GROUP_COLORS[:4], edgecolors=style.INK
    )
    ax.set(title=f"W = {model.inertia_:.0f}", xticks=[], yticks=[], aspect="equal")


best_w = KMeans(4, n_init=20, random_state=0).fit(blobs).inertia_
ends = {
    seed: KMeans(4, n_init=1, init="random", random_state=seed).fit(blobs).inertia_
    for seed in range(60)
}
good = [seed for seed, w in ends.items() if w < 1.01 * best_w][:2]
bad = [seed for seed, w in ends.items() if w > 1.3 * best_w][:2]
fig, axes = style.subplots(4, height=2.8, width=10.5)
for ax, seed in zip(axes, good + bad):
    plot_start(ax, seed)
fig.tight_layout()
plt.show()

Figure 7.2.2: K-means with k = 4 on simulated data with four groups, from four different random starts. Two starts find the four groups (lowest within-cluster variation W); the others end with two groups merged and one split in half, a worse but stable answer.

So run k-means from many starts (20 to 50) and keep the split with the lowest W. Better starting rules, such as k-means++ which spreads the first centres out, make bad endings rarer.

7.2.3 How many groups?

The within-cluster variation always falls as k grows (with k = n every case is its own group and W is 0), so it cannot pick k by itself. Two common guides:

• The elbow: plot W against k and look for a bend where more groups stop helping much.
• The silhouette: for each case, let a be its average distance to the other cases in its group and b its average distance to the cases of the nearest other group. Its silhouette is $s = (b - a) / \max(a, b)$: near 1 if it sits well inside its group, near 0 on a border, negative if it would fit better elsewhere. The best k has the highest average silhouette.

ks = range(1, 11)
within = [KMeans(k, n_init=20, random_state=0).fit(scaled).inertia_ for k in ks]
silhouettes = [
    silhouette_score(scaled, KMeans(k, n_init=20, random_state=0).fit_predict(scaled))
    for k in ks
    if k > 1
]
groups3 = KMeans(3, n_init=20, random_state=0).fit_predict(scaled)
each = silhouette_samples(scaled, groups3)


def plot_elbow(ax: Axes) -> None:
    """Within-cluster variation against k."""
    ax.plot(list(ks), within, "o-", color=style.BLUE, ms=5)
    ax.plot(3, within[2], "o", ms=10, color=style.CORAL, fillstyle="none", mew=2)
    ax.set(
        title="Elbow",
        xlabel="number of groups k",
        ylabel="within-cluster variation W",
        xticks=list(ks),
    )


def plot_silhouette_by_k(ax: Axes) -> None:
    """Average silhouette against k."""
    ax.plot(list(ks)[1:], silhouettes, "o-", color=style.GREEN, ms=5)
    ax.set(
        title="Average silhouette",
        xlabel="number of groups k",
        ylabel="average silhouette",
        xticks=list(ks)[1:],
    )


def plot_silhouette_bars(ax: Axes) -> None:
    """Every wine's silhouette for k = 3, sorted within its group."""
    y = 0
    for g in range(3):
        values = np.sort(each[groups3 == g])
        ax.barh(
            np.arange(y, y + len(values)), values, height=1, color=GROUP_COLORS[g], alpha=0.8
        )
        y += len(values) + 6
    ax.axvline(each.mean(), color=style.INK, ls="--", lw=1)
    ax.set(title="Each wine, k = 3", xlabel="silhouette", yticks=[])


fig, axes = style.subplots(3, height=3.1, width=10.5)
plot_elbow(axes[0])
plot_silhouette_by_k(axes[1])
plot_silhouette_bars(axes[2])
fig.tight_layout()
plt.show()

Figure 7.2.3: Left: within-cluster variation W of the wines by k; the bend at k = 3 suggests three groups. Middle: average silhouette by k; also highest at 3. Right: the silhouette of each wine for k = 3, sorted within each group; most are well above 0, a few sit on a border.

Here three groups is a clear answer, and they match the cultivars closely (the adjusted Rand index, a measure of agreement between two groupings, from 0 for chance to 1 for identical, is 0.90). Real data are rarely this kind: elbows are often vague and silhouettes low. Ask whether the groups make sense, not only whether they exist.

7.2.4 Hierarchical clustering

Hierarchical clustering does not need k up front. It starts with every case as its own group and repeatedly joins the two closest groups, until one group holds everything. The record of the joins is a tree, the dendrogram: the height of each join is how different the joined groups were, and cutting the tree at some height gives the groups.

How “closest” is measured between two groups is the linkage:

• complete: the distance between their farthest members;
• single: between their nearest members, which tends to make long chains;
• average: the average over all pairs;
• Ward’s: the increase in within-cluster variation the join would cause, which gives compact groups like k-means.

sample = np.concatenate(
    [rng.choice(np.where(cultivar == c)[0], 20, replace=False) for c in range(3)]
)
LINKAGES = ["ward", "complete", "average", "single"]


def plot_dendrogram(ax: Axes, method: str) -> None:
    """A dendrogram with leaves coloured by cultivar."""
    tree = linkage(scaled[sample], method=method)
    parts = dendrogram(
        tree, ax=ax, color_threshold=0, above_threshold_color=style.INK, no_labels=True
    )
    leaves = cultivar[sample][parts["leaves"]]
    bottom = -0.04 * tree[:, 2].max()
    ax.scatter(
        np.arange(len(leaves)) * 10 + 5,
        np.full(len(leaves), bottom),
        c=CULTIVAR_COLORS[leaves],
        s=10,
        clip_on=False,
    )
    agreement = adjusted_rand_score(cultivar[sample], fcluster(tree, 3, criterion="maxclust"))
    ax.set(title=f"{method}: agreement {agreement:.2f}", ylim=(bottom * 2, None))
    ax.grid(False)


fig, axes = style.subplots(4, cols=2, height=2.7, width=10)
for ax, method in zip(axes, LINKAGES):
    plot_dendrogram(ax, method)
fig.tight_layout()
plt.show()

Figure 7.2.4: Dendrograms of the same 60 wines (20 per cultivar) with four linkages; the dots under each tree show the cultivars. Ward’s tree splits cleanly into branches that follow the cultivars; complete and average linkage capture them only partly; single linkage chains the wines on one by one and finds no groups at all.

The “agreement” in each title is the adjusted Rand index between three groups cut from the tree and the three cultivars.

local = np.random.default_rng(4)
group_a = local.normal([0, 0], 0.5, (5, 2))
group_b = local.normal([3, 1], 0.5, (5, 2))
pair_distance = np.linalg.norm(group_a[:, None] - group_b[None], axis=2)


def setup_pair(ax: Axes, title: str) -> None:
    """Both groups as dots, with the panel title."""
    ax.scatter(*group_a.T, color=style.BLUE, s=30, zorder=3)
    ax.scatter(*group_b.T, color=style.CORAL, s=30, zorder=3)
    ax.set(
        title=title, xticks=[], yticks=[], aspect="equal", xlim=(-1.5, 4.5), ylim=(-1.6, 2.6)
    )


def join(ax: Axes, i: int, j: int, **kw) -> None:
    """A line from point i of group a to point j of group b."""
    ax.plot(*np.c_[group_a[i], group_b[j]], color=style.INK, **kw)


fig, axes = style.subplots(5, height=2.6, width=11)
i, j = np.unravel_index(pair_distance.argmin(), pair_distance.shape)
setup_pair(axes[0], f"single: {pair_distance.min():.1f}")
join(axes[0], i, j, lw=2)
i, j = np.unravel_index(pair_distance.argmax(), pair_distance.shape)
setup_pair(axes[1], f"complete: {pair_distance.max():.1f}")
join(axes[1], i, j, lw=2)
setup_pair(axes[2], f"average: {pair_distance.mean():.1f}")
for i in range(5):
    for j in range(5):
        join(axes[2], i, j, lw=0.6, alpha=0.5)
centre_a, centre_b = group_a.mean(axis=0), group_b.mean(axis=0)
setup_pair(axes[3], f"centroid: {np.linalg.norm(centre_a - centre_b):.1f}")
axes[3].plot(*np.c_[centre_a, centre_b], color=style.INK, lw=2, marker="x", ms=9)
both = np.r_[group_a, group_b]
centre = both.mean(axis=0)
increase = (
    ((both - centre) ** 2).sum()
    - ((group_a - centre_a) ** 2).sum()
    - ((group_b - centre_b) ** 2).sum()
)
setup_pair(axes[4], f"Ward: +{increase:.1f}")
for point in both:
    axes[4].plot(*np.c_[point, centre], color=style.GREY, lw=0.8)
axes[4].plot(*centre, "x", color=style.INK, ms=10, mew=2)
fig.tight_layout()
plt.show()

Figure 7.2.5: Five ways to measure the distance between two groups (blue and red). Single linkage: the closest pair. Complete: the farthest pair. Average: the mean of all 25 pairs. Centroid: between the two group centres. Ward: how much the within-group sum of squares grows if the two groups are joined, the squares to the joint centre (grey cross) minus the squares to each group’s own centre.

The linkage matters most when the groups are not round blobs:

moons, _ = datasets.make_moons(300, noise=0.06, random_state=0)
shapes = np.r_[moons, local.normal([2.6, 1.6], 0.15, (80, 2))]

fig, axes = style.subplots(4, height=2.7, width=11)
for ax, method in zip(axes, ["single", "complete", "average", "ward"]):
    groups = AgglomerativeClustering(3, linkage=method).fit_predict(shapes)
    ax.scatter(*shapes.T, s=6, c=GROUP_COLORS[groups])
    ax.set(title=method, xticks=[], yticks=[])
fig.tight_layout()
plt.show()

Figure 7.2.6: Two interlocking half-moons and a round blob (simulated), cut into three groups by each linkage. Single linkage follows the curved shapes, because it only needs a chain of close neighbours. Ward, complete and average linkage prefer compact, round groups and cut across the moons.

Single linkage pays for that flexibility: a few points between two groups form a bridge, and it joins them. Centroid linkage can even produce inversions, joins lower than earlier joins, which make the dendrogram hard to read.

The clustered heatmap

A common way to show a clustering of a data table is a clustered heatmap: the table as colours, with rows and columns both reordered by their own dendrograms, so that similar rows sit together and similar columns sit together.

row_tree = linkage(scaled, method="ward")
column_tree = linkage(scaled.T, method="average", metric="correlation")

fig = plt.figure(figsize=(9, 7))
grid = fig.add_gridspec(
    2, 3, width_ratios=[1.2, 6, 0.3], height_ratios=[1.2, 6], wspace=0.02, hspace=0.02
)
top = fig.add_subplot(grid[0, 1])
columns = dendrogram(
    column_tree, ax=top, color_threshold=0, above_threshold_color=style.INK, no_labels=True
)
top.set_axis_off()
left = fig.add_subplot(grid[1, 0])
rows = dendrogram(
    row_tree,
    ax=left,
    orientation="left",
    color_threshold=0,
    above_threshold_color=style.INK,
    no_labels=True,
)
left.set_axis_off()
left.invert_yaxis()
heat = fig.add_subplot(grid[1, 1])
heat.imshow(
    scaled[rows["leaves"]][:, columns["leaves"]],
    aspect="auto",
    interpolation="nearest",
    cmap="RdBu_r",
    vmin=-2.5,
    vmax=2.5,
)
names = [
    wine.feature_names[k].replace("_", " ").replace("od280/od315 of diluted wines", "od280")
    for k in columns["leaves"]
]
heat.set_xticks(range(13), names, rotation=60, ha="right", fontsize=8)
heat.set_yticks([])
heat.grid(False)
strip = fig.add_subplot(grid[1, 2])
strip.imshow(
    cultivar[rows["leaves"]][:, None], aspect="auto", cmap=ListedColormap(CULTIVAR_COLORS)
)
strip.set_axis_off()
plt.show()

Figure 7.2.7: The 178 wines (rows) by their 13 standardised measurements (columns), red above the average and blue below. Rows are ordered by a Ward dendrogram of the wines (left), columns by an average linkage dendrogram of the measurements with correlation distance (top). The strip on the right gives the cultivar, which the clustering did not use: blocks of wines with the same profile line up with the cultivars.

7.2.5 Consensus clustering

K-means depends on its random start (Section 7.2.2), and any clustering depends on which cases happened to be in the sample. Consensus clustering turns that instability into information. Cluster many random subsamples of the data, say 80% of the cases each time, with a fresh random start each time. Then, for every pair of cases, count how often they landed in the same group, among the runs where both were sampled:

$$M_{ij} = \frac{\text{runs where } i \text{ and } j \text{ are in the same group}}{\text{runs where both } i \text{ and } j \text{ were sampled}}$$

The consensus matrix M is 1 for pairs that always go together, 0 for pairs that never do. For a good k, almost all pairs are near 0 or 1: the clustering is the same every time. For a poor k, many pairs sit in between, because cases are shuffled between groups from run to run.

KS = (2, 3, 4, 5)


def consensus_matrix(k: int, runs: int = 200) -> np.ndarray:
    """Share of the runs, among those sampling both, in which each pair shared a k-means group."""
    local = np.random.default_rng(k)
    n = len(scaled)
    together, sampled = np.zeros((n, n)), np.zeros((n, n))
    for run in range(runs):
        rows = local.choice(n, int(0.8 * n), replace=False)
        groups = KMeans(k, n_init=1, init="random", random_state=run).fit_predict(scaled[rows])
        sampled[np.ix_(rows, rows)] += 1
        together[np.ix_(rows, rows)] += groups[:, None] == groups[None, :]
    return together / np.maximum(sampled, 1)


def ordered(matrix: np.ndarray) -> np.ndarray:
    """Rows and columns in the order of an average linkage tree on 1 − consensus."""
    distance = (1 - matrix)[np.triu_indices(len(matrix), 1)]
    order = dendrogram(linkage(distance, method="average"), no_plot=True)["leaves"]
    return matrix[np.ix_(order, order)]


matrices = {k: consensus_matrix(k) for k in KS}
pairs = {k: matrices[k][np.triu_indices(len(scaled), 1)] for k in KS}
ambiguous = {k: np.mean((pairs[k] > 0.1) & (pairs[k] < 0.9)) for k in KS}

fig = plt.figure(figsize=(10.5, 6))
grid = fig.add_gridspec(2, 4, height_ratios=[1, 1])
for col, k in enumerate(KS):
    ax = fig.add_subplot(grid[0, col])
    ax.imshow(ordered(matrices[k]), cmap="Blues", vmin=0, vmax=1)
    ax.set(title=f"k = {k}", xticks=[], yticks=[])
ax = fig.add_subplot(grid[1, :2])
for k, color in zip(KS, (style.GREY, style.BLUE, style.GOLD, style.CORAL)):
    ax.step(
        np.sort(pairs[k]), np.linspace(0, 1, len(pairs[k])), color=color, lw=2, label=f"k = {k}"
    )
ax.yaxis.set_major_formatter(style.PERCENT)
ax.legend(fontsize=8)
ax.set(xlabel="consensus of a pair", ylabel="share of pairs at or below")
ax = fig.add_subplot(grid[1, 2:])
ax.bar([str(k) for k in KS], [ambiguous[k] for k in KS], color=style.BLUE, alpha=0.7)
ax.yaxis.set_major_formatter(style.PERCENT)
ax.set(xlabel="k", ylabel="ambiguous pairs (PAC)")
fig.tight_layout()
plt.show()

Figure 7.2.8: Consensus clustering of the wines with k-means, 200 runs on 80% subsamples, one random start each. Top: the consensus matrix for k = 2 to 5, rows and columns ordered by clustering the matrix itself (white: never together, dark: always together). At k = 3 the matrix is three clean blocks; at k = 4 and 5 the extra groups come and go between runs, leaving grey. Bottom left: the distribution (cumulative) of all pair values; a good k has a flat middle. Bottom right: the proportion of ambiguous pairs, with consensus between 0.1 and 0.9 (PAC); lowest at k = 3.

The final groups come from clustering the consensus matrix itself (here, average linkage on 1 − M, cut at k groups). Cases that keep switching groups can be flagged as uncertain rather than forced into one. The same recipe works with any clustering method in place of k-means, and the same idea, “rerun on resampled data and see what survives”, is the bootstrap of Chapter 3.3 applied to groups.

7.2.6 What counts as similar

All of this rests on a distance. With standardised measurements and Euclidean distance, two cases are similar when their values are close. Sometimes the shape of a profile matters more than its level: two shoppers who buy the same mix of products, one ten times as much as the other. Then a correlation-based distance, 1 − r between the two profiles, groups by shape:

SHAPES = np.array([[5, 4, 1, 1, 0.5, 0.5], [0.5, 1, 1, 4, 5, 3]])
profiles, shape_of, level_of = [], [], []
for shape in range(2):
    for level in (1, 1.3, 4, 5.5, 6, 1.5):
        profiles.append(SHAPES[shape] * level + rng.normal(0, 0.3, 6))
        shape_of.append(shape)
        level_of.append(level)
profiles = np.array(profiles)


def plot_profiles(ax: Axes, metric: str, title: str) -> None:
    """Profiles coloured by the two groups found with one distance."""
    groups = (
        fcluster(linkage(profiles, method="average", metric=metric), 2, criterion="maxclust")
        - 1
    )
    for profile, g in zip(profiles, groups):
        ax.plot(range(1, 7), profile, "o-", color=GROUP_COLORS[g], lw=1.2, ms=3, alpha=0.8)
    ax.set(title=title, xlabel="product type", xticks=range(1, 7))


fig, axes = style.subplots(2, height=3.1, width=9.5, sharey=True)
plot_profiles(axes[0], "euclidean", "Euclidean: by amount")
plot_profiles(axes[1], "correlation", "Correlation: by mix")
axes[0].set_ylabel("spending")
fig.tight_layout()
plt.show()

Figure 7.2.9: Twelve simulated shoppers, each a profile of spending on six product types. Left: grouped by Euclidean distance, the shoppers split by how much they spend. Right: grouped by correlation distance, they split by what they buy, whatever the amount.

The choice is about the question, not the data. The same goes for scaling (Chapter 7.1): with raw units the measurement with the largest numbers decides who is “close”.

7.2.7 Groups of any shape: density-based clustering

K-means and Ward’s method look for round groups and put every case in some group. DBSCAN (density-based clustering) instead defines a group as a region where points are packed densely. Two settings, a radius ε and a minimum count m, define three kinds of points:

• a core point has at least m points within distance ε;
• a border point is within ε of a core point but is not core itself;
• everything else is noise, and belongs to no group.

Core points within ε of each other are joined into one group, with their border points. The number of groups is not chosen; it follows from the density.

noisy_shapes = np.r_[shapes, local.uniform([-1.3, -0.8], [3.2, 2.2], (30, 2))]


def plot_groups(ax: Axes, groups: np.ndarray, title: str) -> None:
    """Points coloured by group, noise (label −1) as grey crosses."""
    palette = np.array([plt.cm.tab10(g % 10) for g in range(max(groups.max() + 1, 1))])
    member = groups >= 0
    ax.scatter(*noisy_shapes[member].T, s=6, c=palette[groups[member]])
    ax.scatter(*noisy_shapes[~member].T, s=14, marker="x", color=style.GREY)
    ax.set(title=title, xticks=[], yticks=[])


fig, axes = style.subplots(4, height=2.7, width=11)
plot_groups(
    axes[0], KMeans(3, n_init=20, random_state=0).fit_predict(noisy_shapes), "k-means, k = 3"
)
for ax, radius in zip(axes[1:], (0.15, 0.06, 0.5)):
    groups = DBSCAN(eps=radius, min_samples=5).fit_predict(noisy_shapes)
    found = len(set(groups) - {-1})
    plot_groups(ax, groups, f"DBSCAN, ε = {radius}: {found} groups")
fig.tight_layout()
plt.show()

Figure 7.2.10: Two half-moons and a blob with 30 scattered noise points (simulated). K-means with k = 3 cuts across the moons and forces the noise into groups. DBSCAN (m = 5) with a suitable radius ε = 0.15 finds the three shapes and leaves the scattered points as noise (grey crosses). Too small a radius breaks the shapes into many pieces and calls much of the data noise; too large a radius merges the two moons into one group.

A useful way to pick ε: for every point, find the distance to its m-th nearest neighbour, sort these distances, and look for the bend where they start to rise steeply. HDBSCAN tries all radii at once and keeps the groups that persist, which also copes with groups of different densities. Spectral clustering is another route to non-round groups: it builds a graph connecting each point to its near neighbours and cuts the graph where it has the fewest connections.

7.2.8 Soft clustering

K-means and hierarchical clustering put every case firmly in one group. A Gaussian mixture (Chapter 2.4) treats the groups as overlapping normal clouds and gives each case a probability of belonging to each group, which shows where the groups are uncertain:

mixture = GaussianMixture(3, n_init=5, random_state=0).fit(plane)
probabilities = mixture.predict_proba(plane)

fig, ax = plt.subplots(figsize=(6.5, 4))
ax.scatter(
    plane[:, 0],
    plane[:, 1],
    c=GROUP_COLORS[probabilities.argmax(axis=1)],
    s=4 + 60 * (probabilities.max(axis=1) - 0.33) ** 2 * 4,
    alpha=0.7,
)
ax.set(xlabel="component 1", ylabel="component 2")
plt.show()

Figure 7.2.11: A mixture of three normal clouds fitted to the wines (on the first two principal components). Each wine is coloured by its most likely group; the size of the dot shows how sure the model is. Wines between the clouds get split probabilities and small dots.

7.2.9 Practical issues

• Stability: rerun the clustering on bootstrap resamples (Chapter 3.3) or subsets, as consensus clustering does; groups that keep appearing are more believable than groups that come and go.
• Outliers pull centres and can form groups of their own; consider removing or down-weighting them, or use methods that allow cases to belong to no group.
• Validation: groups are a hypothesis. Check them against outside information that the clustering did not use, as the cultivars did here.
• Clustering always returns groups, even in data with none.

7.2.10 In one picture

• K-means minimises the within-cluster variation $W$ by alternating “nearest centre” and “move the centre”; run it from many starts.
• Choose k with the elbow, the silhouette $(b - a)/\max(a, b)$, and judgement.
• Hierarchical clustering joins the closest groups step by step; the linkage (single, complete, average, centroid, Ward) defines “closest”. Single linkage follows chains and odd shapes; Ward finds compact groups. A clustered heatmap orders rows and columns by their trees.
• Consensus clustering reruns the clustering on subsamples; the share of runs in which two cases share a group shows which groups, and which k, are stable.
• The distance decides what “similar” means: Euclidean for levels, correlation for shapes.
• DBSCAN finds dense regions of any shape and leaves sparse points as noise.
• Gaussian mixtures give soft memberships.
• Check stability and meaning; unsupervised methods always find something.

Next: finding the cases that do not fit.