Are your clusters real?

The slide says the customer base splits into four segments. There is a scatter plot with four colours, a table of averages per segment, and names: Bargain Hunters, Loyalists, Dormant, Whales. The room accepts it, because the picture shows four groups and the groups have names.

Nothing on that slide is evidence that four groups exist. The algorithm was asked for four and it returned four. It would have returned four from a table of random numbers, with a scatter plot that looks much the same.

The algorithm never refuses

k-means partitions. That is the whole contract. Given k, it draws straight boundaries that minimise within-cluster spread, and it does that on structured data and on noise with equal confidence. So does hierarchical clustering cut at a height, and so does any other method that assigns every point to a group.

Run it on uniform noise and look at what comes back:

import numpy as np
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

rng = np.random.default_rng(0)
X = rng.uniform(size=(2000, 2))              # no structure whatsoever

labels = KMeans(n_clusters=4, n_init=10, random_state=0).fit_predict(X)
print(f"{silhouette_score(X, labels):.3f}")

0.408

0.408 on data with nothing in it, and the plot is four tidy regions with visible borders. Neither number nor picture is embarrassing on its own. That is the problem: they are the same output real structure would give, so seeing them tells you nothing about which case you are in.

Everything below is an attempt to answer one question the algorithm cannot: would this partition have appeared anyway?

Compare against data with the structure removed

The direct way is to build a null. Generate data with no clusters but the same shape as yours, cluster it the same way, and see how much better the real data scores.

Two ways to build the reference set. Sample uniformly inside the bounding box of each feature, which is what the gap statistic does. Or shuffle each column independently, which keeps every feature’s real distribution — the skew, the long tail, the spikes at round numbers — and destroys only the relationships between columns. The second is usually the fairer test, because a heavy-tailed column can produce apparent clusters by itself.

def inertia(X, k, seed=0):
    return KMeans(n_clusters=k, n_init=10, random_state=seed).fit(X).inertia_

def shuffled_reference(X, rng):
    Z = X.copy()
    for j in range(Z.shape[1]):
        rng.shuffle(Z[:, j])
    return Z

k = 4
real = np.log(inertia(X, k))
null = [np.log(inertia(shuffled_reference(X, rng), k)) for _ in range(20)]
print(f"{np.mean(null) - real:.4f} {np.std(null):.4f}")   # gap, and its noise

-0.0008 0.0013

A gap near zero means the partition is doing no better on your data than on data whose structure has been deleted. Run it across a range of k and the curve is more informative than any single value: real structure shows a gap that opens up at the right k. Noise shows a gap that wanders around zero.

This is the idea behind the gap statistic, stated plainly rather than derived. The published version is more careful about the reference distribution and the correction term; the plain version is enough to stop a bad slide.

Now run the same sweep on real data. The features below come from the public dataset: four counts per subscription, of what the customer did in their first thirty days.

import pandas as pd
from sklearn.preprocessing import StandardScaler

subs = pd.read_csv("https://dataacademy.ai/data/subscriptions.csv",
                   parse_dates=["signed_up_at"])
events = pd.read_csv("https://dataacademy.ai/data/events.csv",
                     parse_dates=["occurred_at"])

ev = events.merge(subs[["subscription_id", "customer_id", "signed_up_at"]],
                  on="customer_id")
month1 = ev[(ev.occurred_at - ev.signed_up_at).dt.days.between(0, 29)]
counts = month1.pivot_table(index="subscription_id", columns="event_type",
                            values="event_id", aggfunc="count", fill_value=0)
C = StandardScaler().fit_transform(counts)

rng = np.random.default_rng(0)
for k in range(2, 8):
    real = np.log(inertia(C, k))
    null = [np.log(inertia(shuffled_reference(C, rng), k)) for _ in range(20)]
    print(f"k={k}  gap {np.mean(null) - real:.3f}  sd {np.std(null):.3f}")

k=2  gap 0.467  sd 0.004
k=3  gap 0.596  sd 0.005
k=4  gap 0.552  sd 0.006
k=5  gap 0.447  sd 0.015
k=6  gap 0.472  sd 0.010
k=7  gap 0.483  sd 0.023

Every gap is large and hundreds of standard deviations from zero. None of them is a peak. The curve rises to k=3, falls, and then drifts sideways — there is no k at which the partition suddenly starts earning its keep, only a constant advantage over the shuffled version at every k tried. That pattern has a cause, and it is not clusters: shuffling each column independently destroys the correlations between the four counts, and those four counts are strongly correlated. Any partition of correlated data beats the same partition of decorrelated data. The test measured the correlation, not a gap.

Stability is the test that survives contact with real data

If the four segments are real, they should not depend on which customers happened to be in the extract. Cluster a random 80% of the rows, do it again with a different 80%, and ask whether the same points keep landing together.

Pairwise co-assignment is the practical form. For every pair of points, count the runs where both were sampled, and the runs where both were sampled and given the same label.

def coassignment(X, k, n_runs=50, frac=0.8, seed=0):
    rng = np.random.default_rng(seed)
    n = len(X)
    together = np.zeros((n, n))
    seen = np.zeros((n, n))
    for run in range(n_runs):
        idx = rng.choice(n, size=int(frac * n), replace=False)
        labels = KMeans(n_clusters=k, n_init=10,
                        random_state=run).fit_predict(X[idx])
        same = (labels[:, None] == labels[None, :]).astype(float)
        together[np.ix_(idx, idx)] += same
        seen[np.ix_(idx, idx)] += 1
    return np.divide(together, seen, out=np.zeros_like(together),
                     where=seen > 0)

That matrix is n by n, so it is a few thousand rows before memory bites; subsample first if needed. Read it per cluster: take the points a full-data run assigned to segment 2 and look at their average co-assignment. High means those customers stay together whatever the sample. Around chance means segment 2 is a boundary the algorithm redraws every time it is asked.

On the four counts from before, every segment passes:

from sklearn.decomposition import PCA

M = coassignment(C, 4)
full = KMeans(n_clusters=4, n_init=10, random_state=0).fit_predict(C)
for c in range(4):
    m = full == c
    block = M[np.ix_(m, m)]
    iu = np.triu_indices(m.sum(), 1)
    print(f"cluster {c}  n {m.sum():4d}  mean co-assignment {block[iu].mean():.3f}")

print(PCA().fit(C).explained_variance_ratio_.round(3))

cluster 0  n   95  mean co-assignment 0.891
cluster 1  n  762  mean co-assignment 0.977
cluster 2  n   25  mean co-assignment 0.878
cluster 3  n  316  mean co-assignment 0.917
[0.766 0.14  0.074 0.02 ]

Four segments, co-assignment between 0.878 and 0.977, a large gap against the null — a slide with those numbers on it would clear any review. The last line is why it should not. One direction holds 76.6% of the variance in these four columns, because customers who view a lot also log in, export and buy a lot. The data is a single activity gradient, and k-means cut it into four bands at the same three places every time.

Stability is not proof of structure — a strong gradient with no gaps in it can be cut in the same place repeatedly — but instability is proof of its absence, and it is the cheapest check available.

Internal indices score agreement with an assumption

Silhouette — Rousseeuw’s 1987 index — and Davies-Bouldin measure whether clusters are compact and well separated. Compact and well separated is exactly what k-means is built to produce, so a good score partly confirms that k-means behaved like k-means.

The failure shows up on data where the right answer is not round. Two interlocking crescents, correctly separated by a density-based method, score badly on silhouette, because each crescent is long and parts of it sit closer to the other crescent than to its own far end. Silhouette prefers the wrong answer there.

Before any of this, scaling and the distance metric have already decided most of the outcome. Standardising declares every feature equally important, and one unscaled column in a large unit will otherwise be the only column in the model at all — what k-means is really doing covers why that follows directly from squared distance. Change the scaling, change the metric, and the segments change. Report which choices were made, because they are part of the finding.

What a cluster is allowed to be used for

Report the checks, not just the partition: how k was chosen, the gap against a null, co-assignment per cluster, and whether cluster profiles differ in ways anyone can act on. Two clusters that differ only in a variable nobody controls are one cluster with extra steps.

Then separate two things that get the same word. A discovered cluster is a claim about the data and needs the evidence above. A business segment is a decision — marketing can run four campaigns, so there are four segments — and needs no evidence at all, only a sensible rule and a name.

Both are legitimate. Trouble starts when a decision is presented as a discovery, because discoveries are expected to hold. “We split the base into four groups by spend and recency, and here is what each one costs to serve” is honest and useful. “The data revealed four natural customer types” is a claim that most customer bases, being continuous, cannot support.