Tools of segmentation: k means, latent class analysis and going beyond RFM

The two customers RFM said were the same

Last November two orders landed at Northgate Hampers within an hour of each other. Northgate is Priya's invented online hamper shop in Cork, at northgatehampers.example, the same shop as in the last two chapters.

The first came from an office manager in Cork: twelve of the same hamper, twelve delivery addresses, a company name on the invoice. €540. The second came from a woman in Waterford who sent three hampers to her three sisters, one gift message each. €190. Both were first orders, both recent, and both in the top fifth of order values, because most orders in the shop are a single hamper at around €40 to €70.

The RFM scoring Priya had set up as a quick win put them in the same cell: recency 5, frequency 1, monetary 5. The automation attached to that cell sent both a "welcome, here is 10% off your next hamper" email in December.

The office manager did not need 10% off. She needed someone to ring her next October and ask how many she wants this time. The woman in Waterford needed a reminder in late February that Mother's Day is coming. RFM saw two numbers and one behaviour. There were two behaviours, and frankly two businesses, inside one shop. This chapter is about the tools that tell them apart.

Northgate Hampers, its customers and every number here are made up so that I can show each step. The numbers come from a simulated customer file, and the code further down prints exactly what the tables show.

Where this sits in the series

This is chapter 14 of 21 in my Marketing Analytics series, and the last of part three, the interrelationship techniques. Principal components and factor analysis for marketers boiled 18 behavioural columns for 6,200 Northgate customers down to three component scores: volume, gifting orientation and deal sensitivity. Segmentation: strategy before algorithms argued that you decide what a segment is for before you touch an algorithm. This chapter runs the algorithms on those three scores. Part four opens with modelling the value of marketing communications. Every basic idea used here is built from scratch on the Foundations page.

The strategy chapter worked on two years of orders and 12,400 customers. Here, as in chapter 12, the file holds the 6,200 customers with an order in the last 14 months, so the counts are smaller and the shares differ.

What a good segmentation has to deliver

Four things, and a segmentation that fails any one of them fails.

Separation. Segments sit far apart on the things that drive behaviour. Two segments that differ only slightly get the same marketing, so they are one segment with a line drawn through it.

Cohesion. Customers inside a segment are alike enough that one message, one offer and one landing page work for most of them.

Interpretability. The person who runs the email programme can describe every segment in a sentence without a table. If the sentence needs the words "component two", start again.

Actionability. Each segment needs at least one thing you would do differently for it: channel, offer, timing or product. Otherwise it is a statistical curiosity.

The first two are measurable, and the silhouette further down measures both. The last two are judgement, and they are where most segmentations quietly die: beautifully separated clusters nobody can use.

An example: the RFM cell 515 passes the first two tests on paper, because everyone in it has the same three scores. It fails the third and fourth. Of the 49 Northgate customers who score 515, 40 are corporate buyers and nine are private gift buyers, and no single action suits both.

The toolbox in plain words

Hierarchical clustering starts with every customer as their own cluster and merges the closest pair, over and over, until one cluster remains. You cut the resulting tree, the dendrogram, where it gives the number of groups you want. Ward's linkage merges the pair that adds least to the variance within clusters, so it produces compact groups, and it is my default.¹ The cost is that every step considers all possible merges,¹ fine for 6,000 customers and unworkable for 6 million.

k means picks $k$ centre points, assigns every customer to the nearest centre, moves each centre to the mean of its members, and repeats until nothing changes.¹ It is fast, it scales, and it is what most people mean by "clustering". It also has firm opinions: numeric variables on comparable scales, roughly round groups of similar spread,¹ and it will find $k$ clusters whether or not the data contains any.

Latent class analysis (LCA) assumes the customers are a mixture of a few unobserved classes, each with its own probability of each answer, and that within a class the variables are independent.² Membership is a probability: a customer is 92% one class and 8% another. It handles yes or no answers and categories natively, which k means cannot.

Gaussian mixture models do the same for continuous variables: ellipses instead of circles, and a probability of membership instead of a hard label.³

Decision trees are the profiling tool, not a segmentation method. A shallow tree that predicts the segment from a few plain variables gives the shortest description of each group, in rules you can type into an email tool.

An example: Northgate's three component scores are continuous and already standardised, which is k means territory. Its raw flags, such as "gift message on most orders" or "company on the invoice", are yes or no answers, which is LCA territory. This chapter runs both and compares them.

The RFM score ranks, it does not describe

RFM summarises a customer's past in three numbers: recency (how long since the last order), frequency (how many orders) and monetary value (the average spend per order).⁴ Priya's version ranks every customer into fifths on each and writes the three scores side by side:

$$\text{RFM}_i = 100\,R_i + 10\,F_i + M_i$$

Here $\text{RFM}_i$ is the three digit code for customer $i$, $R_i$ the recency score from 1 to 5 with 5 for the most recent fifth, $F_i$ the frequency score and $M_i$ the monetary score, each from 1 for the lowest fifth to 5 for the highest. A score of 555 is the best cell, a 111 has probably left. Ranking into fifths means sorting and cutting, like finding a median, see Foundations: median and mode.

Worked example: two customers, one cell

In the simulated file the most recent fifth of customers ordered within the last 46 days, so both customers, at 12 days, get $R = 5$.

Frequency is lumpy. 34% of the customers have exactly one order in the 14 months, more than a fifth, so they all share the lowest score. Both customers have one order: $F = 1$.

The top fifth by average order value starts at €69.08. €540 and €190 are both above it: $M = 5$.

So $\text{RFM} = 100 \times 5 + 10 \times 1 + 5 = 515$ for both. In code, with customers 1 and 2 as the two from the opening:

import numpy as np, pandas as pd

c = pd.read_csv("hamper_customers.csv").set_index("customer_id")

def fifth(s):                          # 1 for the lowest fifth, 5 for the highest, ties share a score
    return np.ceil(5 * s.rank(method="min", pct=True)).astype(int)

c["R"] = fifth(-c.days_since_last)     # fewer days since the last order scores higher
c["F"] = fifth(c.orders)
c["M"] = fifth(c.avg_order_value)
c["RFM"] = 100 * c.R + 10 * c.F + c.M
print(c.loc[[1, 2], ["days_since_last", "orders", "avg_order_value", "RFM"]])

It prints 515 twice. The code says "new, high value customer". It cannot say that one is a company buying for twelve offices and the other a sister buying for three sisters, because neither fact is in the three numbers.

How k means works

k means measures how far a customer is from a centre with the ordinary straight line distance, just in more dimensions:

$$d(x_i, \mu_j) = \sqrt{\sum_{v=1}^{V} (x_{iv} - \mu_{jv})^2}$$

Here $d(x_i, \mu_j)$ is the distance between customer $i$ and the centre of cluster $j$, $x_i$ is the customer's row of scores, $\mu_j$ the centre's row, $V$ the number of variables (three at Northgate), $x_{iv}$ the customer's value on variable $v$ and $\mu_{jv}$ the centre's value on it. Squared differences are the raw material of the variance, see Foundations: spread, and a row of numbers is a vector, see Foundations: matrices.

The algorithm minimises the total squared distance between every customer and the centre of their own cluster:¹

$$J = \sum_{j=1}^{k} \sum_{x_i \in C_j} \lVert x_i - \mu_j \rVert^2$$

Here $J$ is the within cluster sum of squares, WCSS for short, $k$ the number of clusters, $C_j$ the set of customers in cluster $j$, $x_i \in C_j$ reads "every customer $i$ in cluster $j$", $\mu_j$ the centre of cluster $j$, which is the mean of its members, and $\lVert x_i - \mu_j \rVert^2$ the squared distance from the formula above. It is the same sum of squared deviations as the variance, only around each cluster's own mean, see Foundations: the mean and Foundations: spread.

The algorithm only guarantees a local minimum, which depends on where the centres start, so you run it from many starting points and keep the best.¹

Worked example: six customers by hand

Six customers, two scores each, written as (gifting, volume): P1 (2, minus 1), P2 (3, minus 1), P3 (2, 0), P4 (minus 1, 2), P5 (0, 2), P6 (minus 1, 1). Two clusters, and a deliberately poor start: the centres sit on P1 and P3.

1. Assign. P1 and P2 are nearest to the first centre (distances 0 and 1, against 1 and 1.41). P3 to P6 are nearest to the second. P4, for example, is $\sqrt{3^2 + 3^2} = 4.24$ from the first and $\sqrt{3^2 + 2^2} = 3.61$ from the second.
2. Update. The first centre moves to the mean of P1 and P2, (2.5, minus 1). The second moves to the mean of P3 to P6, (0, 1.25). With these groups $J = 0.5 + 8.75 = 9.25$.
3. Assign again. P3 is now 1.12 from the first centre and 2.36 from the second, so it switches. Everyone else stays.
4. Update. The centres move to (2.33, minus 0.67) and (minus 0.67, 1.67). Each cluster's squared distances add up to 1.33, so $J = 2.67$.
5. Assign again. Nobody moves, so the algorithm stops.

A poor start cost one extra round here. On real data a poor start can end in a worse local minimum altogether, which is why the Northgate run below uses 50 starts.

Choosing k: the elbow and the silhouette

For choosing $k$ I lean on two measures. The elbow plots $J$ for each $k$ and looks for the point where adding a cluster stops buying much. The silhouette compares, for each customer, how close they sit to their own cluster with how close they sit to the nearest other one:⁵

$$s(i) = \frac{b(i) - a(i)}{\max\{a(i),\, b(i)\}}$$

Here $s(i)$ is the silhouette of customer $i$, $a(i)$ the mean distance from customer $i$ to the other members of their own cluster, $b(i)$ the mean distance to the members of the closest other cluster, and $\max\{a(i), b(i)\}$ the larger of the two, which keeps $s(i)$ between minus 1 and 1. Near 1 means well placed, near 0 means on a boundary, below 0 means probably in the wrong cluster. The mean of $s(i)$ over all customers gives one number per $k$, an ordinary average, see Foundations: the mean and Foundations: ratios and rates.

An example: take P3 in the final two clusters above. Its distances to P1 and P2 are 1 and 1.41, so $a = 1.21$. Its distances to P4, P5 and P6 are 3.61, 2.83 and 3.16, so $b = 3.20$. Then $s = (3.20 - 1.21) / 3.20 = 0.62$. P3 is clearly in the right cluster, though less snugly than P1, whose silhouette is 0.74.

Worked example: k means on Northgate's three components

The input is the three component scores from chapter 12 for all 6,200 customers. They are already standardised and uncorrelated, which removes the first k means trap, unequal scales. The run uses scikit learn,¹ with $k$ from 1 to 8 and 50 random starts each:

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

Z = pd.read_csv("hamper_components.csv").set_index("customer_id")   # volume, gifting, deal
for k in range(1, 9):
    km = KMeans(n_clusters=k, n_init=50, random_state=1).fit(Z)
    sil = silhouette_score(Z, km.labels_) if k > 1 else float("nan")
    print(k, round(km.inertia_), round(sil, 2), np.bincount(km.labels_).min())

It prints the WCSS, the mean silhouette and the size of the smallest cluster for each $k$:

k	WCSS	Mean silhouette	Smallest cluster
1	18,590	none	6,200
2	13,646	0.34	1,263
3	9,025	0.39	1,254
4	5,681	0.45	447
5	4,770	0.38	431
6	4,066	0.36	431
7	3,547	0.34	431
8	3,279	0.31	425

The horizontal axis is the number of clusters, the vertical axis the within cluster sum of squares. Each of the first three extra clusters removes between 3,300 and 5,000. The fifth removes only 911, and every one after that less. The elbow sits at four.

The silhouette and the sizes tell the same story. Two clusters split the regular self buyers from everyone else, true and not much use. Three add a deal driven group but leave the corporate buyers inside the gift buyers' cluster, stretching it. Four separate the corporate buyers, the silhouette peaks at 0.45, and the smallest cluster is 447 customers, 7%, small but real. Five cut the gift buyers in two on deal sensitivity alone, and the silhouette drops. Four it is, which matches what the strategy chapter said the business could act on.

The profile

The centres, in standard units where 0 is the average customer and 1 is one standard deviation above it:

Segment	Volume	Gifting	Deal sensitivity
A Occasion gifters	minus 0.42	0.12	minus 0.82
B Corporate senders	0.59	2.86	minus 0.13
C Treat yourself regulars	1.58	minus 0.82	minus 0.08
D Bargain occasionals	minus 0.58	minus 0.26	1.03

Centres in standard units are for the analyst. The business needs the raw columns, averaged per segment, from the order history the components were built from:

import joblib
km = KMeans(n_clusters=4, n_init=50, random_state=1).fit(Z)
c = pd.read_csv("hamper_customers.csv").set_index("customer_id")
c["cluster"] = km.labels_
g = c.groupby("cluster")
profile = pd.DataFrame({
    "customers": g.size(),
    "orders_per_year": g.orders.mean() * 12 / 14,           # the window is 14 months
    "average_order": g.spend.sum() / g.orders.sum(),
    "gift_message": g.gift_message_share.mean(),
    "company_invoice": g.company_on_invoice.mean(),
    "voucher": g.voucher_share.mean(),
    "december": g.december_share.mean(),
    "revenue_share": g.spend.sum() / c.spend.sum(),
})
print(profile.round(2))
joblib.dump(km, "kmeans_k4.pkl")                            # the scoring script reuses this

Segment	Customers	Share	Orders per year	Average order	Gift message	Company on invoice	Code used	December	Revenue share
A Occasion gifters	2,424	39%	2.4	€63	88%	1%	24%	45%	24%
B Corporate senders	447	7%	1.3	€1,042	87%	91%	12%	79%	41%
C Treat yourself regulars	1,241	20%	6.9	€50	9%	1%	31%	18%	28%
D Bargain occasionals	2,088	34%	1.3	€38	57%	0%	74%	35%	7%

Gift message, code used and December are the average share of each customer's orders, company on invoice the share of customers.

Segment A is the shop's centre of gravity: two fifths of customers, a gift message on nearly every order, two or three occasions a year, few codes. They buy on a date, not on a price. Their gifting centre is only 0.12 because the corporate senders stretch that axis, one more reason to read the raw profile and not just the centres.

Segment B is 7% of customers and 41% of the revenue: 698 orders at an average of €1,042 make €726,992 of the file's €1,770,862. RFM gave 90% of them a frequency score of 1 or 2, because they order once a year, so on a frequency ranking they look like casual buyers. They are the most valuable relationships the shop has.

Segment C almost never sends a gift. They order for themselves nearly seven times a year, the natural audience for a subscription box. Their orders are spread through the year: they are not waiting for the sale.

Segment D is the group the discount automation was built for: a third of the customers, 7% of the revenue, about one order a year and a code on three orders in four. More than half their orders still carry a gift message. The name "bargain occasionals" is a hypothesis, more on that in the pitfalls.

The horizontal axis shows the four segments, the vertical axis each one's share of the 14 months of revenue. The smallest segment by headcount carries the largest share, and the biggest discount audience carries the smallest.

Back to the two customers from the opening. The office manager in Cork lands in B. The woman in Waterford lands in A. RFM saw 515 twice. The segmentation sees a corporate account to phone in October and an occasion gifter to remind in February. Neither should have been sent 10% off.

Latent class analysis: membership as a probability

LCA fits a model and maximises its likelihood. In words: the probability of one customer's pattern of answers is the sum, over all classes, of the probability of being in that class times the probability of that pattern within the class. Within a class the variables are assumed independent, so that second probability is a plain product:²

$$L(\theta) = \prod_{i=1}^{n} \sum_{c=1}^{K} \pi_c \prod_{v=1}^{V} P(x_{iv} \mid c)$$

Here $L(\theta)$ is the likelihood of all the parameters $\theta$ together, $n$ the number of customers, $K$ the number of classes, $\pi_c$ the share of customers in class $c$, $V$ the number of observed variables, $x_{iv}$ the answer customer $i$ gave on variable $v$, and $P(x_{iv} \mid c)$ the probability that a member of class $c$ gives that answer. Multiplying for "and" and adding for "either" are the basic rules of probability, see Foundations: probability. In practice the software maximises the logarithm of $L$, which turns the long product into a sum, see Foundations: logarithms.

After fitting, Bayes' rule gives every customer a posterior probability of belonging to each class:²

$$P(c \mid x_i) = \frac{\pi_c \prod_{v=1}^{V} P(x_{iv} \mid c)}{\sum_{c'=1}^{K} \pi_{c'} \prod_{v=1}^{V} P(x_{iv} \mid c')}$$

Here $P(c \mid x_i)$ is the probability that customer $i$ belongs to class $c$ given their answers $x_i$, the numerator is the share of class $c$ times the probability of those answers within it, and the denominator adds the same quantity over every class $c'$ so that the posteriors of one customer sum to 1. It is Bayes' rule exactly as in Foundations: probability.

That posterior is the practical difference from k means. k means says "you are in cluster D". LCA says "you are 0.92 D and 0.08 A", and the 0.08 is information you can act on.

Worked example: one posterior by hand

Take two classes and two yes or no answers, with illustrative numbers. Class A, the occasion gifters, is 55% of customers: 95% of them mostly send gifts, 5% use a code on most orders. Class D is 45%: 70% mostly send gifts, 95% use a code on most orders. A customer answers yes to both.

1. Class A: $0.55 \times 0.95 \times 0.05 = 0.0261$.
2. Class D: $0.45 \times 0.70 \times 0.95 = 0.2993$.
3. The sum is $0.3254$.
4. Posteriors: $P(A) = 0.0261 / 0.3254 = 0.08$ and $P(D) = 0.2993 / 0.3254 = 0.92$.

A gift buyer who uses codes is most likely a bargain occasional, but one in 12 customers with that pattern belongs to the gifters. That is the "0.92 D, 0.08 A" from above.

Choosing the number of classes with the BIC

More classes always raise the likelihood, so LCA needs a penalty for complexity. The Bayesian information criterion, the BIC, provides one, and the model with the lowest BIC wins:²⁶

$$\text{BIC} = -2 \ln L + p \ln n$$

Here $\ln L$ is the natural logarithm of the maximised likelihood, $p$ the number of free parameters (the class shares plus every answer probability per class), and $n$ the number of customers. The first term rewards fit, the second charges for every parameter, more heavily in bigger files. The logarithm is explained in Foundations: logarithms.

An example, using the four class fit below: the log likelihood is minus 17,632.67, there are 31 free parameters and 6,200 customers. So $\text{BIC} = 2 \times 17{,}632.67 + 31 \times \ln 6{,}200 = 35{,}265.34 + 31 \times 8.732 = 35{,}536.04$.

Worked example: latent class analysis as a second opinion

To check that k means was not drawing tidy circles on noise, I refit the same 6,200 customers with LCA on categorical versions of the raw columns: mostly gift orders, company on the invoice, a code on most orders, mostly December orders, three or more orders (each yes or no), and the average order in three bands. "Delivered to another address" stays out, because it nearly duplicates the gift message and would break the independence assumption.

R users have the poLCA package.² In Python the expectation maximisation loop is short enough to write out, and it shows the two formulas above at work:

c = pd.read_csv("hamper_customers.csv").set_index("customer_id")
items = pd.DataFrame({
    "gift": c.gift_message_share >= 0.5,
    "company": c.company_on_invoice == 1,
    "voucher": c.voucher_share >= 0.5,
    "december": c.december_share >= 0.5,
    "repeat": c.orders >= 3,
    "basket": np.digitize(c.avg_order_value, [50, 250]),   # under €50, €50 to €250, over €250
}).astype(int)
H = [np.eye(items[v].max() + 1)[items[v]] for v in items]   # one column per answer level

def lca(H, K, starts=10, iters=200, seed=0):
    rng = np.random.default_rng(seed)
    best = None
    for _ in range(starts):
        pi = np.full(K, 1 / K)
        theta = [rng.dirichlet(np.ones(h.shape[1]), K) for h in H]   # P(answer | class)
        for _ in range(iters):
            logp = np.log(pi) + sum(h @ np.log(t).T for h, t in zip(H, theta))
            ll = np.logaddexp.reduce(logp, axis=1)                    # log of the sum over classes
            post = np.exp(logp - ll[:, None])                         # Bayes: posterior per class
            pi = post.mean(axis=0)
            theta = [np.clip(post.T @ h / post.sum(axis=0)[:, None], 1e-9, 1) for h in H]
        if best is None or ll.sum() > best[0]:
            best = (ll.sum(), pi, theta, post)
    return best

n, fits = len(items), {}
for K in range(2, 7):
    ll, pi, theta, post = fits[K] = lca(H, K)
    p = (K - 1) + K * sum(h.shape[1] - 1 for h in H)            # free parameters
    print(K, round(-ll), p, round(-2 * ll + p * np.log(n)), np.bincount(post.argmax(1), minlength=K).min())

Like k means, the loop only finds a local maximum, so it starts ten times and keeps the best.² It prints:

Classes	Negative log likelihood	Parameters	BIC	Smallest class
2	20,754	15	41,640	2,108
3	18,701	23	37,603	430
4	17,633	31	35,536	430
5	17,628	39	35,597	430
6	17,628	47	35,666	224

The BIC bottoms out at four classes. A fifth class buys five points of log likelihood and pays 61 points of BIC for it.

Then the comparison. The class numbers LCA hands out are arbitrary, so the classes are matched to the k means clusters on their content, not on their index:

from scipy.optimize import linear_sum_assignment
post = fits[4][3]
table = pd.crosstab(km.labels_, post.argmax(1))
rows, cols = linear_sum_assignment(-table.values)              # best one to one match
print(table, round(table.values[rows, cols].sum() / n, 3))

93.7% of customers, 5,808 of 6,200, land in the matching class. Of the 392 who move, 284 move between A and D: occasion gifters who used codes, or bargain occasionals who wrote gift messages. The posteriors agree. 93.6% of customers have a largest posterior above 0.9, 4.6% between 0.7 and 0.9, and 111 customers, 1.8%, at 0.7 or below. 84 of those 111 sit between A and D. That boundary is real. k means just hides it by drawing a hard line through it.

Which tool for which job

The flow reads from top to bottom: the purpose first, then the inputs decide the method. My rule: k means when the inputs are continuous and already cleaned into components, LCA when they are mostly yes or no answers, survey items or categories, a Gaussian mixture for soft membership on continuous data. Whatever you pick, run one of the others as a check. Agreement above 75 to 80% between two methods is the strongest evidence I know that the segments live in the data rather than in the algorithm. Northgate's 93.7% clears it comfortably.

What RFM cannot see and behavioural segments can

RFM is a good ranking of value at a point in time, the right tool for "who do I contact first" when the budget stretches to one mailing. As a segmentation it fails for three reasons that sit in its definition.

It has no needs in it. Recency, frequency and monetary value say what happened, not why. Identical scores may hide a wedding, a bereavement, a contract renewal or a Tuesday.

It has no cause of purchase. A gift buyer buys because a date is coming, a corporate buyer because a budget year is closing, a self buyer because the fridge is empty. A frequency score flattens all three calendars.

It has no response in it. What did each group do the last time the shop sent a discount, a product story, a reminder? With a segment field on the customer record, that becomes a lookup.

An example, with numbers that are this example's assumptions, not outputs of the simulation: D clicks discount emails three times as often as A, A clicks the "Mother's Day is in three weeks" email more than anyone, and B opens no marketing emails at all, while an October phone call brings a reorder from about 60% of them. None of this can be measured until the segment exists.

Behavioural segments give you the trigger, the need and the response, which is what you need to write different marketing for different people, the point of segmentation is not demographics.

Running it in practice

Data. One row per customer from 12 to 24 months of orders: categories bought, whether delivery and billing addresses differ (not the addresses themselves), gift messages, codes, the company field, and email engagement if you have it. In a Solidus shop all of this sits in the spree_orders, spree_line_items, spree_addresses and spree_adjustments tables.⁷ A few hundred customers is too few for k means to mean anything. I would start from about a thousand, and trust it from about five thousand.

Steps. Build the table with 10 to 20 behavioural variables. Standardise, or reduce to components. Run k means for $k$ from 2 to 8 with many starts. Fit LCA or a Gaussian mixture as a second opinion. Pick $k$ with the business in the room. Profile, name, write the scoring script, and load the segment field into the shop and the email tool.

Checks before you trust it. Refit k means from 20 different seeds and compare each result with the original using the adjusted Rand index, which is 1 for identical groupings and close to 0 for chance agreement.⁸ Below 0.9 the solution is unstable and I would not ship it. Then the bootstrap: resample customers with replacement, refit, score everyone, compare.⁹ I want a median adjusted Rand index above 0.7.

from sklearn.metrics import adjusted_rand_score
ref = km.labels_
seeds = [adjusted_rand_score(ref, KMeans(4, n_init=10, random_state=s).fit(Z).labels_) for s in range(20)]
rng = np.random.default_rng(0)
boot = []
for b in range(100):                                   # resample customers, refit, score everyone
    sample = Z.iloc[rng.integers(0, len(Z), len(Z))]
    boot.append(adjusted_rand_score(ref, KMeans(4, n_init=10, random_state=b).fit(sample).predict(Z)))
print(round(min(seeds), 3), round(np.median(boot), 3), round(np.percentile(boot, 5), 3))

It prints 0.999 0.991 0.978: the worst of 20 seeds agrees almost perfectly, and even the weakest 5% of bootstrap refits stay near 0.98. Two more checks need time rather than code. A holdout in time: fit on the first ten months, score the last four, and check that the profile still reads the same. Predictive validity: does membership predict something the model never saw, such as response to last spring's campaign? Segments that predict nothing outside their own inputs are a description, not a tool.

Scoring. The scoring script is short, and it must reuse the fitted model rather than refit:

km = joblib.load("kmeans_k4.pkl")                      # centroids frozen at fit time
names = {0: "A", 3: "B", 2: "C", 1: "D"}               # matched on the profile, not the index
Z = pd.read_csv("hamper_components.csv").set_index("customer_id")   # this month's scores
out = pd.DataFrame(index=Z.index)
out["segment"] = [names[c] for c in km.predict(Z)]
out["distance"] = km.transform(Z).min(axis=1)          # far from every centre: flag for review
out.to_csv("customer_segments.csv")

On the fitting month it reproduces the profile exactly: 2,424 in A, 447 in B, 1,241 in C, 2,088 in D. The 1% furthest from any centre sit more than 1.98 standard units away, and they get a look by hand. In production, the month's component scores come from the scaler and weights saved in chapter 12, never from a refit.

Rescore monthly, but do not refit monthly: a refit moves the centres and relabels everyone, which breaks every campaign report that names a segment. Refit once or twice a year, as a small project with a migration table. For a shop of this size the whole job is two to three weeks, half of it on the customer table.

Pitfalls

Unscaled inputs. Order value in euro next to a gift message flag of 0 or 1: k means will cluster on order value and nothing else. Standardise, or better, cluster on components.

Circles where the data has none. If the true structure is a long thin group, say self buyers spread along a spend gradient, k means will chop it in two. Mixture models soften this, and a silhouette per cluster usually reveals it.

Rediscovering RFM. If the clustering variables are dominated by spend and frequency, the clusters will be "high, medium, low value" and you have built RFM with more steps. Keep value variables from dominating the clustering and use them to profile afterwards. That is why segment B's revenue share is a finding, not an input.

Too many segments. Every segment needs a creative brief, a landing page variant, an email flow and an owner. Four is grand for a small team. Eight is a spreadsheet nobody opens.

Static membership. People move. A D customer who buys three times at full price is now a C or an A. Score monthly and keep a migration table. A healthy programme moves people from D to A and from A to B.

Names become facts. "Bargain occasionals" is a hypothesis about 2,088 people generated from three numbers. Some are gift buyers who found the shop through a voucher site and would happily pay full price for a hamper they trust. Test the name: send half of D a full price product story instead of a code, and measure.

LCA has two of its own. Label switching: class 1 in one run is class 3 in the next, so match classes on their profiles, never on their index.² Local independence: two near duplicate indicators, like "gift message" and "delivered to another address", break the assumption and tend to invent extra classes. Drop or merge one, as above.

How I do this for clients

The data I need is read access to the order database or an export of it, 12 to 24 months, plus email engagement if it exists. Nothing beyond what the shop already holds, and never the addresses themselves.

It starts with the free workshop, in which we settle what the segments are for and which two or three actions the business could realistically take for each. Then two weeks of real work on me: the customer table, the components, k means and LCA side by side, the stability checks, the profile, the names. You can walk away at any point, and the written plan is yours to keep.

What you get:

• A scored customer table with a segment and a confidence for every customer, loaded into your shop and your email tool.
• A one page profile per segment, written for the people who will use it.
• A decision memo with the number that matters (in the example, 7% of customers bring 41% of the revenue and are treated as one off buyers), with its error margin and risk.
• The scoring script, and a dashboard that answers one question: how many customers are in each segment this month, and where did they move from?

After that, if you want it, I run the activation: one customer list per segment in the shop, the email tool and the ad platforms, with the October call list for B, the February reminder for A, the subscription offer for C, and lookalike audiences seeded from A and B alone, because seeded from the whole file they would mostly find more of D.

Your raw data stays yours, and the terms for code and models are in the analytics and data science terms. The method sits under Data Science, the activation under Performance Marketing, and the fees are on the pricing page.

Questions to ask before you sign off a segmentation

• Which variables were used to cluster, and were spend and frequency among them?
• How was the number of segments chosen, and what did the silhouette and a second method say?
• Can each segment be described in one sentence without looking at a table?
• What will the business do differently for each segment, and who owns that?
• How stable are the segments under resampling and across time?
• What is the plan for rescoring, and how will migration between segments show up?
• Which of the segment names are hypotheses that still have to be tested?

The series is inspired by Mike Grigsby's Marketing Analytics.¹⁰ The explanations, examples and numbers here are my own.

Footnotes

1. scikit learn developers, Clustering (sections on k means, hierarchical clustering and the silhouette coefficient), scikit learn user guide, https://scikit-learn.org/stable/modules/clustering.html

2. Drew A. Linzer and Jeffrey B. Lewis, poLCA: An R Package for Polytomous Variable Latent Class Analysis, Journal of Statistical Software, vol. 42, issue 10, 2011, https://www.jstatsoft.org/article/view/v042i10

3. scikit learn developers, Gaussian mixture models, scikit learn user guide, https://scikit-learn.org/stable/modules/mixture.html

4. Peter S. Fader, Bruce G. S. Hardie and Ka Lok Lee, RFM and CLV: Using Iso value Curves for Customer Base Analysis, Journal of Marketing Research, 2005 (authors' working paper), https://www.brucehardie.com/papers/rfm_clv_2005-02-16.pdf

5. Peter J. Rousseeuw, Silhouettes: a graphical aid to the interpretation and validation of cluster analysis, Journal of Computational and Applied Mathematics, vol. 20, 1987, pp. 53 to 65, https://linkinghub.elsevier.com/retrieve/pii/0377042787901257

6. Gideon Schwarz, Estimating the Dimension of a Model, The Annals of Statistics, vol. 6, no. 2, 1978, pp. 461 to 464, https://projecteuclid.org/journals/annals-of-statistics/volume-6/issue-2/Estimating-the-Dimension-of-a-Model/10.1214/aos/1176344136.full

7. Solidus contributors, core database schema (spree_orders, spree_line_items, spree_addresses, spree_adjustments), Solidus source code on GitHub, https://github.com/solidusio/solidus/blob/main/core/db/migrate/20160101010000_solidus_one_four.rb

8. Lawrence Hubert and Phipps Arabie, Comparing partitions, Journal of Classification, vol. 2, 1985, pp. 193 to 218, https://link.springer.com/article/10.1007/BF01908075

9. Christian Hennig, Cluster wise assessment of cluster stability, Computational Statistics and Data Analysis, 2007, https://linkinghub.elsevier.com/retrieve/pii/S0167947306004622

10. Mike Grigsby, Marketing Analytics: A Practical Guide to Improving Consumer Insights Using Data Techniques, Kogan Page, 2022, ISBN 9781398608191, https://www.netgalley.co.uk/catalog/book/267990