Learning objectives

By the end of this session you should be able to:

• Explain why PCA looks for directions of largest variance, and when it helps to scale. the input first.
• Run PCA with prcomp() and read the parts of the result you actually use: sdev, rotation, x.
• Produce a scree plot and a biplot, in base R and in ggplot2, and say what each one shows.
• Use the common rules of thumb for how many PCs to keep, and their limits.
• Compute distances with dist(), build a tree with hclust(), and cut it into groups with cutree().
• Put the two together on a small transcriptomics-style matrix.

Why PCA?

Real datasets often have many correlated variables. Looking at them two at a time does not scale — 4 features give 6 scatter plots, 20 features give 190.

When two variables are strongly correlated, one direction already carries most of what both of them say. PCA finds such directions, ranks them by how much variance they capture, and lets you plot the first couple on one page.

Three places this comes up in bioinformatics:

• Sample QC — spot an outlier sample on PC1–PC2 before any differential analysis.
• Batch effects — check whether the largest source of variance is your biology or the day the samples were run.
• Compression — pass a few PCs into a downstream model instead of thousands of genes.

What PCA does

• PC1 — the direction of largest variance in the data.
• PC2 — the direction of largest variance at right angles to PC1.
• PC\(k\) — largest variance orthogonal to all earlier PCs.
• There are at most \(\min(n-1,\ p)\) PCs.

The transformation is a rotation of the coordinate axes. Nothing is lost until you decide to drop some components — that is the step that actually reduces dimensionality.

Step 1 — center, and usually scale

iris_num <- iris[, 1:4]
apply(iris_num, 2, var)                 # variances differ across columns

Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
   0.6856935    0.1899794    3.1162779    0.5810063 

iris_z <- scale(iris_num)               # center + divide by SD
round(apply(iris_z, 2, mean), 10)       # ~0

Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
           0            0            0            0 

round(apply(iris_z, 2, sd),   10)       # 1

Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
           1            1            1            1 

• Centering subtracts the column mean; it does not change variances.
• Scaling divides by the SD, so every column has SD 1. PCA then runs on the correlation matrix rather than the covariance matrix.

A rule of thumb: variables in the same units and comparable spread can be left unscaled; variables on different units are usually better scaled, so one large-range column does not dominate PC1 on its own.

Step 2 — prcomp()

pca <- prcomp(iris[, 1:4], center = TRUE, scale. = TRUE)
names(pca)

[1] "sdev"     "rotation" "center"   "scale"    "x"       

The parts you use most:

• pca$sdev — standard deviation of each PC (\(\sqrt{\text{eigenvalue}}\)).
• pca$rotation — the loadings: how variables combine into PCs (one row per variable).
• pca$x — the scores: where each observation lands (one row per observation).

The summary() table

pca <- prcomp(iris[, 1:4], center = TRUE, scale. = TRUE)
summary(pca)

Importance of components:
                          PC1    PC2     PC3     PC4
Standard deviation     1.7084 0.9560 0.38309 0.14393
Proportion of Variance 0.7296 0.2285 0.03669 0.00518
Cumulative Proportion  0.7296 0.9581 0.99482 1.00000

Read the Proportion of Variance row first — it is the headline. For standardised iris, PC1 ≈ 73 % and PC2 ≈ 23 %, so the first two axes together carry about 96 %.

Variance explained — scree plot

pca <- prcomp(iris[, 1:4], scale. = TRUE)
screeplot(pca, type = "lines",
          main = "")
abline(h = 1, lty = 2, col = "#D55E00")

df <- data.frame(
  PC  = factor(paste0("PC", 1:4),
               levels = paste0("PC", 1:4)),
  eig = pca$sdev^2)

ggplot(df, aes(PC, eig, group = 1)) +
  geom_col(fill = "#0072B2", alpha = 0.7) +
  geom_line() + geom_point(size = 2) +
  geom_hline(yintercept = 1, linetype = 2,
             colour = "#D55E00") +
  labs(y = "Eigenvalue")

The dashed line marks eigenvalue = 1 (the Kaiser rule). Two PCs sit above it for standardised iris.

Loadings — what each PC means

pca <- prcomp(iris[, 1:4], scale. = TRUE)
round(pca$rotation, 2)

               PC1   PC2   PC3   PC4
Sepal.Length  0.52 -0.38  0.72  0.26
Sepal.Width  -0.27 -0.92 -0.24 -0.12
Petal.Length  0.58 -0.02 -0.14 -0.80
Petal.Width   0.56 -0.07 -0.63  0.52

• Rows are the original variables, columns are the PCs.
• A large absolute value means the variable contributes strongly to that PC.
• The sign gives the direction of the contribution.

For standardised iris, PC1 is positive on all four measurements — an overall “flower size” axis. PC2 puts Sepal.Width opposite the petal measurements — a shape contrast.

The biplot

pca <- prcomp(iris[, 1:4], scale. = TRUE)
biplot(pca, cex = 0.6,
       col = c("grey60", "#D55E00"))

sc <- as.data.frame(pca$x)
sc$Species <- iris$Species
L <- as.data.frame(pca$rotation) * 2.5
L$var <- rownames(L)

ggplot(sc, aes(PC1, PC2)) +
  geom_point(aes(colour = Species), alpha = .8) +
  geom_segment(data = L,
    aes(0, 0, xend = PC1, yend = PC2),
    arrow = arrow(length = unit(.15, "cm")),
    colour = "#D55E00") +
  geom_text(data = L,
    aes(PC1 * 1.1, PC2 * 1.1, label = var),
    colour = "#D55E00", size = 3)

Points are the scores; arrows are the loadings. Arrows close together suggest correlated variables; an arrow’s length reflects how well that variable is shown in these two PCs.

How many PCs to keep?

Three common rules, none of them definitive:

1. Kaiser — keep PCs with eigenvalue > 1 (meaningful only after scaling).
2. Cumulative variance — keep enough PCs to reach, say, 80–90 %.
3. Scree elbow — the visual bend in the plot.

PCA — things to watch for

Clustering — the goal

Unsupervised: no labels, just the question are there natural groups here?

Two families to know:

• Partitioning (k-means, PAM) — you choose \(k\) first; you get \(k\) flat groups.
• Hierarchical — build a tree of nested groupings, then cut it wherever you like to get any \(k\).

Agglomerative hierarchical clustering works like this:

1. every observation starts as its own cluster,
2. repeatedly merge the two closest clusters,
3. stop when a single cluster remains,
4. cut the tree at a chosen height to read off groups.

Distance — the input to hclust()

iris_z <- scale(iris[, 1:4])
d      <- dist(iris_z, method = "euclidean")

as.matrix(d)[1:5, 1:5] |> round(2)      # first 5 x 5 block

     1    2    3    4    5
1 0.00 1.17 0.84 1.10 0.26
2 1.17 0.00 0.52 0.43 1.38
3 0.84 0.52 0.00 0.28 0.99
4 1.10 0.43 0.28 0.00 1.25
5 0.26 1.38 0.99 1.25 0.00

The same scaling point returns here: with very different units, dist() is dominated by whichever variable has the largest range, so it is usually worth scaling first.

cor(as.vector(dist(iris[, 1:4])),
    as.vector(dist(scale(iris[, 1:4]))))

[1] 0.9529472

Linkage — distance between clusters

Once points form clusters, you need a rule for the distance between two clusters. Four common choices behave differently:

d <- dist(scale(iris[, 1:4]))
op <- par(mfrow = c(1, 2), mar = c(2,3,2,1))
plot(hclust(d, method = "single"),
     labels = FALSE, main = "single",
     xlab = "", sub = "")
plot(hclust(d, method = "complete"),
     labels = FALSE, main = "complete",
     xlab = "", sub = "")

par(op)

op <- par(mfrow = c(1, 2), mar = c(2,3,2,1))
plot(hclust(d, method = "average"),
     labels = FALSE, main = "average",
     xlab = "", sub = "")
plot(hclust(d, method = "ward.D2"),
     labels = FALSE, main = "ward.D2",
     xlab = "", sub = "")

par(op)

• single — minimum pairwise distance; can produce long, chained clusters.
• complete — maximum pairwise distance; compact clusters.
• average — the mean; a middle ground.
• ward.D2 — minimises the increase in within-cluster variance; a common default for continuous data.

Cutting the tree

hclust() returns a tree; cutree() turns it into flat labels.

iris_z <- scale(iris[, 1:4])
hc <- hclust(dist(iris_z), method = "ward.D2")
cl <- cutree(hc, k = 3)

table(cluster = cl, species = iris$Species)

       species
cluster setosa versicolor virginica
      1     49          0         0
      2      1         27         2
      3      0         23        48

table(cluster, truth) is a quick check when the truth is known. For iris with Ward linkage, setosa comes out cleanly; versicolor and virginica overlap in feature space and mix a little.

Checking your clusters

Cutting a dendrogram always returns groups — whether or not real structure exists. Two quick numeric sanity checks:

library(cluster)
d  <- dist(scale(iris[, 1:4]))
hc <- hclust(d, method = "ward.D2")
cl <- cutree(hc, k = 3)

mean(silhouette(cl, d)[, "sil_width"])   # -1 to 1, higher = more separated

[1] 0.446689

cor(d, cophenetic(hc))                    # 0 to 1, tree vs. original distances

[1] 0.8226304

• Silhouette width compares each point’s distance to its own cluster against the next-nearest one. Here it comes out middling (~0.45) — matching the versicolor/virginica overlap seen in table(cluster, species).
• Cophenetic correlation checks how faithfully the tree’s merge heights reproduce the original distances (~0.82 here — a reasonably good summary).

Both are diagnostics, not proof: pair them with table(cluster, truth) when a known grouping exists, and treat any fixed cutoff as a rough guide, not a law.

Colouring the dendrogram

hc <- hclust(dist(scale(iris[, 1:4])), method = "ward.D2")
plot(hc, labels = FALSE, main = "iris — ward.D2",
     xlab = "", sub = "")
rect.hclust(hc, k = 3,
            border = c("#D55E00", "#0072B2", "#009E73"))

rect.hclust() draws the boxes for a chosen \(k\). This three-colour set (orange / blue / green) stays distinguishable under the common forms of colour-blindness.

Heatmap: matrix and dendrogram together

heatmap() clusters rows (and, by default, columns) and draws the data matrix next to the dendrogram — a fast way to show sample/feature structure at a glance, and the natural next step from a plain dendrogram.

iris_z  <- scale(iris[, 1:4])
sp_col  <- c(setosa = "#D55E00", versicolor = "#0072B2", virginica = "#009E73")

heatmap(iris_z,
        hclustfun = function(d) hclust(d, method = "ward.D2"),
        Colv = NA, labRow = NA, margins = c(6, 2),
        RowSideColors = sp_col[iris$Species])

The side colour bar (orange / blue / green, same palette as rect.hclust() above) shows the true species next to the tree heatmap() built on its own — a quick visual check of where the clustering agrees or disagrees with a known label.

PCA + clustering together

The two techniques answer complementary questions — PCA asks which axes carry the variance, clustering asks which observations group together. A natural combination: cluster the data, then show the clusters on PC1–PC2.

iris_z <- scale(iris[, 1:4])
pca <- prcomp(iris_z)
hc  <- hclust(dist(iris_z), method = "ward.D2")
cl  <- factor(cutree(hc, k = 3))

df <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2],
                 cluster = cl, species = iris$Species)

ggplot(df, aes(PC1, PC2)) +
  geom_point(aes(colour = cluster, shape = species), size = 3, alpha = 0.8) +
  stat_ellipse(aes(colour = cluster), level = 0.68) +
  labs(subtitle = "colour = cluster, shape = true species")

Summary

A few more, briefly:

• prcomp() gives you sdev, rotation, and x — most work uses those three.
• Kaiser, cumulative variance, and the scree elbow are guides, not laws.
• hclust() + cutree() + table(cluster, truth) is a fast sanity loop on any small dataset.

Session info

R version 4.4.1 (2024-06-14)
Platform: x86_64-conda-linux-gnu
Running under: Ubuntu 20.04.6 LTS

Matrix products: default
BLAS/LAPACK: /home/aabdallah/anaconda3/envs/R44/lib/libopenblasp-r0.3.30.so;  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=de_DE.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       

time zone: Europe/Berlin
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] cluster_2.1.8.1 ggplot2_4.0.0  

loaded via a namespace (and not attached):
 [1] vctrs_0.6.5        cli_3.6.5          knitr_1.50         rlang_1.1.6       
 [5] xfun_0.52          generics_0.1.4     S7_0.2.0           jsonlite_2.0.0    
 [9] labeling_0.4.3     glue_1.8.0         htmltools_0.5.8.1  scales_1.4.0      
[13] rmarkdown_2.29     grid_4.4.1         evaluate_1.0.4     tibble_3.3.0      
[17] MASS_7.3-64        fastmap_1.2.0      yaml_2.3.10        lifecycle_1.0.4   
[21] compiler_4.4.1     dplyr_1.1.4        RColorBrewer_1.1-3 pkgconfig_2.0.3   
[25] farver_2.1.2       digest_0.6.37      R6_2.6.1           tidyselect_1.2.1  
[29] pillar_1.10.2      magrittr_2.0.3     withr_3.0.2        tools_4.4.1       
[33] gtable_0.3.6