Class 7/Lab 2: k-means on Thyroid Patients

Live in-browser R — edit and re-run any cell

Published

October 16, 2026

Introduction

In class we saw k-means on simulated data: Lloyd’s algorithm step by step, bad starts, scaling, and choosing \(K\). This lab uses real patients. The question is one a clinician could ask:

If we give the algorithm only the lab results, and no diagnoses, does it find the thyroid conditions by itself?

We know the answer for these patients, so we can check. In real unsupervised work there is no answer key. This lab shows what to look at when you cannot check.

The plan:

Look at the data. Five lab tests on very different scales.
Preprocessing changes the answer. k-means on raw, scaled, and transformed values.
Bad starts. Why nstart matters.
Choosing \(K\). The elbow, the silhouette, and the gap statistic do not agree.
Hierarchical clustering. Four linkages, four different answers.

Every code cell runs live in your browser via webR, so you do not need to install R. Packages and data load once when the page opens. Give it a few seconds on first load. Click Run Code to run a cell. Run the cells in order from the top, because later cells use objects that earlier cells make.

A note on `caret`

In earlier labs, caret::train() fit every model. train() tunes a model against an outcome, and clustering has no outcome. So in this lab, caret does the preprocessing with preProcess(), the same tool that train(..., preProcess = ...) used inside earlier labs. The clustering itself uses kmeans() and hclust() from base R, and silhouette() and clusGap() from the cluster package.

The data: thyroid function tests

The data come from 215 patients at a hospital in Belgium (Coomans et al., 1983). Each patient had five blood tests:

Column	Test	What it measures
`RT3U`	T3-resin uptake (%)	how much thyroid hormone binding capacity is free
`T4`	total serum thyroxine	the main hormone the thyroid makes
`T3`	total serum triiodothyronine	the active form of the hormone
`TSH`	basal TSH	the pituitary’s signal to the thyroid (“make more”)
`dTSH`	TSH change after a TRH injection	how strongly the pituitary responds

Doctors gave each patient a diagnosis from the full medical record: Normal (euthyroid), Hyper (overactive thyroid), or Hypo (underactive thyroid). We hide the diagnosis from the algorithm, and we use it only at the end to check.

150 of the 215 patients are normal. The two disease groups are small (35 and 30 patients).

What a clinician expects

A clinician reads these tests together. In hyperthyroidism, T4 and T3 are high, and the pituitary stops sending TSH. In hypothyroidism, T4 is low, and the pituitary sends a lot of TSH to push the thyroid harder. The median values per group show this:

If k-means finds groups like these, it found real clinical structure.

Scales and shapes

Two problems are visible:

The scales are different. The SD of RT3U is 13.2. The SD of T3 is 1.4. k-means uses Euclidean distance, so a difference of 10 in RT3U counts the same as a difference of 10 in T3. But 10 units of T3 is the whole range of the data.
TSH and dTSH are very skewed. Most patients have a TSH near 1, and a few have values above 30. Those few patients are far from everyone else, and they pull the centroids toward them.

Part 1 — Preprocessing changes the answer

k-means on the raw values

kmeans() needs the data and the number of clusters. nstart = 25 runs Lloyd’s algorithm from 25 random starts and keeps the best result. Part 2 explains why.

The cluster numbers are arbitrary labels. k-means does not know which cluster is “Hyper”. Read each column of the table: which patients ended up together?

The result is reasonable: cluster 1 is mostly hypo patients, and cluster 3 is mostly hyper patients. But look at the centroids. The largest gaps between them are in RT3U (90 to 126), the test with the largest SD. On the raw scale, RT3U decides most of the distances.

Scaling: give each test the same weight

preProcess() learns the mean and SD of each column. predict() then applies them. This is the same step that preProcess = c("center", "scale") did inside train() in Labs 3 to 6.

Now all 150 normal patients are together. But the hyper patients split into two clusters (16 and 19). Scaling fixed one problem and showed another: the skew. After scaling, the few patients with very high TSH are still far from everyone, and k-means uses one of its three clusters on them.

Fix the skew too: the Yeo-Johnson transformation

A log transform pulls in long right tails. But dTSH has negative values, so we cannot take its log. The Yeo-Johnson transformation is a cousin of the log (and of the Box-Cox transformation). It also works for zero and negative values, and preProcess() picks the strength of the transform for each column.

This is the best of the three:

Preprocessing	Normal together	Hyper together	Hypo together
Raw	137 / 150	25 / 35	23 / 30
Center + scale	150 / 150	19 / 35	22 / 30
Yeo-Johnson + center + scale	145 / 150	33 / 35	24 / 30

Same algorithm, same data, same \(K\). Three different answers. The preprocessing is part of the model. Choose it with care, and report it.

Do the clusters make clinical sense?

Profile each cluster on the original units, which a clinician can read:

One cluster has high T4 (16.4) and high T3 (3.5), with TSH near 1 and no TSH response. That is the hyperthyroid profile.
One cluster has low T4 (3.1) and very high TSH (11.9), with a large TSH response. That is the hypothyroid profile.
The large cluster (153 patients) has middle values on every test. That is the normal profile.

k-means never saw a diagnosis. It found the three textbook patterns from the numbers alone. This profile table is the most useful output of a cluster analysis. Without an answer key, it is the only way to decide whether clusters mean something.

A picture

Two tests carry most of the clinical signal: T4 (the hormone) and TSH (the signal that controls it). We plot them on log scales. Colour shows the cluster, and shape shows the true diagnosis.

The mistakes sit where the groups touch. A few hypo patients have only a small TSH rise, and they look normal. A few normal patients have a high T4, and they sit with the hyper patients. A clinician would call these borderline cases. k-means must put each patient in exactly one cluster, so it cannot say “borderline”. The Gaussian mixture models in Class 8 can.

Part 2 — Bad starts and `nstart`

Lloyd’s algorithm always lowers the within-cluster sum of squares (WSS). But it stops at the first point where no single step helps, a local minimum. Where it stops depends on where it starts.

With \(K = 3\), these data are easy: nearly every start finds the same answer. With more clusters, there are more ways to go wrong. Run k-means 100 times with \(K = 5\) and one start each, and compare the WSS:

About one run in five stops at a worse solution, and the worst is about 9% worse than the best. Each of those runs gives a different set of clusters, and nothing in the output warns you. So:

Always use many starts (nstart = 25 or more).
Always set a seed, so that you and your readers get the same clusters.

Part 3 — Choosing \(K\)

Here we know there are three diagnoses. In real work you do not know \(K\). We try three tools and see what each one says.

The elbow

WSS always falls when \(K\) goes up: more centroids always sit closer to the points. So we look for the bend, the \(K\) after which the gain gets small. Here, WSS falls by 346 from \(K=1\) to 2, and by 236 from 2 to 3. After that, each step gains 70 or less. The elbow is at \(K = 3\). But “where is the bend” is a judgement, and different people see it at different places.

The silhouette

The silhouette asks one question about each point: is this point closer to its own cluster than to the nearest other cluster?

For each point \(i\):

Calculate \(a(i)\). This is the mean distance from point \(i\) to all other points in its own cluster.
For each other cluster, calculate the mean distance from point \(i\) to all points in that cluster. Find the smallest of these mean distances. This is \(b(i)\). The cluster that gives \(b(i)\) is the nearest other cluster.
Calculate the silhouette width: \[s(i) = \frac{b(i) - a(i)}{\max\{a(i),\; b(i)\}}.\]

The value of \(s(i)\) is always between \(-1\) and \(1\):

If \(s(i)\) is near \(1\), point \(i\) is near its own cluster and far from the nearest other cluster. The point is in the correct cluster.
If \(s(i)\) is near \(0\), point \(i\) is at the same distance from both clusters. The point is on the border between two clusters.
If \(s(i)\) is less than \(0\), point \(i\) is nearer to the other cluster than to its own cluster. The point is probably in the wrong cluster.

To get one value for one clustering, calculate the mean of \(s(i)\) for all points. This is the average silhouette width. Do this for each value of \(K\). Use the \(K\) that gives the largest average silhouette width.

Two warnings:

The silhouette is not defined for \(K = 1\). It cannot tell you “there are no clusters”.
The silhouette prefers clusters that are compact and well separated. A large value does not prove that the clusters have a real meaning.

cluster::silhouette() needs the cluster labels and the distance matrix:

The silhouette picks \(K = 2\) (0.445), with \(K = 3\) close behind (0.406). What does \(K = 2\) look like?

\(K = 2\) puts all the hyper patients with the normal patients. It finds “hypo vs. everyone else”. The hypo patients, with their very high TSH, are the most separated group, so this split gives the most compact, best separated clusters. It is a valid answer to the geometric question. It is the wrong answer to the clinical one.

Silhouette widths for each patient

The average hides detail. Plot \(s(i)\) for every patient in the \(K = 3\) solution, sorted within each cluster, and colour each bar by diagnosis:

No patient has a negative width. The hyper patients have the smallest widths (mean 0.29): they sit close to the normal cluster. This is why the silhouette prefers to merge them with the normal patients. The plot also shows the borderline patients from the scatter plot: they are the short bars at the right end of each cluster.

The gap statistic

The gap statistic compares the WSS of the data with the WSS of reference data that have no clusters: points drawn uniformly over the range of each variable. It is the same idea as a permutation test: build a “nothing here” version of the data, and ask if the real data look different. \(B\) is the number of reference data sets. This cell takes a while.

The rule from the slides (the smallest \(K\) with Gap(\(K\)) \(\geq\) Gap(\(K+1\)) \(-\) SE) gives \(K = 1\): “no clusters”. The curve is almost flat. Gap(3) to Gap(8) are all within about 0.02 of each other.

Why does it fail here? The uniform reference fills a box. These data are one large group of normal patients with two small groups at its edges. Compared with a box, the normal group alone already looks “clustered”. The two small groups then add little extra. The gap statistic works best when the clusters have similar sizes.

Three tools, three answers

Tool	Chosen \(K\)
Elbow	3 (by eye)
Silhouette	2
Gap statistic (1-SE rule)	1
Medical knowledge	3

This happens often in practice. Each tool answers its own question about geometry. None of them knows what the clusters are for. Use the tools to find a short list of \(K\) values. Then profile the clusters for each \(K\) (as in Part 1) and choose with a domain expert.

Part 4 — Hierarchical clustering

Hierarchical clustering does not need \(K\) in advance. It merges the two closest groups, again and again, until one group remains. The linkage defines “closest” for groups of points:

complete: the largest distance between a point in one group and a point in the other
average: the mean of all such distances
single: the smallest such distance
Ward: the merge that adds the least WSS (the k-means criterion)

We use the same distance matrix d as before, so all four methods see the same data.

The boxes show where a cut into 3 groups falls. The trees have very different shapes:

Single linkage makes a long staircase. It adds one patient at a time to one big group. This is called chaining.
Complete and Ward make compact, balanced branches.

Now cut each tree at \(K = 3\) and compare with the diagnosis:

Linkage	Normal together	Hyper together	Hypo together
complete	147 / 150	24 / 35	24 / 30
average	150 / 150	15 / 35	24 / 30
single	150 / 150	—	—
Ward	136 / 150	32 / 35	24 / 30
k-means (Part 1)	145 / 150	33 / 35	24 / 30

Single linkage fails. One cluster holds 211 of the 215 patients. The other two clusters hold 3 patients and 1 patient. A 3-cluster cut of a chained tree only splits off the last few outliers.
Ward is the closest to k-means. This is expected, because Ward merges by the same criterion that k-means minimizes.
Average linkage keeps the normal group whole, but it loses half of the hyper patients into it.

The choice of linkage is a modelling choice, like the preprocessing in Part 1. It can change the answer as much as the choice of algorithm.

Try it yourself

1. Four clusters. Re-run k-means with centers = 4. Which group does the fourth cluster split? Profile the clusters on the original units (copy the profile code from Part 1). Is the new split clinically meaningful, or does k-means cut a continuous group in two?

2. Fewer tests. A clinic often orders only T4 and TSH. Cluster with those two columns only (use z[, c("T4", "TSH")]). How many patients does each diagnosis lose, compared with all five tests?

3. Scaled only, then the silhouette. Re-run the silhouette loop from Part 3 on z_scale (centered and scaled, no Yeo-Johnson). Which \(K\) does it pick now? Compare with the table in Part 1: the preprocessing with the better match to the diagnoses has the lower silhouette. What does this tell you about using the silhouette to choose the preprocessing?

Things worth remembering

Preprocessing is part of the model. Raw, scaled, and transformed data gave three different clusterings. Always scale. Check for skew, and transform when one tail dominates.
Use many starts and set a seed. One start can stop at a worse local minimum, and nothing warns you.
The tools for choosing \(K\) often disagree. The elbow, the silhouette, and the gap statistic each answer a geometric question. Use them to find a short list, then decide with the profiles and with someone who knows the domain.
Profile your clusters. A table of cluster medians on the original units is what turns “cluster 2” into “hypothyroid”.
Linkage matters. Single linkage chains. Ward behaves most like k-means.