Unsupervised Learning: PCA and Clustering

This chapter opens Part IV, unsupervised learning. In Parts II and III every table had an outcome column to predict; here there is none, so there is nothing to score a prediction against and no train-versus-test accuracy to report. We look for structure in the inputs instead, starting with PCA and K-means. Chapters 24 to 26 then survey more methods.

1. Unsupervised does not mean “diagnosis without a lab”

Unsupervised learning looks for patterns in input values without a target label. PCA rotates measurements into directions explaining variation; K-means divides observations into groups by distance. Neither technique discovers that a cluster has diabetes merely because its members look alike on a chart.

Example: six invented patients with only pulse and SBP, and no diagnosis column at all. kmeans() still splits them into two groups, because it only looks at how close the points are.

six <- data.frame(pulse = c(62, 65, 60, 98, 102, 95),
                  sbp   = c(112, 118, 110, 150, 158, 146))
set.seed(1)
kmeans(six, centers = 2, nstart = 10)$cluster

[1] 1 1 1 2 2 2

The first three patients land in one group and the last three in the other, using nothing but the numbers. What the groups mean medically is a separate question that the algorithm cannot answer.

2. Why reduce the number of variables?

Real tables can have dozens or thousands of columns. Having many variables (many dimensions) causes three plain problems, often called the curse of dimensionality:

1. Data gets sparse. Cut each variable into 10 ranges and two variables make 100 boxes, five make 100,000, ten make ten billion. A few hundred patients cannot fill them, so most boxes are empty.
2. Distances look alike. With many variables, the nearest neighbor is almost as far away as the farthest point, so methods built on distance (K-means, kNN) lose their signal.
3. Models overfit and need more data. Each extra variable gives a model one more way to fit noise, so you need more patients to keep the fit honest.

The second point can be seen in a few lines. We draw 50 random points in 2, 10, 100 and 1000 dimensions and compare the closest pair with the most distant pair.

set.seed(19)
ratio <- sapply(c(2, 10, 100, 1000), function(p) {
  pts <- matrix(runif(50 * p), nrow = 50)   # 50 points, p variables
  d <- as.vector(dist(pts))                 # all pairwise distances
  min(d) / max(d)
})
names(ratio) <- c("2 vars", "10 vars", "100 vars", "1000 vars")
round(ratio, 2)

   2 vars   10 vars  100 vars 1000 vars 
     0.01      0.28      0.67      0.89 

The ratio is the closest distance divided by the farthest. It is near 0 (0.01) for 2 variables (some points are very close, some far) and climbs to 0.28, 0.67 and 0.89 for 10, 100 and 1000 variables, which means every pair is about equally far apart. Reducing the number of variables, by choosing some or by combining them, is one way out.

3. Derived variables: the simplest reduction

Before any algorithm, you can combine two columns into one that a clinician already understands. The shock index (pulse divided by systolic blood pressure) is one number that replaces two. BMI is another: weight in kilograms divided by height in metres squared replaces two measurements with one. The ehr table already contains bmi, so we only compute the shock index.

ehr <- readr::read_csv("data/synthetic_ehr.csv", show_col_types = FALSE)
ehr$shock_index <- ehr$pulse / ehr$sbp
summary(ehr$shock_index)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     NAs 
 0.4028  0.5370  0.6151  0.6200  0.6993  0.9891      50 

cor(ehr$shock_index, ehr$sbp, use = "complete.obs")

[1] -0.5470704

Most shock index values sit between about 0.5 and 0.7 here, with 50 missing where pulse or SBP was missing. The correlation with sbp is about -0.55: a moderate negative link, because SBP is the denominator. A derived variable can be partly redundant with its inputs like this, so check it with cor() before assuming it adds new information. It keeps your own judgment in the result, whereas PCA, which comes next, finds combinations from the data alone.

4. Prepare a numerical feature table

For this descriptive illustration, we work on the entire synthetic feature table. Do not do this when PCA/clusters will become predictive model features evaluated on a held-out test; in that case, learn the preprocessing and PCA components within training folds.

library(tidymodels)
features <- ehr |> select(age, bmi, sbp, pulse)
rec <- recipe(~ age + bmi + sbp + pulse, data = features) |>
  step_impute_median(all_numeric_predictors()) |>
  step_normalize(all_numeric_predictors())
processed <- bake(prep(rec, training = features), new_data = features)
x <- as.matrix(processed)
summary(x)

      age                 bmi                 sbp          
 Min.   :-1.706026   Min.   :-1.953756   Min.   :-2.75814  
 1st Qu.:-0.865619   1st Qu.:-0.807234   1st Qu.:-0.66251  
 Median : 0.001051   Median : 0.005491   Median :-0.01619  
 Mean   : 0.000000   Mean   : 0.000000   Mean   : 0.00000  
 3rd Qu.: 0.867720   3rd Qu.: 0.818215   3rd Qu.: 0.68888  
 Max.   : 1.708127   Max.   : 1.964738   Max.   : 2.49074  
     pulse           
 Min.   :-1.5847750  
 1st Qu.:-0.9318365  
 Median : 0.0009328  
 Mean   : 0.0000000  
 3rd Qu.: 0.9337021  
 Max.   : 1.5866406  

Why normalize? Age in years, pulse in beats/minute, and SBP in mmHg have different scales. Distance-based methods could otherwise be dominated by large-unit variables. Mean/SD standardization is a choice, not always the right choice for outliers.

PCA only helps when variables overlap, that is, when some are correlated. Check with cor() first; a small printed matrix is enough (Chapter 7, section 12, draws a similar heatmap; its numbers differ slightly from the ones below because the missing values are handled differently).

round(cor(x), 2)

        age   bmi  sbp pulse
age    1.00 -0.01 0.74  0.06
bmi   -0.01  1.00 0.34  0.02
sbp    0.74  0.34 1.00  0.07
pulse  0.06  0.02 0.07  1.00

Each cell is the correlation between two variables, from -1 to 1, and the diagonal is always 1. age and sbp are strongly correlated (0.74), bmi and sbp weakly (0.34 here, on the median-imputed 600-row table; 0.35 in Try it 8, which uses complete rows only), and pulse hardly correlates with anything. So PCA can combine age and sbp into one direction, but pulse will mostly stay a component of its own.

5. Principal component analysis

pca <- prcomp(x, center = FALSE, scale. = FALSE)
summary(pca)

Importance of components:
                          PC1    PC2    PC3     PC4
Standard deviation     1.3488 1.0070 0.9929 0.42529
Proportion of Variance 0.4548 0.2535 0.2465 0.04522
Cumulative Proportion  0.4548 0.7083 0.9548 1.00000

pca$rotation                  # Which variables contribute to each component?

             PC1         PC2         PC3          PC4
age   -0.6401196 -0.37165057  0.20162530  0.641459231
bmi   -0.2826287  0.87550558 -0.24828676  0.303256752
sbp   -0.7046768  0.03755465  0.07388643 -0.704670866
pulse -0.1174668 -0.30652157 -0.94458547  0.002089977

pca_scores <- as.data.frame(pca$x[, 1:2])
plot(pca_scores$PC1, pca_scores$PC2,
     xlab = "PC1 score", ylab = "PC2 score")

↔ Swipe the plot sideways to see all of it.

Principal components are weighted combinations of original variables. PC1 explains the greatest variance by construction, followed by PC2 orthogonal to PC1. Greatest variance is not the same as greatest disease-predictive power. A low-dimensional visual is not a validated risk model.

PC1 runs along the greatest spread of the points and PC2 runs at right angles to it; the figure shows two of the four variables to keep the idea visible.

Where a score comes from

Above, the data were standardized by the recipe, so we told prcomp() not to center or scale. Normally you give prcomp() the raw columns and let it do both with center = TRUE, scale. = TRUE. It then stores what it used in $center (the column means) and $scale (the column standard deviations), so it can treat new patients the same way.

complete <- na.omit(features)                  # prcomp cannot take missing values
pca_raw <- prcomp(complete, center = TRUE, scale. = TRUE)
pca_raw$center

      age       bmi       sbp     pulse 
 50.46228  30.44990 122.22050  74.99033 

pca_raw$scale

      age       bmi       sbp     pulse 
19.100535  7.003204 12.916162 11.027534 

# Score for patient 1 by hand: standardize, then weight by the loadings
z <- (unlist(complete[1, ]) - pca_raw$center) / pca_raw$scale
by_hand <- as.vector(z %*% pca_raw$rotation)
rbind(by_hand = by_hand, from_prcomp = pca_raw$x[1, ])

                  PC1      PC2         PC3       PC4
by_hand     0.3049897 1.466635 0.008056018 0.1829025
from_prcomp 0.3049897 1.466635 0.008056018 0.1829025

$center holds the four column means and $scale the four standard deviations. The two rows printed at the end are identical: a score is the standardized values multiplied by the loadings and added up, once per component.

Drawing the PC1 direction

PC1 is a line through the middle of the cloud. With two standardized variables you can draw it: the slope is the PC1 loading of the second variable divided by that of the first, and the line passes through the origin because the data are centered.

two <- as.data.frame(scale(complete[, c("age", "sbp")]))
pca_two <- prcomp(two)                         # already standardized
slope <- pca_two$rotation["sbp", "PC1"] / pca_two$rotation["age", "PC1"]
slope

[1] 1

ggplot(two, aes(age, sbp)) +
  geom_point(alpha = 0.4) +
  geom_abline(slope = slope, intercept = 0, colour = "firebrick", linewidth = 1.3) +
  coord_equal() +
  labs(x = "age (standardized)", y = "SBP (standardized)")

↔ Swipe the plot sideways to see all of it.

The slope is 1: for two standardized variables with a positive correlation, PC1 is always the 45-degree diagonal. The line follows the long axis of the oval, and the spread of points across the line is what PC2 describes.

6. K-means clustering

set.seed(500)
clusters <- kmeans(x, centers = 2, nstart = 25)
table(clusters$cluster)

  1   2 
307 293 

plot(pca_scores$PC1, pca_scores$PC2,
     col = clusters$cluster, pch = 19,
     xlab = "PC1", ylab = "PC2")

↔ Swipe the plot sideways to see all of it.

We asked for 2 groups arbitrarily, but the algorithm does not tell us there really are 2 medically meaningful groups. nstart=25 tries many starting configurations and keeps a lower-within-cluster-variation solution. Different k choices and preprocessing may produce different groupings.

K-means alternates between assigning each point to its nearest center and moving each center to the mean of its points, until no point changes cluster.

7. Examine labels only after clustering, and carefully

table(cluster = clusters$cluster, synthetic_label = ehr$diabetes)

       synthetic_label
cluster  No Yes
      1 299   8
      2 247  46

prop.table(table(clusters$cluster, ehr$diabetes), margin = 1)

            No        Yes
  1 0.97394137 0.02605863
  2 0.84300341 0.15699659

A post-hoc enrichment in an invented outcome is descriptive. You may call the groups “Cluster 1” and “Cluster 2,” not “high-risk” and “low-risk” without independent evidence, study design, and validation. Cluster IDs can reverse across runs and have no inherent ordering.

8. Common errors

• Applying PCA with missing values unhandled.
• Forgetting to scale predictors when using Euclidean distance.
• Using test-set measurements to calculate scaling/PCA parameters before model evaluation.
• Assuming the optimal k by one statistic implies a clinical truth.
• Treating a cluster/PC plot as evidence of causal features or diagnosis.

Example: the first two errors can be seen directly. prcomp() refuses a table with missing values, and unscaled variances differ so much that the largest-unit variable would dominate.

msg <- tryCatch(prcomp(features), error = function(e) conditionMessage(e))
msg

[1] "infinite or missing values in 'x'"

round(apply(features, 2, var, na.rm = TRUE), 1)

  age   bmi   sbp pulse 
362.5  50.4 168.6 121.2 

The error message says the input has missing or infinite values, so impute or drop rows first (as sections 4 and 5 did). The variances show why scaling matters: age has about seven times the variance of bmi (362.5 against 50.4) only because of its units, so without scaling it would dominate any distance or PCA.

Try it

1. Remove pulse, re-run PCA, and compare the first two components.

Show solution to question 1

Answer: Dropping pulse barely changes PC1, which is still mostly age and sbp (loadings about -0.64 and -0.71). PC2 is still mostly bmi; its sign flips, which is arbitrary in PCA.

Step by step

1. Keep three of the four standardized columns and run prcomp() again.

pca3 <- prcomp(x[, c("age", "bmi", "sbp")], center = FALSE, scale. = FALSE)
round(pca3$rotation[, 1:2], 2)

Output

      PC1   PC2
age -0.64  0.42
bmi -0.29 -0.91
sbp -0.71 -0.01

round(pca$rotation[, 1:2], 2)

Output

        PC1   PC2
age   -0.64 -0.37
bmi   -0.28  0.88
sbp   -0.70  0.04
pulse -0.12 -0.31

round(summary(pca3)$importance["Proportion of Variance", ], 2)

Output

 PC1  PC2  PC3 
0.60 0.34 0.06 

2. Compare the loadings. PC1 is nearly identical with and without pulse (age -0.64 in both, sbp -0.70 with pulse and -0.71 without). PC2 is dominated by bmi in both versions (0.88 with pulse, -0.91 without); the opposite signs are not a real change, because PCA may flip the sign of any component. pulse only had a small share of PC2 (-0.31). Without it, PC1 and PC2 carry 0.60 and 0.34 of the variance, so the first two components keep 94% of it.

2. Fit k = 3 clusters. How many observations appear in each?

Show solution to question 2

Answer: The 600 observations split into clusters of 189, 237 and 174.

Step by step

1. Run kmeans() on the same standardized matrix x with centers = 3, using a seed and nstart so the result can be repeated.

set.seed(500)
clusters3 <- kmeans(x, centers = 3, nstart = 25)
table(clusters3$cluster)

Output

  1   2   3 
189 237 174 

2. Read the table: the three counts add up to 600 (189 + 237 + 174) and no cluster is tiny. The sizes belong to this particular cut of the cloud; a different k, seed or preprocessing can give other sizes.

3. Explain why cluster numbers 1 and 2 have no clinical meaning.

Show solution to question 3

Answer: The numbers are only labels that the algorithm hands out. They carry no order and no medical meaning, and they can swap from one run to the next.

Step by step

1. Run K-means twice with different seeds and compare the two labelings.

set.seed(1); run_a <- kmeans(x, centers = 2, nstart = 1)$cluster
set.seed(2); run_b <- kmeans(x, centers = 2, nstart = 1)$cluster
table(run_a, run_b)

Output

     run_b
run_a   1   2
    1   8 285
    2 305   2

2. Run A’s cluster 1 is almost the same set of patients as run B’s cluster 2 (285 patients), and run A’s cluster 2 matches run B’s cluster 1 (305). The groups are nearly the same, but the numbers are swapped, and only 8 + 2 patients moved. The labels 1 and 2 changed while nothing about the patients did. 3. K-means has no rule such as “cluster 1 is the sicker group”. Any meaning has to come from looking at the members afterwards and from outside evidence.

4. Explain why a full-data PCA would leak information if its components were passed into a supervised model tested on that same dataset.

Show solution to question 4

Answer: A full-data PCA learns its components using the rows that will later be used as the test set, so test-set information flows into the model’s inputs and the test score looks better than it should.

Step by step

1. PCA picks directions from the variances and covariances of whatever rows you give it. If those include the test rows, the directions are partly shaped by them.
2. See the effect by fitting PCA on training rows only and on all rows.

set.seed(19)
train_rows <- sample(nrow(complete), 0.8 * nrow(complete))
pca_train <- prcomp(complete[train_rows, ], center = TRUE, scale. = TRUE)
pca_all <- prcomp(complete, center = TRUE, scale. = TRUE)
round(cbind(train_only = pca_train$rotation[, 1], all_rows = pca_all$rotation[, 1]), 3)

Output

      train_only all_rows
age       -0.644   -0.639
bmi       -0.255   -0.272
sbp       -0.712   -0.704
pulse     -0.114   -0.149

3. The two columns of loadings differ a little (for example pulse -0.114 against -0.149 and bmi -0.255 against -0.272): the all-rows PCA was shaped by the held-out rows as well. The difference is small here, but the procedure is still wrong because the test rows were used to build the model’s inputs. The correct order is: split first, fit the PCA on the training rows, then project the test rows onto those fixed loadings with predict(pca_train, newdata = ...).

5. Repeat the distance-ratio experiment with 5 and 20 variables. Where does the ratio start to approach 1?

Show solution to question 5

Answer: With 5 variables the ratio is 0.11 and with 20 it is 0.37, still far from 1. It reaches 0.69 at 100 variables and reaches 0.88 at about 1000.

Step by step

1. Repeat the experiment with the variable counts 5 and 20 added.

set.seed(19)
ratio2 <- sapply(c(2, 5, 10, 20, 100, 1000), function(p) {
  d <- as.vector(dist(matrix(runif(50 * p), nrow = 50)))
  min(d) / max(d)
})
names(ratio2) <- c("2", "5", "10", "20", "100", "1000")
round(ratio2, 2)

Output

   2    5   10   20  100 1000 
0.01 0.11 0.12 0.37 0.69 0.88 

2. Read the ratios: the ratio (closest distance divided by farthest) rises from 0.01 for 2 variables to 0.11, 0.12 and 0.37 for 5, 10 and 20, then 0.69 for 100 and 0.88 for 1000. A ratio near 1 means every pair of points is about equally far apart. The values for 2, 10, 100 and 1000 differ from section 2 (for example 0.12 here against 0.28 for 10 variables) because adding the counts 5 and 20 changes which random numbers each experiment receives. With only 50 random points per setting, the ratio varies noticeably from draw to draw, so the trend (rising with the number of variables) is what to read, not the exact values.

6. Compute the shock index and its correlation with sbp (use = "complete.obs"). Why is it so strongly tied to SBP?

Show solution to question 6

Answer: The correlation is about -0.55. The shock index is pulse / sbp, so it falls when SBP rises: SBP is in the denominator.

Step by step

1. Divide pulse by sbp and correlate the result with sbp, dropping missing pairs.

si <- ehr$pulse / ehr$sbp
cor(si, ehr$sbp, use = "complete.obs")

Output

[1] -0.5470704

2. The printed value is -0.547, a moderate negative link. A larger denominator makes the ratio smaller, and pulse hardly correlates with the other variables (section 4), so most of the ratio’s variation follows sbp.

7. Compute the PC2 score of patient 2 by hand from pca_raw$center, pca_raw$scale and pca_raw$rotation, and compare it with pca_raw$x.

Show solution to question 7

Answer: The PC2 score of row 2 is 0.99 by hand and from prcomp() (0.98988 both ways).

Step by step

1. Take row 2 of the complete-case table, subtract the stored means, and divide by the stored standard deviations.
2. Multiply by the PC2 column of the loadings and add up.

z2 <- (unlist(complete[2, ]) - pca_raw$center) / pca_raw$scale
c(by_hand = sum(z2 * pca_raw$rotation[, "PC2"]),
  from_prcomp = unname(pca_raw$x[2, "PC2"]))

Output

    by_hand from_prcomp 
    0.98988     0.98988 

3. The two numbers are the same, so a score is the standardized values times the loadings, added up.

8. Draw the PC1 line for bmi and sbp instead of age and sbp. Is the slope still 1? Why or why not?

Show solution to question 8

Answer: Yes, the slope is still 1. With two standardized variables PC1 is always the diagonal; a weaker correlation (0.35 for bmi and sbp) only makes the cloud rounder around the line.

Step by step

1. Standardize bmi and sbp, run prcomp(), and divide the PC1 loadings as before.

two_b <- as.data.frame(scale(complete[, c("bmi", "sbp")]))
pc_b <- prcomp(two_b)
slope_b <- pc_b$rotation["sbp", "PC1"] / pc_b$rotation["bmi", "PC1"]
slope_b

Output

[1] 1

cor(two_b)[1, 2]

Output

[1] 0.3476845

2. Draw the line on the standardized points.

ggplot(two_b, aes(bmi, sbp)) +
  geom_point(alpha = 0.4) +
  geom_abline(slope = slope_b, intercept = 0, colour = "firebrick", linewidth = 1.3) +
  coord_equal() +
  labs(x = "BMI (standardized)", y = "SBP (standardized)")

↔ Swipe the plot sideways to see all of it.

3. The line is the same 45-degree diagonal, because for two standardized variables with a positive correlation PC1 is always the diagonal. The correlation here is only 0.35, so the cloud is rounder around the line than for age and SBP.