Two clusterings of the same 486 customers came back with silhouette scores of 0.5918 and 0.2395. Higher is supposed to be better, so the first one wins. It is also completely useless, and understanding why that is true is most of what this part has to teach.

Every model in this series so far has had an answer key. Churn was labelled, so a score told us whether we were right. Clustering removes the key. Nothing in the maths knows what a good segment is, which means the algorithm will always return an answer and will never tell you the answer is meaningless. Judgement moves from the tuning loop into the setup, and the two decisions that matter most are made before you fit anything: what you feed it, and on what scale.

Where the churn project stands after Part 12

Last part we finished the supervised leg. A tuned random forest reached roughly 0.836 AUC against a logistic regression at roughly 0.833 on the same folds, and I argued that the linear model was the one to ship. This part uses the identical feature frame, 486 rows and 15 columns built in Part 6, but throws the label away. What we add is a segmentation the retention team can read, and a compressed representation of the same customers in two dimensions instead of fifteen.

Both techniques answer questions a classifier cannot. A churn probability tells you who to call. A segment tells you what to say when they pick up, and whether the twelve people in the retention team should be running four campaigns or one. That distinction is worth holding onto, because it explains why a clustering that adds nothing to your AUC can still be the most valuable thing you produce that quarter.

Scaling decides what a cluster means

K means groups points by squared Euclidean distance, so a column with a wide numeric range contributes more to that distance than a column with a narrow one. Our frame has tenure running 1 to 72, MonthlyCharges 18.95 to 116.05, and TotalCharges 19.15 to 8468.20, alongside twelve dummy columns that are only ever 0 or 1. TotalCharges spans a range roughly seven thousand times wider than a dummy. Feed that in raw and the algorithm is not clustering customers, it is cutting one column into slices.

# tested against: Python 3.10, numpy 2.2.6, pandas 2.3.3, scipy 1.15.3, scikit-learn 1.7.2
import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

df = pd.read_csv('Telco-Customer-Churn.csv')
df['TotalCharges'] = pd.to_numeric(df['TotalCharges'], errors='coerce')
df = df.dropna(subset=['TotalCharges'])

num = ['tenure', 'MonthlyCharges', 'TotalCharges']
cat = ['Contract', 'InternetService', 'PaymentMethod',
       'PaperlessBilling', 'TechSupport', 'OnlineSecurity']
X = pd.get_dummies(df[num + cat], columns=cat, drop_first=True).astype(float)
Xs = StandardScaler().fit_transform(X)
print('X shape', X.shape)

raw = KMeans(n_clusters=4, n_init=25, random_state=42).fit(X)
sca = KMeans(n_clusters=4, n_init=25, random_state=42).fit(Xs)
print('unscaled silhouette', round(silhouette_score(X, raw.labels_), 4))
print('scaled   silhouette', round(silhouette_score(Xs, sca.labels_), 4))
print('unscaled sizes', np.bincount(raw.labels_))
print('scaled   sizes', np.bincount(sca.labels_))
print(pd.DataFrame(raw.cluster_centers_, columns=X.columns)[num].round(1))

X shape (486, 15)
unscaled silhouette 0.5918
scaled   silhouette 0.2395
unscaled sizes [ 55 112 238  81]
scaled   sizes [164  94 127 101]
   tenure  MonthlyCharges  TotalCharges
0    66.1           102.5        6795.7
1    35.5            70.8        2155.9
2    14.3            49.6         512.8
3    53.8            82.9        4299.8

Those four unscaled centroids differ on TotalCharges and on essentially nothing else. Contract type, internet product and payment method all had a vote, and all of them were drowned out. Silhouette rewards this because cutting a single continuous variable into four contiguous bands genuinely does produce compact, well separated groups in that variable. Compactness is not relevance. Every dummy column in the frame was effectively ignored, and no retention manager has ever asked to see customers sorted by lifetime spend into quartiles.

Choosing a value for k

K means needs the number of clusters up front. Two standard aids exist. Inertia, the sum of squared distances from each point to its own centroid, falls monotonically as k rises, so you look for a bend rather than a minimum. Silhouette compares how close each point is to its own cluster against the nearest rival cluster, and ranges from minus one to one. Neither is a decision procedure. Both are inputs to one.

for k in range(2, 9):
    km = KMeans(n_clusters=k, n_init=25, random_state=42).fit(Xs)
    s = silhouette_score(Xs, km.labels_)
    print(f'k={k}  inertia {km.inertia_:.1f}  silhouette {s:.4f}')

Read that table as a maximiser and you pick k equals 2, which separates customers with internet service from customers without it. True, obvious, and useless: you already have that column. Read it as an engineer and you notice inertia falls hard from 5152 to 4175 to 3722 and then settles into a steady grind of roughly 200 per extra cluster. Structure worth having is spent by about k equals 4. Silhouette bottoms out at exactly that point, at 0.2395, because four real segments overlap at their boundaries in a way two blunt halves do not. Overlap is what a genuine population looks like.

My rule after enough of these: pick k from the business constraint first, then check the maths has not vetoed it. If a retention team can run four distinct campaigns and not eight, the candidate range is three to five, and the job of the metrics is to say whether anything in that range is a disaster. Nothing here was. Going the other way, letting silhouette nominate k and then asking the business to live with it, has produced two segments and a wasted fortnight for me more than once.

Four customer segments, profiled

A cluster label is a number until you profile it. Profiling means joining the labels back onto the original untransformed columns and describing each group in units a human recognises. Here is where the churn label comes back, not to fit anything, but to check whether the segments the features found happen to line up with the outcome we care about.

lab = KMeans(n_clusters=4, n_init=25, random_state=42).fit(Xs).labels_
prof = df.assign(cluster=lab).groupby('cluster').agg(
    n=('customerID', 'size'),
    tenure=('tenure', 'mean'),
    monthly=('MonthlyCharges', 'mean'),
    churn_rate=('Churn', lambda s: (s == 'Yes').mean()))
print(prof.round(3))
print('overall churn rate', round((df['Churn'] == 'Yes').mean(), 4))

Cluster 0 is the interesting one. Short tenure at 19.2 months, high monthly bill at 81.11, 82.9 percent on fiber, almost nobody on a two year contract, and 47.6 percent churn against a base rate of 25.3 percent. That is a recognisable person: someone who bought the expensive product recently, has no contractual reason to stay, and is deciding month by month. Cluster 2 is the mirror image, similar bills but 58 months of tenure and half on two year terms, churning at 4.7 percent. Same spend, opposite behaviour, and the difference is commitment rather than price.

None of that came from the label. K means was shown fifteen feature columns and no outcome, and it still sorted the population into groups whose churn rates differ by a factor of ten. When that happens, your features carry real signal about the outcome, which is a useful and free diagnostic. When it does not happen, and every cluster sits within a point or two of the base rate, treat that as a warning about your feature set well before you start blaming the model.

Principal component analysis on the same frame

PCA finds new axes that are linear combinations of your columns, ordered so the first captures the most variance, the second the most of what is left, and so on. Each new axis is orthogonal to the ones before it. You keep the first few and discard the rest, trading a little information for a lot fewer dimensions. Like k means it is distance based, so it needs standardised inputs for the same reason and with the same failure mode if you skip that step.

from sklearn.decomposition import PCA

p = PCA(random_state=42).fit(Xs)
ev = p.explained_variance_ratio_
print('first six', np.round(ev[:6], 4))
print('cumulative', np.round(np.cumsum(ev[:6]), 4))
print('components for 90 percent', int(np.argmax(np.cumsum(ev) >= 0.90) + 1), 'of', Xs.shape[1])

p2 = PCA(n_components=2, random_state=42)
Z = p2.fit_transform(Xs)
print('Z shape', Z.shape, 'kept', round(p2.explained_variance_ratio_.sum(), 4))
print(pd.Series(p2.components_[0], index=X.columns).sort_values(key=abs, ascending=False).head(4).round(3))

first six [0.3337 0.1857 0.098  0.0838 0.0647 0.0606]
cumulative [0.3337 0.5194 0.6175 0.7012 0.7659 0.8265]
components for 90 percent 8 of 15
Z shape (486, 2) kept 0.5194
OnlineSecurity_No internet service    0.406
InternetService_No                    0.406
TechSupport_No internet service       0.406
MonthlyCharges                       -0.405
dtype: float64

Those three columns load at 0.406 apiece because they are the same fact recorded three times. A customer with no internet service is flagged as having no internet service under OnlineSecurity, under TechSupport, and under InternetService itself. Our one hot encoding from Part 6 duplicated that indicator across every internet dependent product. PCA found the redundancy instantly, and PC1 is essentially a no internet axis with monthly spend loading in the opposite direction at minus 0.405. PC2 is a tenure and commitment axis, led by tenure at 0.504 and two year contract at 0.412.

Two components hold 51.9 percent of the variance, eight are needed for 90 percent, and nine for 95 percent. Components 14 and 15 explain nothing at all, to four decimal places, because those directions are exactly determined by the others. Fifteen columns of real rank about nine. That is not a criticism of the encoding, it is what one hot encoding of correlated products does, and knowing the number is more useful than pretending you have fifteen independent measurements.

Plotting the clusters in PC space explains the silhouette score better than any number could. Cluster 1, the no internet group, is genuinely isolated on the right of PC1. Clusters 0, 2 and 3 occupy the same region and separate mainly along PC2 and along axes the plot cannot show. That is a fair picture of a customer base. People do not fall into four tidy islands, and a segmentation that claimed they did would be describing an artefact.

Cluster labels as model features

A tempting move at this point is to feed the cluster label back into the classifier as a categorical feature. It sounds free. It is not, and the reason is the same leakage argument from Part 11. If you fit k means on the whole dataset and then cross validate, every training fold carries labels that were computed with the held out rows in the room. K means saw no churn values, so it feels safe, but it saw the held out feature vectors, and that is enough to inflate a score.

from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedKFold, cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.metrics import roc_auc_score

y = (df['Churn'] == 'Yes').astype(int).values
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
base = make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000))
print('baseline    ', round(cross_val_score(base, X, y, cv=cv, scoring='roc_auc').mean(), 4))

dummies = pd.get_dummies(pd.Series(lab, name='cluster'), prefix='c', drop_first=True).astype(float)
Xc = pd.concat([X.reset_index(drop=True), dummies], axis=1)
print('leaky labels', round(cross_val_score(base, Xc, y, cv=cv, scoring='roc_auc').mean(), 4))

aucs = []
for tr, te in cv.split(X, y):
    sc = StandardScaler().fit(X.iloc[tr])
    A, B = sc.transform(X.iloc[tr]), sc.transform(X.iloc[te])
    k = KMeans(n_clusters=4, n_init=25, random_state=42).fit(A)   # fit on train only
    A2 = np.hstack([A, np.eye(4)[k.labels_][:, 1:]])
    B2 = np.hstack([B, np.eye(4)[k.predict(B)][:, 1:]])
    m = LogisticRegression(max_iter=1000).fit(A2, y[tr])
    aucs.append(roc_auc_score(y[te], m.predict_proba(B2)[:, 1]))
print('honest      ', round(np.mean(aucs), 4))

baseline     0.8326
leaky labels 0.8367
honest       0.8323

Why does the honest version add nothing? Because a k means label is a coarse summary of columns the classifier already has, and logistic regression can read those columns directly at full resolution. Compressing tenure, contract and product into a four valued category throws information away rather than creating it. Cluster features earn their place when they encode something the model cannot see, a behavioural sequence for instance, or when the downstream model is too simple to find interactions on its own. Neither applies here.

PCA components as features tell a slightly different story. Two components, holding just over half the variance, scored 0.8406 AUC against 0.8326 for all fifteen raw columns, with fold standard deviations of 0.0259 and 0.0304 respectively. Discarding half the information made the model marginally better, which is a mild regularisation effect on a small sample rather than a discovery. On 486 rows and five folds, a difference of 0.008 sits inside the noise, and I would not claim it. What I would claim is that fifteen correlated columns were not buying anything over two clean ones, and that is worth knowing before someone proposes adding thirty more.

Two errors you will hit in the first hour

Both of these caught me on the way to the numbers above. First, sweeping k from 1 rather than 2, which is a natural thing to write and immediately fatal, because a silhouette is undefined when every point belongs to the same cluster.

Traceback (most recent call last):
  File 'p13.py', line 16, in <module>
    silhouette_score(Xs, km.labels_)
  File '.../sklearn/metrics/cluster/_unsupervised.py', line 35, in check_number_of_labels
    raise ValueError(
ValueError: Number of labels is 1. Valid values are 2 to n_samples - 1 (inclusive)

Second, and much nastier in production, calling predict on a batch of new customers whose dummy columns do not match the ones the model was fitted on. If a scoring batch happens to contain no fiber customers, get_dummies simply does not create that column, and you arrive here.

Traceback (most recent call last):
  File 'p13.py', line 26, in <module>
    km.predict(StandardScaler().fit_transform(Xn))
  File '.../sklearn/cluster/_kmeans.py', line 1085, in predict
    X = self._check_test_data(X)
ValueError: X has 13 features, but KMeans is expecting 15 features as input.

Fix is to stop using get_dummies for anything that will be scored later and switch to OneHotEncoder with handle_unknown set to ignore, inside a ColumnTransformer that remembers the training categories. Ad hoc dummy encoding is fine for a notebook and a liability the moment a second dataset touches the model. Part 21 rebuilds this properly as a package, and this is one of the concrete reasons it needs doing.

DBSCAN, t-SNE and where each one fits

DBSCAN groups points by density instead of distance to a centroid, finds the number of clusters itself, and can label points as noise. That sounds like an upgrade until you run it on one hot encoded data, where most of the space is empty and density is close to meaningless. At eps 2.0 it returned 26 clusters and marked 230 of 486 customers as noise. At eps 2.5 it collapsed to 3 clusters and 17 noise points. Nothing in between was usable. Density methods want continuous, roughly uniform feature spaces, and a frame that is 80 percent binary is not that.

t-SNE and UMAP are the other two names you will meet. Both produce striking two dimensional pictures, and both are for visual inspection only. Neither preserves global distances, so the gap between two visual blobs carries no reliable meaning, and t-SNE in particular will produce apparent clusters from pure noise if you set perplexity badly. Use them to look, never to decide, and never as features for a downstream model.

What I would use clustering for on the churn project

Ship the four segments as a campaign artefact, not as model features. Cluster 0 is 164 customers at 47.6 percent churn, recent buyers of the expensive product with no contract holding them, and the obvious offer is a discounted twelve or twenty four month term. Cluster 2 proves the mechanism works, since the same spend on a two year contract churns at 4.7 percent. That is a specific, testable retention hypothesis, and it came out of an algorithm that never saw a single churn value.

Keep the classifier from Part 12 exactly as it is. Cluster labels bought 0.8323 against a baseline of 0.8326, so they cost complexity and returned nothing measurable. Keep PCA in the toolkit for two jobs: reading loadings as a redundancy audit on any new feature frame, and producing the two dimensional plot when you need to show non technical stakeholders that the segments are real. Run the segment offer as a controlled experiment rather than a blanket rollout, using the design from A/B testing for analysts, because a 47.6 percent churn group will regress toward the mean on its own and you need a control arm to tell that apart from an effect.

Phase 2 closes here. Part 14 opens the advanced leg with imbalanced data, and our frame is already mildly imbalanced at roughly one churner in four, which is gentle compared to the fraud and defect problems where resampling decisions start to hurt. Before you move on, rerun the k sweep on your own copy with the scaler removed and watch silhouette climb while the segments turn to sand. Seeing a better score attached to a worse answer, on your own screen, is the lesson that sticks.

References

• scikit-learn user guide, Clustering, including the silhouette coefficient and the properties of k means
• sklearn.cluster.KMeans API reference, for the n_init auto default introduced in version 1.4
• scikit-learn user guide, Decomposing signals in components, for PCA and explained variance
• Telco Customer Churn dataset, the public IBM sample file this series works from

About The Author

Dr. Pranay Jha

Dr. Pranay Jha is a Cloud and AI Consultant with 18+ years of experience in hybrid cloud, virtualization, and enterprise infrastructure transformation. He specializes in VMware technologies, multi-cloud strategy, and Generative AI solutions. He holds a PhD in Computer Applications with research focused on Cloud and AI, has published multiple research papers, and has been a VMware vExpert since 2016 and a VMUG Community Leader.

See author's posts