Cluster Evaluation and Profiling

Your k-Means run on 40,000 customers returns a silhouette of 0.62 and a Calinski-Harabasz score of 996. You paste both into Slack, and the head of CRM asks the only question that matters: "Great — what do I do differently on Monday?" If your answer is "the clusters are compact and well separated," the project dies right there. Internal metrics prove the algorithm found structure; they never prove the structure is worth money. A defensible evaluation climbs three rungs — internal, external, business — and only the third one gets budget.

Rung one: the internal arithmetic. Interview panels like to ask you to compute these scores from cluster statistics, so know the two formulas cold.

The silhouette coefficient for a point i is s(i) = (b(i) − a(i)) / max(a(i), b(i)), where a(i) is its mean distance to its own cluster and b(i) the mean distance to the nearest other cluster. A point with a(i) = 1.6 and b(i) = 3.2 scores (3.2 − 1.6)/3.2 = 0.50. The dataset score is the mean over all points: above 0.5 signals real structure, above 0.7 strong separation, near 0 means the algorithm drew boundaries through soup, and negative values mean points are closer to a neighbor cluster than their own. It is undefined at k = 1 and O(n²) in memory, so on 40,000 rows you compute it on a 5,000-row sample and say so in the writeup.

The Calinski-Harabasz score has the shape of an F-test: CH = (SSB / (k − 1)) / (SSW / (n − k)), where SSW is exactly the k-Means inertia, Σ_c Σ_{x∈c} ‖x − μ_c‖², and SSB is the between-cluster dispersion, Σ_c n_c‖μ_c − μ‖², with μ the global mean. Plug in real numbers: n = 1000, k = 4, SSW = 8200, SSB = 24600 gives CH = (24600/3) / (8200/996) = 8200 / 8.23 ≈ 996. Higher is better, but the absolute value means nothing alone — it is comparable only across k values on the same scaled data, because both dispersions shrink and grow with the dataset.

from sklearn.metrics import silhouette_score, calinski_harabasz_score

sample_idx = rng.choice(len(X_scaled), 5000, replace=False)
print(silhouette_score(X_scaled[sample_idx], labels[sample_idx]))  # O(n^2), so sample
print(calinski_harabasz_score(X_scaled, labels))                   # cheap enough for all rows

Rung two: external labels. If a prior segmentation exists — a manual taxonomy, last year's clusters, domain annotations — score agreement with adjusted Rand index or normalized mutual information. An ARI near 0 means your clusters carry no information about the old ones, which is a finding rather than a failure: either the old taxonomy was arbitrary or your features tell a different story. Either way you now know whether you replaced a process or reproduced it.

Rung three: the business test. This is where clustering projects live or die, and it has two flavors:

The rule that prevents self-deception: choose the business metric before you cluster. Churn per cluster with overlapping confidence intervals means your features do not carry churn signal — no amount of re-running k-Means fixes that; you need different features. And the experiment no stakeholder can argue with is the randomized one: serve the "likely to lapse" cluster a win-back offer, hold out a random 10% of that cluster, and report the conversion delta with a p-value.

One mechanics note before you score anything: DBSCAN output contains the −1 noise label, and treating noise as an ordinary cluster inflates or deflates every metric above depending on its size. Decide explicitly whether noise is excluded or reported as its own segment — the mechanics are covered in DBSCAN Density Clustering.

The centroid table says C2 has recency mean 24.8, orders mean 6.3, average-order-value mean 182.4. The CRM team's response is "so what?" Profiling is the translation layer: it converts cluster centers into named segments with sizes, spreads, and behaviors that map to distinct actions. Done well, the profile table is the deliverable — the clustering was preprocessing.

Pass one: centers in original units. If you clustered on standardized features — and you did, after k-Means Clustering drilled that in — each centroid lives in z-space. Invert the scaler before anyone reads it: a mean recency of −0.7 z means nothing to a marketer; 24 days means "about three weeks." Report medians alongside means, because customer features are almost always right-skewed: the mean revenue of a cluster with two whale accounts describes nobody in it, while the median describes almost everyone.

Pass two: the distinguishing-feature table. For each cluster, compute how far its feature mean sits from the global mean, in standard deviations: z_c,f = (mean_c,f − global_mean_f) / std_f. The feature with the largest |z| is what the cluster is for. A real profile from a retail dataset:

Four rows, four different actions — that is the acceptance test for a profile: if two rows imply the same action for the customer, merge the clusters or raise k. Name segments from the distinguishing column, not from hope, and hold off on celebratory names like "champions" until the churn/LTV separation test from the previous chapter passes.

import pandas as pd
import matplotlib.pyplot as plt

feats = ['recency_days', 'orders_12m', 'avg_order_value']
df['cluster'] = labels

med = df.groupby('cluster')[feats].median().round(1)
z = ((df.groupby('cluster')[feats].mean() - df[feats].mean())
     / df[feats].std()).round(2)
sizes = df.groupby('cluster').size().rename('n')
print(pd.concat([sizes, med, z.add_suffix('_z')], axis=1))

fig, axes = plt.subplots(1, len(feats), figsize=(14, 4))
for ax, f in zip(axes, feats):
    df.boxplot(column=f, by='cluster', ax=ax)
    ax.set_title(f)
plt.tight_layout()

Pass three: the boxplots, because means hide spread. A cluster with mean AOV $95 could be a tight $90–100 group you can sell to with one offer, or a bimodal pile of $30 and $400 orders that is really two segments the algorithm failed to split. The boxplot shows median, IQR, whiskers, and outliers per cluster per feature in one glance. Read three signals off it: a tall box means the segment is loose and targeting it will be noisy; two clusters with boxes identical in position and spread mean that feature separates nothing (check your scaling — on unscaled features the biggest-magnitude column dominates both the clustering and the profile); and a bimodal box is a split waiting to happen.

Finally, profile the noise set. If DBSCAN's −1 group carries 4× the median order value of every cluster, "noise" was the most interesting segment all along — often fraud rings or bulk buyers, which is exactly the finding worth presenting to a fraud team.

You rerun the notebook with random_state=7 instead of 42, and the "Lapsed low-value" segment loses 40% of its members to other clusters. Was the segmentation real, or a coin flip with extra steps? k-Means guarantees convergence, not stability: n_init=10 keeps the best of ten restarts, but that is still one draw from a distribution of plausible partitions. Stability checks measure how much of your answer survives perturbation, and they separate segments you can build on from artifacts you cannot.

Test one: the seed sweep. Fit k-Means 20 times with different random_state values and score each run against a reference fit with the adjusted Rand index — pair-counting agreement between two labelings, corrected for the chance agreement that cluster sizes alone would produce. ARI is 1.0 for identical partitions, ~0 for chance, and can go negative.

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

ref = KMeans(n_clusters=4, n_init=10, random_state=42).fit_predict(X_scaled)
aris = [adjusted_rand_score(ref, KMeans(n_clusters=4, n_init=10,
                            random_state=s).fit_predict(X_scaled))
        for s in range(20)]
print(f"ARI mean={np.mean(aris):.2f} std={np.std(aris):.2f}")

Test two: subsample stability. Seeds change the optimizer; subsampling changes the data — the harsher test. Draw 80% of rows without replacement, cluster the subset, and compute ARI against the reference labels restricted to the sampled rows:

rng = np.random.default_rng(42)
aris = []
for _ in range(20):
    idx = rng.choice(len(X_scaled), int(0.8 * len(X_scaled)), replace=False)
    sub = KMeans(n_clusters=4, n_init=10, random_state=42).fit_predict(X_scaled[idx])
    aris.append(adjusted_rand_score(ref[idx], sub))

Read the numbers against rough thresholds: ARI ≥ 0.8 — the partition is reproducible; 0.6–0.8 — usable, but report the unstable clusters explicitly; below 0.6 — you are segmenting noise. For per-cluster detail, match each cluster to its best-overlap counterpart across runs and compute Jaccard = |A∩B| / |A∪B|. One cluster at Jaccard 0.45 while the rest sit at 0.9 is a localized failure, and the fix is usually local: merge it, or accept that this subpopulation has no stable boundary at this k.

Method-specific wrinkles: with DBSCAN, compute ARI over non-noise points only and track the −1 fraction separately, since a wandering noise fraction invalidates every other number; with hierarchical clustering, re-cut the tree at the same k — the algorithm is deterministic, so any instability you see comes from the data, not the optimizer.

The last check is the one production teams forget: temporal holdout. Fit centroids on week 1, assign week 2's customers by nearest centroid, and compare per-cluster proportions and profiles. Segments that evaporate out of sample were tuned to last week's noise.

When an interviewer asks "your clusters changed between runs — how do you know they're not noise?", the protocol answers itself: seed sweep, subsample ARI against a reference, per-cluster Jaccard to localize the failure, then either lower k or report confidence honestly. And if a segmentation survives all of that in 40 dimensions and the scores still look worse than they should, the reason is Curse of Dimensionality — your next topic.