Predictive Modeling and Machine Learning Fundamentals Questions

Direct answer

Aggregate the raw events into one row per user of behaviour summaries (how often, how intensely, how broadly, how recently) over a shared observation window, put the columns on a comparable scale, cluster with a method that fits the shape of the data (k-means as a first pass), pick the number of segments with several measures plus a stability check and a usefulness test, profile each segment in the original units, give each a plain name, and validate the segments on something the clustering did not use. Segments are only useful if a product team can say what it would do differently for each one.

Terms

Segment: a group of users with similar behaviour. k-means: an algorithm that puts k centre points in the data, assigns each user to the nearest centre, moves each centre to the mean of its users, and repeats. Inertia: the total squared distance from users to their centres, lower means tighter clusters but it always falls as k grows. Silhouette: for each user, how much closer they are to their own cluster than to the next closest one, from -1 to 1. ARI (adjusted Rand index): agreement between two clusterings, 1 for identical and about 0 for chance.

From raw events to features

1. Decide the unit and window. One row per user, a fixed window (for example 30 days), and an agreed set of users (new accounts with too little history are handled separately).
2. Clean events. Remove bots and internal accounts, and keep users with no events (they are a real segment).
3. Summaries a person can read. Active days, events per active day, distinct features used, days since last activity.
4. Transform and scale. Counts are skewed, so take a log, then standardize so no column dominates only because of its units. Without scaling, k-means measures distance mostly in the largest-unit column.
5. Cluster, choose k, check, profile, name, validate.

Worked example

The simulated clickstream has 1,800 users over 30 days generated from three hidden habits (casual, steady, power). The hidden habit is never used for clustering; it is printed at the end only to check the demo, which is possible because the data is simulated. Real data has no such answer key, so the substitute is outcomes the clustering did not use.

import numpy as np
import pandas as pd
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score, adjusted_rand_score

rng = np.random.default_rng(0)
# ---- 1. raw events: one row per (user, day, feature used); 30 days of data, 1,800 users from 3 hidden habits ----
N, DAYS = 1800, 30
habit = rng.choice(3, N, p=[0.5, 0.3, 0.2])                       # 0 = casual, 1 = steady, 2 = power (hidden; never used below)
active_prob = np.array([0.08, 0.35, 0.80])[habit]
events_per_day = np.array([1.5, 3.0, 9.0])[habit]
n_features = np.array([2, 5, 12])[habit]
rows = []
for u in range(N):
    days = np.nonzero(rng.random(DAYS) < active_prob[u])[0]
    for dd in days:
        for _ in range(max(1, rng.poisson(events_per_day[u]))):
            rows.append((u, dd, rng.integers(0, n_features[u])))
ev = pd.DataFrame(rows, columns=["user", "day", "feature"])
print("events:", len(ev), "users with at least one event:", ev.user.nunique())

# ---- 2. one row per user: behaviour summaries, with a shared observation window ----
g = ev.groupby("user")
f = pd.DataFrame({
    "active_days": g.day.nunique(),
    "events_per_active_day": g.size() / g.day.nunique(),
    "distinct_features": g.feature.nunique(),
    "days_since_last": DAYS - 1 - g.day.max(),
}).reindex(range(N))
f["active_days"] = f.active_days.fillna(0)
f = f.fillna({"events_per_active_day": 0, "distinct_features": 0, "days_since_last": DAYS})
X = StandardScaler().fit_transform(np.log1p(f.values))              # log tames the long right tail, then put columns on one scale

# ---- 3. choose k with several views, none of which uses the hidden habit ----
print("\nk  silhouette  inertia   stability (mean ARI between 5 reseeded runs of 10 starts each)")
for k in range(2, 7):
    runs = [KMeans(k, n_init=10, random_state=s).fit(X).labels_ for s in range(5)]
    ari = np.mean([adjusted_rand_score(runs[i], runs[j]) for i in range(5) for j in range(i + 1, 5)])
    km = KMeans(k, n_init=10, random_state=0).fit(X)
    print(f"{k}  {silhouette_score(X, km.labels_):.3f}       {km.inertia_:9.1f}  {ari:.3f}")

# ---- 4. profile the chosen solution in the original units ----
km = KMeans(3, n_init=10, random_state=0).fit(X)
f["segment"] = km.labels_
order = f.groupby("segment").active_days.mean().sort_values().index          # label segments by activity so the numbering is not arbitrary
f["segment"] = f.segment.map({old: new for new, old in enumerate(order)})
prof = f.groupby("segment").agg(users=("active_days", "size"), active_days=("active_days", "mean"),
                                events_per_day=("events_per_active_day", "mean"),
                                features=("distinct_features", "mean"), days_since_last=("days_since_last", "mean")).round(1)
print("\n", prof.to_string())
print("\nshare of each hidden habit inside each segment (used only to check the demo):")
print(pd.crosstab(f.segment, pd.Series(habit, name="hidden_habit"), normalize="index").round(2).to_string())

events: 103727 users with at least one event: 1750

k  silhouette  inertia   stability (mean ARI between 5 reseeded runs of 10 starts each)
2  0.530          2644.5  1.000
3  0.536          1516.9  1.000
4  0.502          1067.6  1.000
5  0.532           787.9  0.992
6  0.507           637.7  0.995

          users  active_days  events_per_day  features  days_since_last
segment                                                               
0          778          2.2             1.6       1.6             11.5
1          640          9.4             2.9       4.4              1.6
2          382         24.1             9.0      12.0              0.2

share of each hidden habit inside each segment (used only to check the demo):
hidden_habit     0     1    2
segment                      
0             0.99  0.01  0.0
1             0.18  0.82  0.0
2             0.00  0.00  1.0

Reading the output:

• 103,727 events, and 1,750 of the 1,800 users appear in them; the other 50 have no events, so the code gives them 0 for active days, events per active day and distinct features, and 30 (the whole window) for days since last activity, because they were last seen at least 30 days ago. That last value is deliberate: a zero there would wrongly say they were active today. These 50 land in the least active segment.
• Choosing k. Silhouette is nearly flat from k = 2 to 5 (0.530, 0.536, 0.502, 0.532), so the 0.006 gap between k = 2 and k = 3 is noise, not a winner. Inertia keeps falling with every k, as it must, and shows no sharp elbow: it drops 43% from k = 2 to 3, 30% from 3 to 4 and 26% from 4 to 5. Stability (mean ARI between five reseeded runs, each the best of 10 starts) is between 0.992 and 1.000 for every k from 2 to 6, so each of these segmentations reproduces and stability does not separate them either. Its job is to warn when a candidate k is not reproducible, and none is a problem here. A single start per run would add a different kind of noise, the luck of the starting centroids, which is why each run uses 10 starts. When the three views agree that several values of k are acceptable, the choice is made on use: k = 3 gives three profiles that differ on every column and suggest three different actions, whereas k = 2 can only split the three profiles into two groups, so two of them share a segment. The demo data was generated from three habits, so 3 is the expected answer here; on real data the same reading applies, and the usefulness test and the outcome validation below carry the decision. What would change it is a product need for finer groups, accepting that each extra group is less distinct.
• Profile in original units. Segment 0 (778 users): about 2.2 active days, 1.6 events per active day, 1.6 features, last seen 11.5 days ago. Segment 1 (640 users): 9.4 active days, 2.9 events per day, 4.4 features, last seen 1.6 days ago. Segment 2 (382 users): 24.1 active days, 9.0 events per day, 12.0 features, last seen 0.2 days ago. Natural names: occasional visitors, regular users, power users.
• Check against the answer key. Segment 0 is 99% casual-habit users, segment 1 is 82% steady, segment 2 is 100% power. Segment numbers are arbitrary, so the code orders them by average active days to keep them reproducible.

Making segments actionable

Give each segment a rule a system can apply (the profile above translated into thresholds, or a saved model assigning new users), a name, a size, and one action: for example onboarding nudges for occasional visitors, feature discovery prompts for regular users, early access or feedback invitations for power users. Validate on outcomes the clustering did not use, such as 90-day retention, and re-run the clustering on a later window to see whether the same segments reappear. Refresh on a schedule so assignments do not drift silently.

When k-means is not the right tool

K-means assumes roughly round, similar-sized groups and is pulled by outliers. Density-based clustering (DBSCAN) finds irregular shapes and labels noise, hierarchical clustering gives a tree of merges to inspect, and Gaussian mixtures give soft memberships. In practice k-means (often after trying hierarchical clustering on a sample) is the usual first pass, and DBSCAN and Gaussian mixtures are second-line tools that you reach for only if profiles look smeared or a few extreme users dominate.

Pitfalls

• Clustering raw counts: heavy users dominate the distance.
• Choosing k with only one score, or from the elbow of inertia alone.
• Presenting segments without an action, a size, or a validation outcome.
• Treating cluster numbers as meaningful or stable across re-runs.

Direct answer

K-means is fast, simple and gives each segment a centroid (its mean customer) that is easy to describe, but it assumes clusters are compact, roundish, of similar spread, and it needs $k$ (the number of clusters) up front, depends on the starting centroids and is pulled by outliers and by feature units. For a customer segmentation I would start with k-means on standardised, sensibly transformed features as the baseline, because stakeholders can read a centroid. I would switch method only when diagnostics show the shapes violate its assumptions.

Strengths and weaknesses

Evidence on five datasets

A single outlier shows the mechanism first: the mean of 1, 2, 3, 4, 5 is 3, and adding one value of 100 moves it to 19.17. A centroid is a mean, so far points drag it in the same way. The program then generates five labelled datasets: three overlapping round blobs of equal spread, two crescents, three blobs on a line with one much wider than the others, and two blobs plus 30 outlier points in two ways (scattered uniformly over a wide box, or packed into one far spot). It scores k-means, a Gaussian mixture (GMM: each cluster is a Gaussian bell shape with its own centre and its own covariance, meaning its own spread and tilt, set here by covariance_type="full"; it gives soft membership, a probability of belonging to each cluster rather than a single label), and DBSCAN (density-based: a point with at least min_samples points (itself included) inside radius eps starts a dense region, regions grow into clusters, and sparse points are noise, labelled -1) against the planted labels with the adjusted Rand index (ARI: 1 is a perfect match, about 0 a random one). The outlier rows are scored on the 600 real points only, so a perfect method can reach 1.0. Features are standardised first. DBSCAN's eps is the best of six values, chosen with the labels, so its column is an optimistic upper bound that a real project cannot compute; the last column is how many clusters DBSCAN formed at that eps. Without labels, a common way to pick eps is a k-distance plot: sort every point's distance to its min_samples-th nearest neighbour (counting the point itself, as in the rule above) and set eps near the bend where distances start to climb steeply. After the table the program prints how many clusters DBSCAN formed at each eps on the round blobs, where the real points of the scattered-outlier row went, and a redraw of both outlier rows with 20 other seeds.

import numpy as np
from sklearn.cluster import KMeans, DBSCAN
from sklearn.datasets import make_blobs, make_moons
from sklearn.metrics import adjusted_rand_score
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler

# One outlier drags a mean: five customers, then the same five plus one extreme value
five = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
print("mean of 1..5:", five.mean(), "  mean after adding a value of 100:", round(float(np.append(five, 100).mean()), 2))

rng = np.random.default_rng(0)
cases = {}
cases["round blobs, equal size"] = make_blobs(600, centers=3, cluster_std=1.0, random_state=0)
cases["two moons"] = make_moons(600, noise=0.07, random_state=0)
Xa, ya = make_blobs(900, centers=[[0, 0], [5, 0], [10, 0]], cluster_std=[0.6, 0.6, 2.2], random_state=0)
cases["unequal spread"] = (Xa, ya)
Xb, yb = make_blobs(600, centers=[[0, 0], [8, 0]], cluster_std=1.0, random_state=0)
# Outlier rows are scored on the 600 real points only, so a perfect method can reach 1.0
cases["+30 outliers, scattered"] = (np.vstack([Xb, rng.uniform(-30, 40, (30, 2))]), yb)
cases["+30 outliers, one far spot"] = (np.vstack([Xb, rng.normal([40, 40], 0.5, (30, 2))]), yb)
grid = [0.1, 0.15, 0.2, 0.3, 0.4, 0.5]   # DBSCAN eps values tried, after standardising; min_samples=5

blob_dbscan, where = {}, []
print(f"{'dataset':28s} {'k-means':>8s} {'GMM':>8s} {'DBSCAN best':>12s} {'at eps':>7s} {'groups':>7s}")
for name, (X, y) in cases.items():
    k = len(set(y))
    real = np.arange(len(y))                 # the points that have a planted label
    Z = StandardScaler().fit_transform(X)
    km = KMeans(k, n_init=10, random_state=0).fit_predict(Z)
    gm = GaussianMixture(k, covariance_type="full", n_init=5, random_state=0).fit_predict(Z)
    # ARI counts DBSCAN's noise label (-1) as one more group. "Best over the grid" uses the labels,
    # so it is an optimistic upper bound that nobody has in practice.
    labs = {e: DBSCAN(eps=e, min_samples=5).fit_predict(Z) for e in grid}
    db = {e: adjusted_rand_score(y, labs[e][real]) for e in grid}
    e_best = max(db, key=db.get)
    if name.startswith("round"):
        blob_dbscan = {e: len(set(l) - {-1}) for e, l in labs.items()}
    if name.endswith("scattered"):
        for who, lab in (("k-means", km), ("GMM", gm)):
            where.append(f"scattered outliers, {who}: real points per cluster {np.bincount(lab[real], minlength=k).tolist()}, "
                         f"the 30 added points per cluster {np.bincount(lab[len(y):], minlength=k).tolist()}")
    n_groups = len(set(labs[e_best]) - {-1})
    print(f"{name:28s} {adjusted_rand_score(y, km[real]):8.3f} {adjusted_rand_score(y, gm[real]):8.3f} "
          f"{db[e_best]:12.3f} {e_best:7.2f} {n_groups:7d}")

print()
print("round blobs, DBSCAN clusters by eps:", blob_dbscan)
print(*where, sep="\n")

# Are the two outlier rows typical? Redraw the blobs and the 30 added points with 20 other seeds.
for label, draw in (("scattered", lambda r: r.uniform(-30, 40, (30, 2))),
                    ("one far spot", lambda r: r.normal([40, 40], 0.5, (30, 2)))):
    km_ok = km_one = gm_ok = 0
    for s in range(1, 21):
        r = np.random.default_rng(s)
        Xs, ys = make_blobs(600, centers=[[0, 0], [8, 0]], cluster_std=1.0, random_state=s)
        Zs = StandardScaler().fit_transform(np.vstack([Xs, draw(r)]))
        kl = KMeans(2, n_init=10, random_state=0).fit_predict(Zs)[:600]
        gl = GaussianMixture(2, covariance_type="full", n_init=5, random_state=0).fit_predict(Zs)[:600]
        km_ok += adjusted_rand_score(ys, kl) > 0.9
        km_one += np.bincount(kl, minlength=2).max() == 600     # all 600 real points in one cluster
        gm_ok += adjusted_rand_score(ys, gl) > 0.9
    print(f"{label:12s}, 20 redraws: k-means separated the real groups in {km_ok}, "
          f"put all 600 real points in one cluster in {km_one}; GMM separated them in {gm_ok}")

Output:

mean of 1..5: 3.0   mean after adding a value of 100: 19.17
dataset                       k-means      GMM  DBSCAN best  at eps  groups
round blobs, equal size         0.794    0.800        0.353    0.20       8
two moons                       0.480    0.494        1.000    0.20       2
unequal spread                  0.404    0.967        0.802    0.10       6
+30 outliers, scattered         1.000    0.000        1.000    0.30       2
+30 outliers, one far spot      0.000    0.000        1.000    0.20       3

round blobs, DBSCAN clusters by eps: {0.1: 24, 0.15: 16, 0.2: 8, 0.3: 1, 0.4: 1, 0.5: 1}
scattered outliers, k-means: real points per cluster [300, 300], the 30 added points per cluster [19, 11]
scattered outliers, GMM: real points per cluster [600, 0], the 30 added points per cluster [2, 28]
scattered   , 20 redraws: k-means separated the real groups in 7, put all 600 real points in one cluster in 13; GMM separated them in 0
one far spot, 20 redraws: k-means separated the real groups in 0, put all 600 real points in one cluster in 20; GMM separated them in 0

How to read it:

• Round blobs: k-means and the GMM tie (0.794 and 0.800; the blobs overlap, so neither reaches 1). DBSCAN's best is 0.353 because the overlapping blobs have no gap for it to find: at small eps it shatters them into many small dense pieces (24, 16 and 8 clusters at 0.10, 0.15 and 0.20, where 3 are planted, as the line under the table shows), and at eps 0.3 and above it merges everything into one cluster (ARI 0.000). The assumption fits k-means, so the simple method is enough and the extra machinery buys nothing.
• Crescents: k-means scores 0.480 because it can only cut with a straight boundary between centroids, while DBSCAN recovers them (1.000).
• Unequal spread: k-means 0.404, GMM 0.967. A mixture model gives each cluster its own spread, which k-means cannot.
• Outliers scattered over a wide box: the GMM scores 0.000. The printed cluster contents show why: it puts all 600 real points in one component, which also takes 2 of the 30 added points, and gives a second, very wide component to the other 28, so the real groups are never separated. A model that must explain every point can spend a component on the outliers. The same happened on all 20 redraws (GMM separated the real groups in 0 of 20). k-means scores 1.000 on the shipped dataset, but that is the lucky draw, not the rule: on the 20 redraws it separates the real groups in only 7. In the other 13 it puts all 600 real points in one cluster and spends the second centroid on the added points, because cutting the wide box of outliers removes more squared distance than separating two tight blobs. So scattered outliers can wreck k-means as well, and whether they do depends on the draw.
• Outliers packed in one far spot: now k-means scores 0.000 for the opposite reason: with $k = 2$ it spends one centroid on the distant clump and splits nothing among the real points. This is not a draw-dependent result: on the 20 redraws k-means put all 600 real points in one cluster every time, and the GMM separated the real groups in none. DBSCAN scores 1.000 by forming a third cluster for the clump and keeping the two real groups intact. So the outlier weakness is real and depends on where the outliers are: they hurt when they are far, numerous or clumped enough to win a centroid, and the diagnostic is to look at cluster sizes and the farthest points. Apart from the two redraw lines, each row of the table comes from one generated dataset, so it illustrates the pattern and does not rank the methods.

How to decide for a segmentation problem

1. Ask what the segments are for. If marketing needs a handful of nameable profiles, a centroid-based method is a natural fit.
2. Look at the data. Standardise, take logs of skewed money and count features, and plot a 2-D projection. Crescents, rings and very uneven spread argue against k-means.
3. Run k-means anyway as the baseline. Check silhouette, stability across seeds and subsamples, and cluster sizes.
4. Compare one alternative that fixes the specific symptom: a GMM for unequal spread or overlapping groups, DBSCAN for arbitrary shapes plus noise, hierarchical clustering (repeatedly merging the two closest groups into a tree, which you can cut at any level) when a tree of nested segments is useful.
5. Choose on usefulness. Keep the simplest method whose segments differ in behaviour you can act on.

What flips the choice: many outliers or unknown shape (density-based), the need for membership probabilities (mixture), or a requirement to scale to very many points with minimal cost (k-means or its mini-batch variant, which updates centroids from small random batches of points instead of all points each pass, trading a little accuracy for speed).

Common mistakes

• Treating the output as proof that clusters exist; k-means returns $k$ clusters on any data.
• Skipping scaling.
• Switching to a fancier method without a diagnostic showing the baseline is failing.

The idea in plain words

A clustering is statistically good when the groups are compact and well separated in the numbers you fed it. It is actionable when the groups differ in something the business can act on: a different offer, price, channel, service level or product plan that would change an outcome. The first is a property of the geometry of the data. The second is a property of the decision. Nothing forces them to agree, because an internal measure (a score such as the silhouette, which compares each point's distance to its own cluster with its distance to the nearest other cluster, or inertia, the total squared distance of points to their cluster centres) never sees the business outcome. It only sees the feature table. The algorithm used below is k-means: choose k centres, assign every customer to the nearest centre, move each centre to the average of its customers, and repeat until nothing moves. A "typical customer" in a segment is that centre. The silhouette runs from -1 to 1: near 1 means points sit much closer to their own cluster than to any other, near 0 means they sit on the boundary, and negative means they look closer to another cluster. A common rule of thumb reads roughly above 0.7 as strong structure, 0.5 to 0.7 as reasonable, 0.25 to 0.5 as weak and below 0.25 as little structure.

So a clustering can score very well and be useless whenever the strongest structure in the features is not the structure that matters for the decision. Typical causes: one feature dominates by scale or units, the structure is an artefact (currency, data source, a time-zone export), or the groups are real but identical in behaviour that you care about.

Statistical quality versus actionability

Actionability checks, in the order to run them: (1) do segments differ on a held-out outcome by more than noise; (2) can a new customer be assigned from data available at decision time; (3) is each segment large enough to treat and different enough to treat differently; (4) does the segmentation stay the same under a resample or next month's refresh; (5) is there a distinct action per segment.

A case where the score is high and the clustering is useless

The data below are synthetic: 3,000 customers. Forty percent are billed in a currency whose amounts are about 100 times larger, so spend is around 10,000 for those accounts and around 100 for the rest (the clusters should therefore split roughly 60/40 if they just follow the currency). The chance of buying after a 10% offer depends only on disc_share (the share of past orders that used a discount) and visits (visits per month), never on the billing currency. Clustering A feeds raw units to k-means with $k=2$. Clustering B uses only the two habit features, standardised (shifted to mean 0 and scaled to standard deviation 1). The silhouette values are on different feature spaces, so compare the response gaps, not the silhouettes, across A and B.

import numpy as np
from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score, silhouette_score
from sklearn.preprocessing import StandardScaler

rng = np.random.default_rng(0)
n = 3000
# Synthetic customers. Forty percent of the accounts are billed in a currency whose numbers are about 100x larger.
big_units = rng.random(n) < 0.4
spend = np.where(big_units, rng.normal(10000, 2500, n), rng.normal(100, 25, n))
spend = np.clip(spend, 1, None)
visits = rng.gamma(3.0, 2.0, n)           # visits per month
disc_share = rng.beta(2, 5, n)            # share of past orders that used a discount

# Response to a 10% offer depends on discount habit and visits, not on billing currency.
logit = -3.0 + 5.0 * disc_share + 0.10 * visits
buys = rng.random(n) < 1 / (1 + np.exp(-logit))

def report(name, labels, X):
    print(name, "k=2  silhouette=%.3f" % silhouette_score(X, labels))
    for c in (0, 1):
        m = labels == c
        print("  cluster %d: size=%4d  offer response=%.3f  mean discount share=%.2f  mean visits=%.1f"
              % (c, m.sum(), buys[m].mean(), disc_share[m].mean(), visits[m].mean()))

X_raw = np.column_stack([spend, visits, disc_share])
lab_a = KMeans(2, n_init=10, random_state=0).fit_predict(X_raw)
report("A raw units   ", lab_a, X_raw)

# Behavioural clustering: only habits, standardised.
X_beh = StandardScaler().fit_transform(np.column_stack([visits, disc_share]))
lab_b = KMeans(2, n_init=10, random_state=0).fit_predict(X_beh)
report("B habits only ", lab_b, X_beh)

# Decision value: how far apart are the two clusters on the outcome the business cares about?
for name, lab in (("A", lab_a), ("B", lab_b)):
    rates = [buys[lab == c].mean() for c in (0, 1)]
    print(name, "response-rate gap between clusters = %.3f" % abs(rates[0] - rates[1]))
print("overall response =", round(buys.mean(), 3))

# Two checks from the list below. Bootstrap = redraw n customers with replacement from the data, many times.
def gap_interval(lab, reps=500):          # 95% interval for (response of cluster 0) - (response of cluster 1)
    g = []
    for b in range(reps):
        i = np.random.default_rng(b).integers(0, n, n)
        g.append(buys[i][lab[i] == 0].mean() - buys[i][lab[i] == 1].mean())
    return np.percentile(g, [2.5, 97.5])

def stability(X, lab, reps=10):           # refit on each resample, assign all customers to its centres, compare labels
    out = []
    for b in range(reps):
        i = np.random.default_rng(b).integers(0, n, n)
        c = KMeans(2, n_init=5, random_state=b).fit(X[i]).cluster_centers_
        out.append(adjusted_rand_score(lab, np.linalg.norm(X[:, None] - c[None], axis=2).argmin(1)))
    return np.mean(out)

for name, lab, X in (("A", lab_a, X_raw), ("B", lab_b, X_beh)):
    lo, hi = gap_interval(lab)
    print(name, "gap (cluster 0 minus cluster 1) 95%% interval %.3f to %.3f, bootstrap ARI %.2f" % (lo, hi, stability(X, lab)))

A raw units    k=2  silhouette=0.875
  cluster 0: size=1809  offer response=0.311  mean discount share=0.28  mean visits=6.1
  cluster 1: size=1191  offer response=0.298  mean discount share=0.29  mean visits=6.1
B habits only  k=2  silhouette=0.346
  cluster 0: size=1225  offer response=0.464  mean discount share=0.44  mean visits=6.9
  cluster 1: size=1775  offer response=0.197  mean discount share=0.18  mean visits=5.6
A response-rate gap between clusters = 0.013
B response-rate gap between clusters = 0.267
overall response = 0.306
A gap (cluster 0 minus cluster 1) 95% interval -0.021 to 0.049, bootstrap ARI 1.00
B gap (cluster 0 minus cluster 1) 95% interval 0.234 to 0.299, bootstrap ARI 0.84

Reading it:

• Clustering A has silhouette 0.875, which the rule of thumb above calls strong structure. The two clusters are simply the two billing currencies (their sizes, 1,809 and 1,191, are close to the 60/40 split built into the data). The offer response is 0.311 against 0.298, a gap of 0.013, and discount share and visits are almost identical across the two clusters. A marketer handed this segmentation would send both groups the same message, so it has no use.
• Clustering B has a much weaker silhouette (0.346, in the weak band) because habits form a continuum, not two blobs. But its clusters respond at 0.464 and 0.197, a gap of 0.267. Sending the offer only to the first would reach customers who respond more than twice as often. That is the actionable one.

The last two lines of the output show the two checks from the list below on these clusters. A bootstrap redraws the same number of customers with replacement, many times, and recomputes a number each time; the spread of those numbers is the interval. For A the interval for the response gap runs from -0.021 to 0.049, which includes zero, so the clusters are not shown to differ on the outcome at all, even though A is perfectly stable (bootstrap ARI 1.00: refitting on resampled customers puts every customer in the same cluster again, ARI being 1 for identical groupings). For B the interval is 0.234 to 0.299, clear of zero, though B is less stable (ARI 0.84) because habits form a continuum and the boundary moves a little with the sample. Stability says the clusters reproduce; only the interval on an outcome says they matter.

The same lesson holds in real work: the silhouette tells you the grouping is clear in the features you chose, never that the features were the right ones.

Making a clustering actionable

1. Start from the decision (which lever, which outcome), then choose features that plausibly move it.
2. Put features on a common scale, or convert units (currency) before clustering. Scale decides which feature dominates the distance.
3. Profile every cluster on an outcome that was not used for clustering (response, churn, margin). Report the between-cluster gap with an interval (the bootstrap interval in the code above), not only internal scores.
4. Check stability (re-cluster on a bootstrap sample and compare assignments with an agreement score such as ARI, as in the code above) and that a new customer can be placed from data available at decision time.
5. Keep the number of segments to what the team can run differently.

If no clustering separates the outcome, say so: a supervised model or a simple rule on the outcome may serve the decision better than segments.

Common mistakes

• Choosing k by silhouette or elbow alone and presenting the result as "the customer segments".
• Treating a high score as evidence the segments mean something commercially.
• Clustering on identifiers, units or collection artefacts (currency, source system, region of export) and then interpreting them as behaviour.
• Never checking a segment's size, stability or reachability.

Clusters rediscover a tier when the tier variable alone is enough to reproduce the cluster labels. So the check is to ask three things in order: how much do the clusters agree with the tier, can the tier variable alone predict the cluster, and does any structure survive once the tier variable's influence is removed from the features? If the answer to the first two is "almost entirely" and to the third is "no", you have rebuilt the tier with extra steps. If clusters agree little with the tier and structure survives, they carry information the business does not already track.

The checks, in the order to run them

1. Agreement with the tier. Compare the two partitions (two ways of splitting the same customers into groups) with the adjusted Rand index, ARI. ARI is 1 when the partitions are identical and close to 0 when they agree no more than two random partitions would. Normalised mutual information (NMI) is a similar score: it measures how much knowing one customer's tier tells you about their cluster, scaled from 0 (nothing) to 1 (everything). Also look at the plain cross-table of tier by cluster: it shows which tier each cluster spills across.
2. Predict the cluster from the tier variable alone. Fit a very small model (here a decision tree with at most two levels of splits) on the single spend column to predict the cluster label, with cross-validation (accuracy measured on rows the tree was not fitted on), and compare it with always guessing the biggest cluster. ARI is harsh when a tier boundary falls somewhere other than where k-means put its boundary, even though both are cutting the same one-dimensional line. This check does not have that problem.
3. Remove the tier variable's influence and re-cluster. Regress each standardised feature (rescaled to mean 0 and standard deviation 1) on spend, meaning fit the straight line that best predicts that feature from spend, and keep the residuals (what is left of each customer's value after subtracting that line, the part spend does not explain), cluster those, and compare with the original clusters. If the new clusters bear no resemblance to the old ones, the old ones were spend. If they match the old ones, the clusters are driven by something spend does not explain.
4. Outcome value, which needs real data. Take a business outcome (retention, margin, response to an offer) and ask whether cluster membership predicts it better than the tier does, for example by comparing held-out AUC (the chance that the model ranks a random customer who had the outcome above one who did not, measured on rows the model was not fitted on; 0.5 is a coin flip, 1 is perfect) of a model with the tier only against a model with tier plus cluster. The last two lines of each block in the run below do exactly this. If the clusters add nothing and would not change any action, they are not worth a new name.

Also keep the spend column out of the clustering inputs if the goal is to find something other than spend. Clustering on spend and then noticing the clusters track spend proves nothing.

A computed example on two synthetic customer sets

Two datasets of 1,200 customers each. In spend_only orders and discount share are functions of yearly spend plus noise, while basket size and days between orders are pure noise. In cross_cutting three behaviours (steady, discount hunter, burst buyer) are assigned independently of spend and change orders, discount share and days between orders. The tier is the tercile of spend (the bottom, middle and top third of customers by spend). K-means is the method that places 3 centres and assigns every customer to the nearest one. The program also makes a repeat-purchase outcome that depends on the tier and, in cross_cutting, on being a burst buyer. K-means with k=3 runs on the standardised features, including spend. The data is synthetic and the seed is fixed, so the numbers show how the checks behave, not what real data looks like.

import numpy as np
import pandas as pd
from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.tree import DecisionTreeClassifier

rng = np.random.default_rng(0)

def make(kind, n=1200):
    spend = np.exp(rng.normal(5, 0.7, n))                    # yearly spend
    if kind == "spend_only":                                  # every feature is driven by spend
        orders = spend / 40 * np.exp(rng.normal(0, 0.15, n))
        disc = 0.3 - 0.0005 * spend + rng.normal(0, 0.03, n)
        gap = 20 + rng.normal(0, 3, n)
    else:                                                     # three behaviours cut across spend
        g = rng.integers(0, 3, n)                             # 0 steady, 1 discount hunter, 2 burst buyer
        orders = spend / 40 * np.exp(rng.normal(0, 0.15, n)) * np.array([1, 1.6, 0.5])[g]
        disc = np.array([0.1, 0.6, 0.2])[g] + rng.normal(0, 0.05, n)
        gap = np.array([20, 12, 60])[g] + rng.normal(0, 3, n)
    basket = spend / orders
    burst = (g == 2) if kind != "spend_only" else np.zeros(n, bool)
    return burst, pd.DataFrame(dict(spend=np.log(spend), orders=np.log(orders), basket=np.log(basket), disc=disc, gap=gap))

for kind in ["spend_only", "cross_cutting"]:
    burst, df = make(kind)
    tier = pd.qcut(df.spend, 3, labels=False).to_numpy()      # business tier: terciles of spend
    X = StandardScaler().fit_transform(df)                    # column 0 is spend
    lab = KMeans(3, n_init=10, random_state=0).fit_predict(X)
    print(kind)
    print("  ARI of clusters vs tier:", round(adjusted_rand_score(tier, lab), 3))
    acc = cross_val_score(DecisionTreeClassifier(max_depth=2, random_state=0), X[:, [0]], lab, cv=5).mean()
    print("  cluster predicted from spend alone (depth-2 tree, 5-fold accuracy):", round(acc, 3),
          " | always guess the biggest cluster:", round(np.bincount(lab).max() / len(lab), 3))
    s = X[:, [0]]                                             # remove what spend explains from every feature
    Xr = np.delete(X - s @ np.linalg.lstsq(s, X, rcond=None)[0], 0, axis=1)
    lab_r = KMeans(3, n_init=10, random_state=0).fit_predict(Xr)
    print("  clusters of the spend-free features: ARI vs tier", round(adjusted_rand_score(tier, lab_r), 3),
          ", ARI vs original clusters", round(adjusted_rand_score(lab, lab_r), 3))
    # Check 4: does the cluster predict an outcome (repeat purchase) better than the tier alone? Held-out AUC, 5-fold.
    rng_y = np.random.default_rng(1)
    logit = 0.9 * (tier - 1) - 1.5 * burst              # the outcome depends on the tier and, in cross_cutting, on burst buyers
    y = (rng_y.random(len(df)) < 1 / (1 + np.exp(-logit))).astype(int)
    tier_oh = np.eye(3)[tier]
    clus_oh = np.eye(3)[lab]
    for name, F in [("tier only", tier_oh), ("tier + cluster", np.c_[tier_oh, clus_oh])]:
        auc = cross_val_score(LogisticRegression(max_iter=1000), F, y, cv=5, scoring="roc_auc").mean()
        print(f"  held-out AUC, {name:15s}: {auc:.3f}")

Output:

spend_only
  ARI of clusters vs tier: 0.459
  cluster predicted from spend alone (depth-2 tree, 5-fold accuracy): 0.875  | always guess the biggest cluster: 0.472
  clusters of the spend-free features: ARI vs tier -0.002 , ARI vs original clusters 0.014
  held-out AUC, tier only      : 0.680
  held-out AUC, tier + cluster : 0.686
cross_cutting
  ARI of clusters vs tier: -0.0
  cluster predicted from spend alone (depth-2 tree, 5-fold accuracy): 0.309  | always guess the biggest cluster: 0.336
  clusters of the spend-free features: ARI vs tier -0.0 , ARI vs original clusters 1.0
  held-out AUC, tier only      : 0.660
  held-out AUC, tier + cluster : 0.737

Reading the run:

• spend_only: ARI against the tier is 0.459, moderate, because the tercile cut points and the k-means cut points are in different places. Check 2 shows the real picture: spend alone predicts the cluster label with accuracy 0.875 against 0.472 for guessing the biggest cluster. After spend is removed, the residual clusters have ARI 0.014 with the original clusters, so nothing else was holding them together. These clusters are a re-cut spend scale.
• cross_cutting: ARI against the tier prints as -0.0 (about zero, a rounded tiny negative value), spend predicts the cluster with accuracy 0.309 against 0.336 for guessing the biggest cluster, and re-clustering the spend-free features reproduces the original clusters (ARI 1.0 on this run). Here the clusters are not the tier.
• Check 4, the outcome test, in the last two lines of each block. In spend_only the tier alone gives held-out AUC 0.680 and adding the cluster gives 0.686: the cluster is the tier re-cut, so it adds almost nothing a decision could use. In cross_cutting the tier alone gives 0.660 and tier plus cluster gives 0.737, because the clusters pick up the burst buyers, who behave differently whatever their spend. On this run the 0.077 gain (its size varies with the sample) is what "carries information the business does not already track" looks like in numbers.
• A low ARI alone is not proof of novelty and a moderate ARI alone is not proof of redundancy. The first three checks together decide; the cross-table and a profile of each cluster in business terms make the answer readable for stakeholders.

Walking through the residual line

The line Xr = np.delete(X - s @ np.linalg.lstsq(s, X, rcond=None)[0], 0, axis=1) is check 3. s is the standardised spend column (1,200 rows by 1 column). np.linalg.lstsq(s, X)[0] fits, for every column of X, the one slope that best predicts that column from spend, so it returns one slope per column. s @ slopes is then the part of every feature that spend alone predicts, and X - ... subtracts it, leaving the residuals. Column 0 is spend itself, which becomes exactly zero after subtracting its own fit, so np.delete(..., 0, axis=1) drops it before clustering.

A three-customer picture of one feature. Standardised spend is $-1, 0, 1$ and standardised orders is $-1.2, 0.3, 0.9$. The best-fit slope through the origin is $(-1)(-1.2)+(0)(0.3)+(1)(0.9)=2.1$ divided by $(-1)^2+0^2+1^2=2$, which is 1.05. Spend predicts orders of $-1.05, 0, 1.05$, and the residuals are $-1.2-(-1.05)=-0.15$, $0.3-0=0.3$ and $0.9-1.05=-0.15$. The first and third customers look ordinary once spend is accounted for; the middle one orders more than their spend suggests. Clustering these leftovers asks whether customers still group together after spend is out of the picture.

What these checks cannot see

• They detect redundancy with the one variable you pass in. Clusters could still mostly restate a different variable (account age, plan type, region), so repeat with each existing business segmentation.
• Residualising removes only the linear dependence on spend. If the relationship is curved, transform spend first (the example uses logs) or use a flexible model for the residuals.
• A cluster can be new but useless. Check 4 is what answers whether it matters; stability of the clusters across resamples (re-run on bootstrap samples and compare labels with ARI) tells you whether it is real.

Common mistakes

• Reporting a single ARI with no cross-table, so a moderate value gets read as either "different" or "same".
• Including the tier variable in the features and then treating the correlation as a finding.
• Treating "clusters differ in mean spend" as a discovery: any partition of customers on any correlated feature will.
• Skipping the outcome check and judging only by shape.

Direct answer

The elbow method fits k-means for a range of cluster counts $k$, plots the inertia (the sum of squared distances from each point to its cluster's centroid) against $k$, and picks the $k$ where the curve stops dropping steeply and flattens, like a bent arm. Inertia always falls as $k$ grows (at $k = n$ it is zero), so the bend marks the point where an extra cluster stops buying much. When there is no clear elbow, that is information: either the data has no strong cluster structure, or the structure is not the compact round kind inertia rewards, and you switch to other evidence rather than forcing a pick.

How to build the plot

1. Standardise the features, since inertia uses Euclidean distance and an unscaled column dominates it.
2. For each $k = 1, \dots, K$ (for example $K = 10$), fit k-means with several restarts (n_init) and a fixed seed, and record inertia_. Restarts matter because one unlucky start gives a noisy point and a false kink.
3. Plot $k$ on the horizontal axis and inertia on the vertical axis, joining the points.
4. Read the shape. A useful number to print next to the plot is the percentage drop from $k$ to $k+1$, since the eye is poor at judging slope.
5. To automate the choice, one heuristic finds the point farthest from the straight line joining the first and last points after rescaling both axes to $[0, 1]$; published methods build on similar ideas (for example Kneedle, which looks for the point where a rising or falling curve bends most sharply away from a straight line), and all are heuristics. They only need a curve that bends one way, a convex shape like a ski slope that starts steep and flattens.

Worked example: a clear elbow and a curve with none

The program builds two datasets: four Gaussian blobs, and one Gaussian blob with no real groups. For each it prints the inertia for $k = 1$ to $10$, the percentage drop to the next $k$, and the result of the chord heuristic above, also prints the rescaled coordinates behind that heuristic, and a second rule built on the drop column (the first $k$ whose next drop is under 40% of the drop before it), and saves the plot (elbow.png) with both curves scaled by their value at $k = 1$. A final pair of lines tallies what both rules pick on 20 other draws of each dataset.

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs

ks = np.arange(1, 11)


def inertia_curve(X):
    return np.array([KMeans(int(k), n_init=10, random_state=0).fit(X).inertia_ for k in ks])


def chord_scores(ks, w):
    """Rescale both axes to [0, 1]. The first point becomes (0, 1) and the last (1, 0), so the
    straight line joining them is x + y = 1, and |x + y - 1| is proportional to the distance from it."""
    x = (ks - ks[0]) / (ks[-1] - ks[0])
    y = (w - w[-1]) / (w[0] - w[-1])
    return x, y, np.abs(x + y - 1)


def chord_elbow(ks, w):
    return int(ks[np.argmax(chord_scores(ks, w)[2])])


def cliff_elbow(w, frac=0.4):
    """First k whose next drop collapses to under frac of the drop before it; None if no step does."""
    d = -np.diff(w) / w[:-1] * 100            # d[i] is the % drop from k = i + 1 to k = i + 2
    hits = [int(ks[i]) for i in range(1, len(d)) if d[i] < frac * d[i - 1]]
    return hits[0] if hits else None


X_clear, _ = make_blobs(n_samples=400, centers=4, cluster_std=0.8, random_state=2)
X_none = np.random.default_rng(0).normal(size=(400, 2))     # one blob, no real groups

fig, ax = plt.subplots(figsize=(5, 3))
for name, X in [("four blobs", X_clear), ("single gaussian blob", X_none)]:
    w = inertia_curve(X)
    drop = -np.diff(w) / w[:-1] * 100                         # % fall from k to k+1
    print(f"--- {name}")
    print("k            :", " ".join(f"{k:6d}" for k in ks))
    print("inertia      :", " ".join(f"{v:6.0f}" for v in w))
    print("% drop to k+1:", "      " + " ".join(f"{d:6.1f}" for d in drop))
    x, y, s = chord_scores(ks, w)
    print("rescaled x   :", " ".join(f"{v:6.2f}" for v in x))
    print("rescaled y   :", " ".join(f"{v:6.2f}" for v in y))
    print("|x + y - 1|  :", " ".join(f"{v:6.2f}" for v in s))
    print("chord elbow k:", chord_elbow(ks, w))
    print("cliff rule k :", cliff_elbow(w))
    ax.plot(ks, w / w[0], marker="o", label=name)
ax.set_xlabel("k"); ax.set_ylabel("inertia / inertia at k=1"); ax.legend()
fig.tight_layout(); fig.savefig("elbow.png", dpi=80)
print("plot saved:", "elbow.png")

# Do the two rules behave the same on 20 other draws of each dataset?
from collections import Counter
chord_blobs, cliff_blobs, chord_one, cliff_one = Counter(), Counter(), Counter(), Counter()
for s in range(1, 21):
    wb = inertia_curve(make_blobs(n_samples=400, centers=4, cluster_std=0.8, random_state=100 + s)[0])
    wo = inertia_curve(np.random.default_rng(100 + s).normal(size=(400, 2)))
    chord_blobs[chord_elbow(ks, wb)] += 1; cliff_blobs[cliff_elbow(wb)] += 1
    chord_one[chord_elbow(ks, wo)] += 1; cliff_one[cliff_elbow(wo)] += 1
print("20 redraws, four blobs : chord picks", dict(sorted(chord_blobs.items())), "| cliff rule picks", dict(cliff_blobs))
print("20 redraws, one blob   : chord picks", dict(sorted(chord_one.items())), "| cliff rule picks", dict(cliff_one))

Output:

--- four blobs
k            :      1      2      3      4      5      6      7      8      9     10
inertia      :  10284   4439   1144    520    473    427    383    345    309    284
% drop to k+1:         56.8   74.2   54.5    9.1    9.6   10.4    9.9   10.3    8.1
rescaled x   :   0.00   0.11   0.22   0.33   0.44   0.56   0.67   0.78   0.89   1.00
rescaled y   :   1.00   0.42   0.09   0.02   0.02   0.01   0.01   0.01   0.00   0.00
|x + y - 1|  :   0.00   0.47   0.69   0.64   0.54   0.43   0.32   0.22   0.11   0.00
chord elbow k: 3
cliff rule k : 4
--- single gaussian blob
k            :      1      2      3      4      5      6      7      8      9     10
inertia      :    799    533    348    280    236    196    174    152    137    124
% drop to k+1:         33.2   34.7   19.5   15.7   17.3   10.9   12.7    9.8    9.4
rescaled x   :   0.00   0.11   0.22   0.33   0.44   0.56   0.67   0.78   0.89   1.00
rescaled y   :   1.00   0.61   0.33   0.23   0.17   0.11   0.07   0.04   0.02   0.00
|x + y - 1|  :   0.00   0.28   0.45   0.44   0.39   0.34   0.26   0.18   0.09   0.00
chord elbow k: 3
cliff rule k : None
plot saved: elbow.png
20 redraws, four blobs : chord picks {2: 3, 3: 16, 4: 1} | cliff rule picks {4: 17, 3: 3}
20 redraws, one blob   : chord picks {3: 7, 4: 13} | cliff rule picks {None: 20}

Four blobs. The drops are 56.8%, 74.2% and 54.5% going from $k = 1$ to 2, 2 to 3 and 3 to 4, and then only 8.1 to 10.4% per step from $k = 4$ onward (9.1, 9.6, 10.4, 9.9, 10.3, 8.1). A practical reading rule compares each drop with the one before it. The cliff steps are 56.8%, 74.2% and 54.5%, and then the drop collapses to 9.1%, under a fifth of the 54.5% before it (9.1 / 54.5 is about 0.17). The elbow is the $k$ where that collapse begins, which is $k = 4$ and matches the four blobs that generated the data. The program applies this as the first $k$ whose next drop is under 40% of the drop before it, and prints 4.

The automatic chord heuristic printed 3, one below the correct 4. The trace shows why. After rescaling, the first point is $(0, 1)$ and the last is $(1, 0)$, so the chord is the line $x + y = 1$, and the distance of a point from that line is proportional to $|x + y - 1|$. For $k = 3$ the rescaled point is $(0.22, 0.09)$, giving $|0.22 + 0.09 - 1| = 0.69$; for $k = 4$ it is $(0.33, 0.02)$, giving $0.64$; for $k = 2$ it is $(0.11, 0.42)$, giving $0.47$. The score peaks at the largest single bend, and because the heuristic measures how far a point bends away from the chord rather than where the cliff ends, it lands on $k = 3$. This is not a one-off: on 20 other draws of the four-blob data the chord heuristic picked 3 sixteen times, 2 three times and 4 once, so it usually lands below the planted count, while the drop-ratio rule picked 4 seventeen times and 3 three times. Neither rule is a verdict; read the collapse in the drop column and use the automatic value as a prompt to look.

One blob. The drops are 33.2%, 34.7%, 19.5%, 15.7%, 17.3%, 10.9% and so on, a smooth decline with no step from "large" to "small". The chord heuristic still printed 3, since any decreasing curve that bends one way has a point farthest from its chord; its top score is 0.45 at $k = 3$, but the score only ranks the points of one curve. The drop-ratio rule printed None: no drop falls below 40% of the one before it (the sharpest is 19.5% after 34.7%, which is 56%). On 20 other draws of a single blob the chord heuristic picked 3 seven times and 4 thirteen times, and the drop-ratio rule found no collapse in all 20. That is the failure mode to remember: the heuristic always returns a number, including when there is nothing there. In an interview, the trustworthy evidence is the drop column (is there a cliff followed by a plateau?) plus a second measure, and the automatic pick is only a convenience.

What to do when there is no clear elbow

1. Use a second, differently-built measure: the silhouette score across the same $k$ range, the gap statistic, or the BIC of a Gaussian mixture. Agreement of two or three is stronger than one bend.
2. Test stability: if the best $k$ or the partitions change a lot across subsamples, the structure is weak.
3. Question the method: inertia rewards round, equally spread clusters. Elongated, ring-shaped or very unequal groups can hide the elbow; a density-based or mixture method may suit the data.
4. Let the use case choose: if the business can run at most five campaigns, pick the best $k$ at or below five and report how much inertia that leaves.
5. Report honestly: "no natural number of clusters; $k = 4$ chosen for usability" is a valid result.

Pitfalls

• Using one run per $k$ without restarts.
• Reading the plot with the vertical axis starting above zero, which exaggerates bends; scaling by the $k = 1$ value, as the program does, makes curves comparable.
• Treating the elbow as proof that $k$ natural clusters exist.