Every exercise uses the same iris dataset (150 flowers, 3 species: setosa, versicolor, virginica; 4 measurements in cm: Sepal.Length, Sepal.Width, Petal.Length, Petal.Width), but each one asks something new: comparing variable subsets, projecting a brand-new flower, picking the number of clusters without knowing the species count, changing the distance metric, and finding exactly which flowers a clustering gets “wrong”. Try each task before you open the solution.
Not everyone works through exercises at the same pace, so each one below is flagged — but the exercises build on each other (later ones reuse objects like hc and cl from earlier ones), so everyone still works through them in order, together as a group. The flags aren’t an invitation to skip ahead; they tell you what to expect, and what to do if you’re falling behind:
🟢 Easy — everyone should get through these by writing the code themselves.
🟡 Medium — a bit harder, or introduces something new; keep typing if you’re keeping pace, but copying from the solution is fine too.
🔴 Challenge — the most demanding exercises, including today’s only new programming idea (sapply()). If you’re falling behind, copy the code from the solution instead of typing it yourself — you’ll need the objects it creates for the exercises that follow — then move on together with the group.
2.1 Session 5
2.1.1 Exercise 1: How many variables do you actually need? 🟢 Core
The slides ran PCA on all four iris measurements together. Here, compare that to PCA on only the two petal measurements (Petal.Length, Petal.Width).
Run prcomp(..., center = TRUE, scale. = TRUE) on iris[, 1:4] (all four measurements) and again on iris[, 3:4] (petals only).
For each one, work out the percentage of variance PC1 alone explains.
Why does PC1 explain so much more when you only use the two petal columns?
Solution 1.
pca_all <-prcomp(iris[, 1:4], center =TRUE, scale. =TRUE)pca_petal <-prcomp(iris[, 3:4], center =TRUE, scale. =TRUE)round(100* pca_all$sdev[1]^2/sum(pca_all$sdev^2), 1)
73% versus 98.1%. Petal.Length and Petal.Width are almost collinear (their biplot arrows point nearly the same way), so one axis already captures almost all of their combined spread. Add the two sepal measurements, which correlate with the petals less, and PC1 has to share variance with the later PCs instead.
We keep using pca_all (the four-variable version) for the rest of this session.
2.1.2 Exercise 2: Where would a new flower land? 🟡 Stretch
predict() also works on a prcomp object: it projects a new observation onto the PC axes you already learned from iris, without re-running the PCA. The slides never showed this, since they only ever plotted the 150 training flowers.
Make up one new flower: Sepal.Length = 6.1, Sepal.Width = 2.8, Petal.Length = 4.7, Petal.Width = 1.2 (all in cm, same units as iris). Project it into pca_all’s PC1–PC2 space with predict(pca_all, newdata = ...). Then compare its PC1 score to the typical PC1 range for each species (group the scores by species and summarise the min and max of PC1 with dplyr). Which species does it look most like?
Hint:predict() needs a data.frame (or tibble) with column names that match the ones you used in prcomp() exactly (Sepal.Length, Sepal.Width, Petal.Length, Petal.Width). Once you have proj <- predict(pca_all, newdata = new_flower), look at proj[, "PC1"] and compare it to each species’ PC1 range from scores %>% group_by(species) %>% summarise(min = min(PC1), max = max(PC1)). The species whose range is closest wins.
This should make the exercise a lot easier — both function calls are already written above. You just fill in new_flower’s four numbers and read off the comparison.
predict() centers and scales the new row using the same mean and SD that prcomp() learned from iris, then multiplies it by the same rotation matrix. So the new point lands on the same PC axes, without touching the original 150 flowers. Its PC1 score (~0.63) falls in the versicolor/virginica overlap zone, well away from setosa’s range (roughly -2.8 to -1.8). That fits its intermediate-to-large petal measurements.
2.1.3 Exercise 3: Picking k without knowing the species count 🔴 Challenge
The slides always cut the tree at k = 3, because we already know iris has three species. In practice you rarely know the “true” number of groups ahead of time.
Build hc <- hclust(dist(scale(iris[, 1:4])), method = "ward.D2"), then work out the mean silhouette width (same as in the slides) for every k = 2, 3, 4, 5, 6. Which k do the numbers alone recommend? Does it match the known 3 species? If not, why, given how separated the three species looked in Exercise 1–2’s PCA?
Hint: You already know cutree(hc, k = 3) and silhouette(cl, d)[, "sil_width"] from the slides, just for one fixed k = 3. Here you repeat those same two steps for several values of k and take the mean silhouette width each time. sapply(2:6, function(k) ...) just runs whatever is inside it once for every k in 2:6 and collects the results into one vector. Fill in the blank below and you’re done:
This should make the exercise a lot easier — the only new thing is sapply() itself. Everything inside it is cutree() and silhouette(), which you already used with a fixed k = 3.
The best score is at k = 2, not 3. That’s a real limitation of silhouette-only model selection: it rewards the split with the widest gap, and that’s not always the split you actually care about biologically. But it’s not just a miss, either. It’s telling you something true about the data: setosa is very distinct from the other two species, and versicolor and virginica are similar enough to each other that splitting them barely helps the average silhouette at all. So the k = 2 result really means: one species stands apart, and the other two look alike. That’s the same overlap Exercise 5 counts directly later (most of its 26 mismatches sit right in this versicolor/virginica overlap). We still cut hc at k = 3 for the rest of this session, since here we do know the species count.
cl <-cutree(hc, k =3)
2.1.4 Exercise 4: Does the distance metric matter? 🟢 Core
Every distance in the slides was Euclidean (dist()’s default): straight- line distance through 4D space. dist(x, method = "manhattan") instead adds up the absolute differences on each variable, like walking a city grid instead of cutting straight across a block. The slides never used this one.
Re-cluster the scaled iris data with Manhattan distance and ward.D2 linkage, cut at k = 3, and cross-tabulate it against the Euclidean-based cl from Exercise 3. Do the two metrics agree?
Hint:dist() takes a method argument: "euclidean" (the default you’ve used so far) or "manhattan". Everything after that, hclust(), cutree(), table(), stays exactly the same. Only the distance matrix changes:
Both metrics agree almost perfectly on setosa. But around a third of the harder-to-separate versicolor/virginica points switch clusters depending on the metric. The choice of distance can matter just as much as the choice of linkage, especially where groups overlap.
2.1.5 Exercise 5: Capstone — which flowers does the clustering get “wrong”? 🔴 Challenge
The slides’ combined PCA + clustering plot shows clusters and species side by side, but it never asks which individual flowers disagree between the two. For each cluster, find its majority species (the species most of its flowers belong to). Then flag every flower whose own species isn’t that majority — these are the points the clustering effectively “gets wrong” compared to the known labels.
Build a PCA plot (PC1 vs PC2, colour = cluster), circle the misclassified points in black so they stand out, then save it as a PNG at 300 dpi and as an RDS file, the same way the slides showed for ggsave() and saveRDS().
Hint: Break it into two steps with dplyr. First, count(scores, cluster, species) counts how many flowers of each species sit in each cluster; grouping by cluster and taking slice_max(n, n = 1) keeps just the top row per cluster, its majority species. Second, join that back onto scores and compare each flower’s own species to its cluster’s majority with !=. Fill in the blanks:
cluster_majority <- scores %>%count(cluster, species) %>%group_by(____) %>%slice_max(n, n =1, with_ties =FALSE) %>%ungroup() %>%select(cluster, expected = species)scores <- scores %>%left_join(cluster_majority, by ="cluster") %>%mutate(mismatch = species ____ expected)
Then filter(scores, mismatch) pulls out just those rows for the extra black-circle geom_point() layer.
This should make the exercise a lot easier — the two blanks are cluster (what you’re grouping by) and != (not equal to). Everything else is count(), slice_max(), and left_join(), which you already know from dplyr.
Solution 5.
scores <-as_tibble(pca_all$x[, 1:2]) %>%mutate(cluster =factor(cl), species = iris$Species)cluster_majority <- scores %>%count(cluster, species) %>%group_by(cluster) %>%slice_max(n, n =1, with_ties =FALSE) %>%ungroup() %>%select(cluster, expected = species)scores <- scores %>%left_join(cluster_majority, by ="cluster") %>%mutate(mismatch = species != expected)count(scores, mismatch) # 26 mismatched out of 150
# A tibble: 2 × 2
mismatch n
<lgl> <int>
1 FALSE 124
2 TRUE 26
p_final <-ggplot(scores, aes(PC1, PC2, colour = cluster)) +geom_point(size =2.5) +geom_point(data =filter(scores, mismatch),shape =1, size =4.5, colour ="black", stroke =1 ) +labs(title ="Where clustering disagrees with the known species",subtitle ="circled: species differs from its cluster's majority species",colour ="Cluster" ) +theme_minimal(base_size =12)print(p_final)
count() + slice_max() finds, for each cluster, the species with the most members. That’s the “expected” label if the clustering were perfect. left_join() then attaches that expected label back onto every flower’s own row, so mismatch <- species != expected can compare them directly. 26 of the 150 flowers land in a cluster where their own species isn’t the majority, and almost all of them sit in the versicolor/virginica overlap. That’s the same region Exercise 3’s k-scan flagged as hard to separate, and where Exercise 2’s made-up flower landed too.
2.2 Exit ticket
What does predict() do differently from re-running prcomp() on new data?
Why might the “best” k by silhouette not match the known number of species?
Name one thing that changed when you swapped Euclidean for Manhattan distance.
Solution 6.
predict() reuses the centering, scaling, and rotation matrix that prcomp() already learned from the original 150 flowers, and applies them to the new row. So it lands on the same PC axes. Re-running prcomp() on the new data alone would give you a brand-new set of axes that you couldn’t compare to the old ones.
Silhouette rewards the split with the clearest gap. Setosa is very distinct from the other two species, but versicolor and virginica overlap a lot. So the numerically cleanest split is “setosa vs. everything else” (k = 2), even though there are really three species. That’s not just a limitation of the method, though — it’s also telling you something true: versicolor and virginica really are similar to each other, and setosa really is different from both.
Setosa’s cluster stayed the same either way, but quite a few of the harder-to-separate versicolor/virginica points switched clusters when Manhattan distance replaced Euclidean distance.
