←Back to the bookChapter 195 · Clustering: Choosing k Honestly

Chapter 195 · Capstone 33 · Clustering

Choosing k Honestly

Silhouette said two segments. Stability said two segments. The answer was four, and more than half the customers were not in any of them. The thing that found it was a variable the algorithm never saw.

import numpy as np, pandas as pd, warnings
warnings.filterwarnings("ignore")
import matplotlib.pyplot as plt
from scipy.stats import chi2_contingency
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans, DBSCAN, AgglomerativeClustering
from sklearn.metrics import silhouette_score, adjusted_rand_score
plt.rcParams.update({"figure.dpi":110,"font.size":11,"axes.spines.top":False,"axes.spines.right":False,
    "axes.grid":True,"grid.alpha":0.22,"axes.titleweight":"bold","axes.titlesize":12.5,"axes.titlelocation":"left"})
VI, DK, LT, MUT, GD, RD = "#6d28d9", "#4c1d95", "#a78bfa", "#94a3b8", "#047857", "#dc2626"

BASE = "https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
fn = "capstone-clustering-choosing-k.xlsx"
def load(sheet):
    try:    return pd.read_excel("../../data/" + fn, sheet_name=sheet)
    except FileNotFoundError: return pd.read_excel(BASE + fn, sheet_name=sheet)
raw   = load("Data")
notes = load("Notes")
print("rows in the export:", f"{len(raw):,}")
raw.head()

rows in the export: 6,190

	customer_id	annual_spend	visits_per_year	avg_basket	tenure_months	discount_share	returns_rate	online_share	categories_bought	campaign_response
0	404953	1039.66	44	25.88	68	0.879	0.062	0.811	5	1
1	404234	476.66	12	39.27	80	0.472	0.196	0.170	3	0
2	401041	1350.16	42	33.56	2	0.790	0.133	0.910	6	1
3	402793	149.40	4	45.53	-1	0.324	0.061	0.796	6	0
4	400741	490.85	5	84.25	12	0.182	0.151	0.215	6	0

Step 1 · The brief

for line in notes.Notes.fillna(""):
    print(line)

CAPSTONE 33 - CUSTOMER SEGMENTATION, AND HOW MANY SEGMENTS THERE REALLY ARE

Twelve months of loyalty-account behavior from a general retailer, 6,000 accounts.

THE REQUEST. Marketing wants customer segments for next year's plan. The brief says
'four or five segments, with names'. Nobody has asked whether there are any.

THE JOB. Decide how many groups the behavior actually supports, decide whether they
are stable enough to build a year of spending on, and be willing to answer 'fewer
than you asked for, and most customers are not in any of them'.

IMPORTANT. campaign_response describes a campaign run the QUARTER AFTER this window.
It must NOT be used to build the segments. It is held back so the segments can be
checked against something the algorithm never saw.

A NOTE ON UNITS. annual_spend runs into the thousands, discount_share runs from 0 to 1.
Any distance-based method will be dominated by whichever column happens to have the
largest numbers unless something is done about it. What is done about it is a decision,
not a formality, and it changes the answer.

KNOWN FAULTS IN THIS EXPORT, left in deliberately:
  - a block of 190 accounts was exported twice
  - 84 accounts have a negative annual_spend: returns exceeded purchases
  - 233 accounts have tenure_months = -1 because the join date was never migrated

THERE IS AN ANSWER KEY in the Truth sheet, because this dataset was generated. Do not
open it until you have committed to an answer. Read Truth_ReadMe first.

Two things in there set up everything that follows.

The brief asks for four or five segments. It does not ask whether there are any. A clustering algorithm will never push back on that, because every clustering algorithm returns exactly as many groups as it is asked for.

campaign_response is held back. It describes a campaign run the quarter after this window, so it cannot be used to build the segments. It can be used to judge them, and by the end of this notebook it is the only thing that does the job.

Step 2 · Cleaning

d = raw.drop_duplicates(subset="customer_id").copy()
print(f"raw rows            {len(raw):,}")
print(f"after duplicates    {len(d):,}")
d = d[d.annual_spend > 0]
print(f"positive spend      {len(d):,}   (net-refund accounts removed)")
d = d[d.tenure_months > 0].reset_index(drop=True)
print(f"tenure known        {len(d):,}   (-1 means the join date was never migrated)")

F = ["annual_spend","visits_per_year","avg_basket","tenure_months",
     "discount_share","returns_rate","online_share","categories_bought"]
X = d[F].values.astype(float)
y = d.campaign_response.values
print(f"\nclustering on {len(F)} behavioral features, {len(d):,} accounts")
print(f"campaign response rate, whole base: {y.mean():.3f}")

raw rows            6,190
after duplicates    6,000
positive spend      5,919   (net-refund accounts removed)
tenure known        5,693   (-1 means the join date was never migrated)

clustering on 8 behavioral features, 5,693 accounts
campaign response rate, whole base: 0.185

Step 3 · First look

k-means measures distance, so the first question is not how many clusters there are. It is what a unit of distance means in each of these eight columns.

rng = d[F].max() - d[F].min()
C = d[F].corr()
pairs = sorted([(abs(C.loc[a, b]), a, b) for i, a in enumerate(F) for b in F[i+1:]], reverse=True)

print(f"{len(F)} behavioral features, {len(d):,} accounts")
print(f"  the widest column spans {rng.max()/rng.min():,.0f} times the range of the narrowest "
      f"({rng.idxmax()} against {rng.idxmin()})")
print(f"  most skewed: {d[F].skew().idxmax()} at {d[F].skew().max():+.2f}")
print("  strongest pairs: " + "; ".join(f"{a} and {b} {C.loc[a,b]:+.2f}" for _, a, b in pairs[:3]))
print("\nThe five-number summary of every column is the left panel; the whole "
      "correlation matrix is the right one.")

8 behavioral features, 5,693 accounts
  the widest column spans 24,673 times the range of the narrowest (annual_spend against returns_rate)
  most skewed: annual_spend at +2.15
  strongest pairs: annual_spend and visits_per_year +0.72; annual_spend and avg_basket +0.65; visits_per_year and categories_bought +0.56

The five-number summary of every column is the left panel; the whole correlation matrix is the right one.

What this shows. Annual spend runs into the thousands of dollars and returns rate is a proportion of one, so on raw units a fifty dollar difference in spend outweighs the entire range of the returns column. The widest column here spans about twenty-five thousand times the range of the narrowest. Any distance-based method run on these columns as they stand is really a method run on spend alone.

The correlation list is a second warning of a different kind. Spend, visits and basket size already move together, above 0.6 in places, so eight columns is not eight independent things to segment on.

fig, axes = plt.subplots(1, 2, figsize=(12.6, 4.6))

ax = axes[0]
data = [d[f].values for f in F]
bp = ax.boxplot(data, vert=False, patch_artist=True, widths=0.6, showfliers=False,
                tick_labels=[f.replace("_", " ") for f in F])
for patch in bp["boxes"]:
    patch.set_facecolor(LT); patch.set_alpha(0.65)
for med in bp["medians"]:
    med.set_color(DK); med.set_linewidth(2)
ax.set_xscale("log")
ax.set_xlabel("value, log scale (the only way all eight fit on one axis)")
ax.set_title("Eight columns, and no shared unit between them")
ax.tick_params(axis="y", labelsize=9)
ax.grid(axis="y", alpha=0)

ax = axes[1]
C = d[F].corr()
im = ax.imshow(C.values, cmap="RdBu_r", vmin=-1, vmax=1)
lbl = [f.replace("_", " ") for f in F]
ax.set_xticks(range(len(F))); ax.set_xticklabels(lbl, fontsize=8, rotation=90)
ax.set_yticks(range(len(F))); ax.set_yticklabels(lbl, fontsize=8)
for i in range(len(F)):
    for j in range(len(F)):
        if i != j and abs(C.values[i, j]) > 0.35:
            ax.text(j, i, f"{C.values[i,j]:.2f}", ha="center", va="center",
                    fontsize=7.5, fontweight="bold",
                    color="white" if abs(C.values[i, j]) > 0.65 else DK)
ax.set_title("Some of them are already saying the same thing")
ax.grid(False)
fig.colorbar(im, ax=ax, fraction=0.046, pad=0.03).ax.tick_params(labelsize=8)
plt.tight_layout(); plt.show()

Left: a log axis is the only way to get all eight columns onto one picture, and that is itself the finding. On the linear scale k-means actually works with, the top row dwarfs the bottom four. Right: spend, basket size and visits carry much of the same information. Standardizing fixes the units problem in the left panel and does nothing at all about the right one, and the next step shows what the first failure costs.

Step 4 · The units decide the answer

Before any method question, a scaling question. Look at what these columns actually are.

print(f"  {'feature':22s} {'mean':>12s} {'std dev':>13s}")
for f, m, s in zip(F, X.mean(0), X.std(0)):
    print(f"  {f:22s} {m:12,.3f} {s:13,.3f}")
print(f"\nannual_spend varies {X[:,0].std()/X[:,4].std():,.0f} times more than discount_share,")
print("purely because one is measured in dollars and the other is a proportion.")

  feature                        mean       std dev
  annual_spend              1,508.009     2,134.688
  visits_per_year              19.631        15.067
  avg_basket                   72.338        51.737
  tenure_months                57.098        33.350
  discount_share                0.396         0.248
  returns_rate                  0.099         0.083
  online_share                  0.583         0.237
  categories_bought             5.873         2.904

annual_spend varies 8,609 times more than discount_share,
purely because one is measured in dollars and the other is a proportion.

km_raw = KMeans(4, n_init=10, random_state=1).fit(X)
lab_raw = km_raw.labels_
print("K-MEANS ON RAW UNITS, k = 4")
print(f"  {'cluster':>8s} {'n':>6s}" + "".join(f"{f[:11]:>13s}" for f in F))
for c in range(4):
    m = lab_raw == c
    print(f"  {c:8d} {m.sum():6,}" + "".join(f"{v:13.2f}" for v in X[m].mean(0)))

def between_share(X, lab):
    out = []
    for j in range(X.shape[1]):
        g = pd.Series(X[:, j]).groupby(lab)
        bet = ((g.mean() - X[:, j].mean())**2 * g.size()).sum()
        out.append(bet/((X[:, j] - X[:, j].mean())**2).sum())
    return np.array(out)

print("\nSHARE OF EACH FEATURE'S VARIANCE THESE CLUSTERS ACCOUNT FOR")
for f, v in zip(F, between_share(X, lab_raw)):
    bar = "#"*int(round(v*40))
    print(f"  {f:22s} {v:6.3f}  {bar}")

K-MEANS ON RAW UNITS, k = 4
   cluster      n  annual_spen  visits_per_   avg_basket  tenure_mont  discount_sh  returns_rat  online_shar  categories_
         0  3,675       442.44        12.09        57.22        50.45         0.42         0.12         0.61         5.05
         1    520      5016.70        37.29       138.04        85.83         0.15         0.06         0.44         9.04
         2    347      7875.40        48.00       164.36        86.78         0.14         0.06         0.44         9.40
         3  1,151      1405.44        27.19        63.20        56.41         0.50         0.08         0.59         6.00

SHARE OF EACH FEATURE'S VARIANCE THESE CLUSTERS ACCOUNT FOR
  annual_spend            0.950  ######################################
  visits_per_year         0.554  ######################
  avg_basket              0.402  ################
  tenure_months           0.142  ######
  discount_share          0.200  ########
  returns_rate            0.072  ###
  online_share            0.068  ###
  categories_bought       0.251  ##########

The four clusters have mean annual spends of roughly 440, 1,400, 5,000 and 7,900 dollars, and they account for 95 percent of the variation in spend and 7 percent of the variation in how much a customer shops online.

This is not a segmentation. It is a ranking of customers by how much they spend, with four names on it, and it would have been quicker to produce with a sort. Euclidean distance adds up squared differences, so a column measured in thousands drowns out one measured in tenths.

Standardizing each column to mean zero and unit variance is the usual fix, and it is worth being clear that it is a decision and not a formality: it asserts that one standard deviation of discount rate is as important as one standard deviation of spend. That is an assumption about the business, and a different one, such as weighting by margin, would give different segments.

Z = StandardScaler().fit_transform(X)
print("standardized: every column now has mean 0 and standard deviation 1")
print(np.round(Z.mean(0), 6), np.round(Z.std(0), 6))

standardized: every column now has mean 0 and standard deviation 1
[-0. -0. -0.  0.  0.  0.  0.  0.] [1. 1. 1. 1. 1. 1. 1. 1.]

Step 5 · How many clusters? The internal answer

The silhouette score measures how much closer each point is to its own cluster than to the next nearest, averaged over points. It runs from -1 to 1 and higher is better.

rs = np.random.default_rng(1)
sub = rs.choice(len(Z), 3000, replace=False)      # silhouette is O(n^2), so score a subsample
sil = {}
for k in range(2, 11):
    lab = KMeans(k, n_init=10, random_state=1).fit_predict(Z)
    sil[k] = silhouette_score(Z[sub], lab[sub])
print("SILHOUETTE BY k")
for k, v in sil.items():
    print(f"  k = {k:2d}   {v:.4f}   {'#'*int(round(v*90))}")
print(f"\nbest k by silhouette: {max(sil, key=sil.get)}")

SILHOUETTE BY k
  k =  2   0.3670   #################################
  k =  3   0.3147   ############################
  k =  4   0.3242   #############################
  k =  5   0.2932   ##########################
  k =  6   0.2790   #########################
  k =  7   0.2739   #########################
  k =  8   0.2604   #######################
  k =  9   0.2207   ####################
  k = 10   0.1841   #################

best k by silhouette: 2

Two. Hold that, and run the identical procedure on data with no structure whatsoever: the same number of rows, the same number of columns, drawn independently from a uniform distribution.

U  = rs.uniform(size=Z.shape)
Un = StandardScaler().fit_transform(U)
siln = {}
for k in range(2, 11):
    lab = KMeans(k, n_init=10, random_state=1).fit_predict(Un)
    siln[k] = silhouette_score(Un[sub], lab[sub])
print("SILHOUETTE BY k, ON STRUCTURELESS DATA")
for k, v in siln.items():
    print(f"  k = {k:2d}   {v:.4f}   {'#'*int(round(v*90))}")
print(f"\nbest k on data with no clusters in it: {max(siln, key=siln.get)}")

SILHOUETTE BY k, ON STRUCTURELESS DATA
  k =  2   0.0950   #########
  k =  3   0.0814   #######
  k =  4   0.0858   ########
  k =  5   0.0862   ########
  k =  6   0.0855   ########
  k =  7   0.0877   ########
  k =  8   0.0900   ########
  k =  9   0.0914   ########
  k = 10   0.0925   ########

best k on data with no clusters in it: 2

There is nothing there, and the procedure still names a best k. It has to: k-means partitions whatever it is given, and silhouette compares partitions.

What the comparison does give you is a scale. The real data scores around 0.37 and structureless data of the same shape scores around 0.09, a factor of four. That gap is evidence there is structure. The location of the maximum is a much weaker signal than the height of it, and it is the location that everybody quotes.

This is the same argument as parallel analysis in Chapter 194: an index of fit means nothing until you know what it scores on noise.

Step 6 · Stability under resampling

A different question, and a better one: if these segments are real, they should survive being re-estimated on a different sample of customers. Refit on 80 percent subsamples and measure how closely each refit agrees with the full solution, using the adjusted Rand index.

def stability(M, k, reps=25, frac=0.80, seed=7):
    r = np.random.default_rng(seed)
    base = KMeans(k, n_init=10, random_state=1).fit(M).labels_
    out = []
    for _ in range(reps):
        idx = r.choice(len(M), int(frac*len(M)), replace=False)
        lab = KMeans(k, n_init=10, random_state=int(r.integers(1e6))).fit_predict(M[idx])
        out.append(adjusted_rand_score(base[idx], lab))
    return float(np.mean(out)), float(np.std(out))

print(f"  {'k':>3s} {'real data':>24s} {'structureless data':>26s}")
stab = {}
for k in range(2, 9):
    a, sa = stability(Z, k)
    b, sb = stability(Un, k)
    stab[k] = a
    print(f"  {k:3d}    {a:.3f} (sd {sa:.3f})      {b:.3f} (sd {sb:.3f})")
print(f"\nmost stable k on the real data: {max(stab, key=stab.get)}")

    k                real data         structureless data

    2    0.998 (sd 0.002)      0.230 (sd 0.400)

    3    0.997 (sd 0.002)      0.162 (sd 0.187)

    4    0.994 (sd 0.003)      0.151 (sd 0.133)

    5    0.987 (sd 0.005)      0.198 (sd 0.086)

    6    0.944 (sd 0.078)      0.223 (sd 0.090)

    7    0.823 (sd 0.111)      0.226 (sd 0.051)

    8    0.915 (sd 0.081)      0.220 (sd 0.054)

most stable k on the real data: 2

The contrast between the two columns is enormous. On the real data a refit on a different 80 percent of customers reproduces the same partition almost exactly, at 0.99 and above for the first few k. On structureless data the same procedure agrees with itself at around 0.15 to 0.23, which is to say barely at all: the boundaries land somewhere different every time because there is nothing for them to land on.

Stability is a good test of whether structure exists. It is a weak test of how much of it there is: here it is above 0.98 for every k from 2 to 5, so it cannot separate them, and its maximum is again at 2.

So two internal criteria, both well regarded, both point at two segments. Marketing asked for four or five. Time to find out who is right.

Step 7 · The external answer

campaign_response was never used to build any of these partitions. If a partition is picking up something real about how customers behave, the clusters should differ in how they responded to a campaign run the following quarter. If it is slicing a continuum, they will differ a little and smoothly.

print(f"overall response rate {y.mean():.3f}")
print(f"  {'k':>3s} {'lowest cluster':>16s} {'highest cluster':>17s} {'spread':>9s} {'chi-square p':>14s}")
spread = {}
for k in range(2, 9):
    lab = KMeans(k, n_init=10, random_state=1).fit_predict(Z)
    r = pd.Series(y).groupby(lab).mean()
    spread[k] = r.max() - r.min()
    p = chi2_contingency(pd.crosstab(lab, y))[1]
    print(f"  {k:3d} {r.min():16.3f} {r.max():17.3f} {spread[k]:9.3f} {p:14.3g}")

overall response rate 0.185
    k   lowest cluster   highest cluster    spread   chi-square p
    2            0.055             0.209     0.154       3.63e-27
    3            0.054             0.235     0.181       2.67e-46
    4            0.053             0.476     0.423      7.66e-179
    5            0.054             0.488     0.434      3.95e-184

    6            0.054             0.489     0.435      1.51e-181
    7            0.054             0.491     0.438      3.73e-180

    8            0.054             0.495     0.441      6.07e-180

There is the elbow, and it is not subtle. Going from three clusters to four more than doubles the spread in response, from 0.18 to 0.42. Going from four to five adds 0.011, and every further split adds a fraction of a point.

Four segments. Not because a fit index peaked there, but because the fourth split separates customers who respond to a price-off campaign at 48 percent from customers who respond at 5 percent, and the fifth split separates nothing.

That is the honest procedure for choosing k, and it needs something the clustering did not see. If no such variable exists, the number of clusters is not identified by the data, and the right answer to "how many segments are there" is a question about what the segments are going to be used for.

Step 8 · What k-means cannot say

Every point gets a cluster. There is no way for k-means to report that a customer sits in the middle of the cloud and belongs to nothing in particular, which is a problem if that describes most of the customer base.

DBSCAN can. It grows clusters out of dense regions and labels everything else noise.

print(f"  {'eps':>5s} {'clusters found':>16s} {'labeled noise':>16s}")
for eps in [0.7, 0.9, 1.1, 1.3, 1.5]:
    lab = DBSCAN(eps=eps, min_samples=25).fit_predict(Z)
    nc = len(set(lab)) - (1 if -1 in lab else 0)
    print(f"  {eps:5.1f} {nc:16d} {(lab == -1).mean():15.1%}")

    eps   clusters found    labeled noise
    0.7                2           91.9%
    0.9                4           55.4%

    1.1                2           15.4%
    1.3                2            3.6%

    1.5                1            1.0%

db = DBSCAN(eps=0.9, min_samples=25).fit_predict(Z)
print(f"DBSCAN at eps = 0.9: {len(set(db)) - (1 if -1 in db else 0)} clusters, "
      f"{(db == -1).mean():.1%} of accounts labeled noise")
print()
print(f"  {'cluster':>10s} {'n':>7s} {'response':>10s}" + "".join(f"{f[:10]:>12s}" for f in F[1:]))
for c in sorted(set(db)):
    m = db == c
    nm = "noise" if c == -1 else str(c)
    print(f"  {nm:>10s} {m.sum():7,} {y[m].mean():10.3f}" + "".join(f"{v:12.2f}" for v in X[m][:,1:].mean(0)))

DBSCAN at eps = 0.9: 4 clusters, 55.4% of accounts labeled noise

     cluster       n   response  visits_per  avg_basket  tenure_mon  discount_s  returns_ra  online_sha  categories
       noise   3,155      0.144       20.78       76.26       64.86        0.34        0.09        0.49        6.20
           0     877      0.483       33.53       26.78       52.10        0.79        0.07        0.79        7.46
           1     661      0.100        3.51      130.01       19.08        0.30        0.24        0.85        2.73
           2     819      0.126        9.27       42.45       56.66        0.33        0.06        0.52        4.74
           3     181      0.039       38.08      149.31       86.88        0.14        0.04        0.46        9.07

More than half the customer base is not in any segment. DBSCAN finds four dense regions covering about 45 percent of accounts and declines to classify the rest.

That is an unwelcome answer and it is a useful one. A customer picked at random is more likely than not to belong to none of these groups, so a plan that assigns every customer to a segment and a treatment is planning for a base that does not exist.

Three of the four dense regions are immediately recognizable from their profiles: heavy discount users who shop online and respond at 48 percent, infrequent buyers with large baskets and a 24 percent return rate, and a small group of long-tenured high spenders who respond at 4 percent. The fourth is a dense pocket of low-activity accounts, and the crosstab below shows what it really is.

Note also how sensitive DBSCAN is to eps. At 0.7 it finds two clusters and calls 92 percent noise; at 1.5 it merges everything into one. That parameter is doing the same job k did in k-means, and it is chosen the same way: by what the resulting partition separates.

truth = load("Truth")
dd = d.merge(truth, on="customer_id", how="left")
T = dd.true_segment.values
ct = pd.crosstab(pd.Series(db).replace({-1: "noise"}), T)
print("DBSCAN CLUSTERS AGAINST THE ANSWER KEY")
print(ct.to_string())

DBSCAN CLUSTERS AGAINST THE ANSWER KEY
col_0  bargain hunters  gift occasionals  loyal high spenders  unstructured middle
row_0                                                                             
0                  855                 0                    0                   22
1                    0               653                    0                    8
2                    0                 0                    0                  819
3                    0                 0                  180                    1
noise              123                99                  682                 2251

Three of the four dense regions are one true segment each, almost cleanly. The fourth, cluster 2, is drawn entirely from the unstructured middle: it is a dense patch of a continuum rather than a group, and it is exactly the kind of thing that gets a name and a campaign budget if nobody checks. Its response rate of 0.126 is indistinguishable from the base rate, which is the clue that was available without the answer key.

Step 9 · The answer key

This dataset was generated, so the truth is known. On real data this section does not exist, which is the reason for building the argument the way we did before opening it.

print("HOW THE BASE ACTUALLY BREAKS DOWN")
print(dd.true_segment.value_counts().to_string())
print(f"\nthe 'unstructured middle' is {(T == 'unstructured middle').mean():.1%} of the clean sample.")
print("Those accounts were drawn from broad uncorrelated distributions. They are not a segment.")
print("\nRESPONSE RATE BY TRUE SEGMENT")
print(pd.Series(y).groupby(T).mean().round(3).to_string())

HOW THE BASE ACTUALLY BREAKS DOWN
true_segment
unstructured middle    3101
bargain hunters         978
loyal high spenders     862
gift occasionals        752

the 'unstructured middle' is 54.5% of the clean sample.
Those accounts were drawn from broad uncorrelated distributions. They are not a segment.

RESPONSE RATE BY TRUE SEGMENT
bargain hunters        0.500
gift occasionals       0.094
loyal high spenders    0.049
unstructured middle    0.145

print("ADJUSTED RAND AGAINST THE ANSWER KEY")
print(f"  {'k':>3s} {'k-means':>10s} {'Ward':>10s}")
for k in range(2, 8):
    a = adjusted_rand_score(T, KMeans(k, n_init=10, random_state=1).fit_predict(Z))
    b = adjusted_rand_score(T, AgglomerativeClustering(k, linkage="ward").fit_predict(Z))
    print(f"  {k:3d} {a:10.3f} {b:10.3f}")
m = db != -1
print(f"\n  DBSCAN, on the {m.sum():,} accounts it was willing to classify: "
      f"{adjusted_rand_score(T[m], db[m]):.3f}")

ADJUSTED RAND AGAINST THE ANSWER KEY
    k    k-means       Ward

    2      0.321      0.316

    3      0.568      0.586

    4      0.848      0.924

    5      0.577      0.601

    6      0.471      0.527

    7      0.432      0.428

  DBSCAN, on the 2,538 accounts it was willing to classify: 0.963

Everything the external check said turns out to be right, and everything the internal indices said turns out to be wrong.

• k = 4 recovers the structure, at an adjusted Rand of 0.85 for k-means and 0.92 for Ward. k = 2, which both internal criteria preferred, scores 0.32.
• DBSCAN scores 0.96 on the accounts it was prepared to classify, and it independently identified that about 55 percent of the base belongs to nothing, against a true figure of 55 percent.
• Ward beat k-means at every k from three to six, which is worth a moment: the segments are not equally sized or equally spread, and k-means prefers clusters that are.

The general lesson is not that silhouette and stability are useless. They established that there was structure at all, which is a real question and was worth answering. They could not count it. Counting it needed evidence from outside the clustering.

Step 10 · The pictures

fig, axes = plt.subplots(1, 2, figsize=(12.6, 4.4))

ax = axes[0]
kk = sorted(sil)
ax.plot(kk, [sil[k] for k in kk],  color=VI, lw=2.8, marker="o", ms=5, label="this customer base")
ax.plot(kk, [siln[k] for k in kk], color=MUT, lw=2.2, marker="s", ms=4, ls="--", label="uniform random data")
ax.fill_between(kk, [siln[k] for k in kk], [sil[k] for k in kk], color=LT, alpha=0.20)
ax.set_ylim(0, 0.44); ax.set_xticks(kk)
ax.set_xlabel("number of clusters, k"); ax.set_ylabel("silhouette")
ax.set_title("Both curves peak; only one has a height")
ax.legend(fontsize=9, loc="upper right")

ax = axes[1]
ks2 = sorted(spread)
ax.plot(ks2, [spread[k]*100 for k in ks2], color=VI, lw=2.8, marker="o", ms=6)
ax.axvline(4, color=RD, lw=1.8, ls=":")
ax.text(4.15, min(spread.values())*100 + 2, "four", color=RD, fontsize=10, fontweight="bold")
ax.set_xticks(ks2)
ax.set_xlabel("number of clusters, k")
ax.set_ylabel("spread in campaign response (pp)")
ax.set_title("The fourth split does the work")
plt.tight_layout(); plt.show()

Left: the height of the gap between the two curves is the evidence that structure exists, and the position of the maximum is not evidence of how much. Right: the held-out campaign result, which is the only thing here that identifies four.

Step 11 · What this does not settle

Segments are a description of one year. Nothing here says a customer stays in their segment. A retention plan built on these groups needs to know the transition rates between them, which requires a second year of data and is a different analysis.

The scaling choice is still an assumption. Standardizing says one standard deviation of returns rate matters as much as one standard deviation of spend. Weighting by contribution to margin would be defensible and would give different boundaries, and there is no test that settles which is right, only a business argument.

The external variable validates the partition for one purpose. These four segments separate response to a price-off campaign very well. They may separate response to a loyalty program or a service change badly. A segmentation is not true or false, it is useful or not, and useful is always useful for something.

And the unstructured middle is not a failure of the analysis. It is a finding. Most customer bases are mostly continuum, and the interesting question about the 55 percent is not which segment to force them into but whether there is anything worth doing for them at all.

---

What to take away

1. Clustering on raw units clusters on whatever has the largest numbers. Ninety-five percent of the variation captured was spend; seven percent was online share.
2. Every method returns clusters, including on data that has none. Silhouette named a best k for uniform random data without hesitation.
3. Read the height of an index, not the position of its maximum. 0.37 against 0.09 is evidence; the argmax at 2 was wrong.
4. Stability detects structure and does not count it. 0.99 on the real data against 0.15 on noise, and above 0.98 for every k from 2 to 5.
5. The number of clusters came from a variable the algorithm never saw. The response spread doubled at the fourth split and flattened after it.
6. k-means must classify everyone. DBSCAN found that 55 percent of the base belongs to no segment, which was the truth and was not in the brief.

←Back to the bookChapter 195 · Clustering: Choosing k Honestly