Activity 3.2 - PC concepts

SUBMISSION INSTRUCTIONS

Render to html
Publish your html to RPubs
Submit a link to your published solutions

Question 1

Make sure to use SVD to find the principal components

data set we will analyze for this question are on the 10 events in the Men’s 2024 Olympic Decathlon in Paris.

decathlon <- read.csv('Data/mens_decathlon_paris2024.csv')
head(decathlon)

         Athlete  Medal      Nation Overall X100m LongJump ShotPut HighJump
1   Markus Rooth   Gold      Norway    8796 10.71     7.80   15.25     1.99
2 Leo Neugebauer Silver     Germany    8748 10.67     7.98   16.55     2.05
3  Lindon Victor Bronze     Grenada    8711 10.56     7.48   15.71     2.02
4    Sven Roosen   None Netherlands    8607 10.52     7.56   15.10     1.87
5  Janek Õiglane   None     Estonia    8572 10.89     7.25   14.58     1.99
6   Johannes Erm   None     Estonia    8569 10.64     7.66   14.61     2.08
  X400m X110mHurdle Discus PoleVault Javelin X1500m
1 47.69       14.25  49.80       5.3   66.87  279.6
2 47.70       14.51  53.33       5.0   56.64  284.7
3 47.84       14.62  53.91       4.9   68.22  283.5
4 46.40       13.99  46.88       4.7   63.72  258.5
5 48.02       14.45  43.49       5.3   71.89  265.6
6 47.19       14.35  46.29       4.6   59.58  259.7

For the purposes of this question, assume we have 10-dimensional data - that is, ignore the Overall column.

A)

Explain why we need to scale this data set before performing PCA

Because, the variables or columns are not on the same scale. X1500m is on a much larger scale compared to HighJump.

B)

Use svd() to find the first 2 principal component scores and their loadings. Full credit will only be granted if you use the svd() ingredients u, d, and v. What percent of the overall variability do the first two PCs explain?

library(tidyverse)

Warning: package 'tidyverse' was built under R version 4.4.3

Warning: package 'tibble' was built under R version 4.4.3

Warning: package 'tidyr' was built under R version 4.4.3

Warning: package 'readr' was built under R version 4.4.3

Warning: package 'purrr' was built under R version 4.4.3

Warning: package 'dplyr' was built under R version 4.4.3

Warning: package 'stringr' was built under R version 4.4.3

Warning: package 'forcats' was built under R version 4.4.3

Warning: package 'lubridate' was built under R version 4.4.3

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.2     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.4     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

selected_data <- decathlon %>% select(X100m:X1500m)
X <- scale(selected_data)
X_svd <- svd(X)

U <- X_svd$u
D <- diag(X_svd$d)
V <- X_svd$v

PCs <- U %*% D
PCs_2 <- PCs[, 1:2] 
head(PCs_2)

           [,1]        [,2]
[1,] -0.9696335  0.90753969
[2,] -2.1338950  2.24014162
[3,] -0.7323012  2.05468273
[4,] -1.4369409 -0.05634227
[5,]  1.0587944 -0.11932168
[6,] -1.1274914 -0.28469476

(PCs_2
 %>% data.frame()
 %>% summarize(across(everything(), var))
 %>% mutate(across(everything(), \(x) x/10))
)

         X1        X2
1 0.2913494 0.1928006

PC1 alone explains 29.1% of the total variability in these 10 variables and PC2 additional 19.2 % of the total variability.

C)

Find and print the loadings. Based on the loadings alone, if the first two PCs are plotted in a 2D plane as shown below, which of the four quadrants will the medalists be in? Explain your reasoning.

loadings <- X_svd$v
rownames(loadings) <- colnames(selected_data)
loadings[,1:2]

                    [,1]        [,2]
X100m        0.446612410  0.01678592
LongJump    -0.493037935 -0.06565628
ShotPut     -0.309431542  0.53134100
HighJump    -0.150386291 -0.05665757
X400m        0.485973869  0.14847674
X110mHurdle  0.345270647  0.31685991
Discus      -0.139145881  0.59570522
PoleVault    0.019829831  0.48083261
Javelin      0.252821246  0.03137029
X1500m       0.005588625 -0.01948870

Running events have positive loadings and field events have negative loadings. PC1 separates runners form throwers or jumpers.

Since we are looking for smaller number or time for runners, the PC1 score will be more negative. Fast or top runners get negative PC1 score.

PC2 represents throwing and upper-body power. Athletes who are good at those will have positive or high PC2 score.

good in running (negative PC1) and good in throwing (positive PC2)

Medalists will be in the 2nd quadrant.

D)

Add the PCs to the decathlon data set and create a scatterplot of these PCs, with the points labeled by the athletes’ names. Color-code the points on whether or not the athlete was a medalist. Use the ggrepel package for better labeling. Verify that your intuition from C) is correct.

decathlon$PC1 <- PCs_2[,1]
decathlon$PC2 <- PCs_2[,2]

library(ggrepel)

Warning: package 'ggrepel' was built under R version 4.4.3

ggplot(decathlon, aes(x = PC1, y = PC2, color = Medal, label = Athlete)) +
  geom_point(size = 3) +
  geom_text_repel(size = 3, show.legend = FALSE) +
  labs(
    x = "(PC1)",
    y = "(PC2)",
    color = "Medal"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, face = "bold")
   )

Yes, it is correct the medalists (Leo, Lindon and Markus) appear in the 2nd quadrant.

E)

Canadian Damian Warner won the gold medal in the decathlon in the 2020 Tokyo games. He began the 2024 decathlon but bowed out after three straight missed pole vault attempts.

These are his results in the 10 events in 2020:

warner <- c(10.12, 8.24, 14.8, 2.02, 47.48, 13.46, 48.67, 4.9, 63.44, 271.08)

Would this have won a medal if it had happened in 2024? To answer this, we will compute his PCs with respect to the 2024 athletes and add it to the plot to see where his 2020 gold-medal performance compares to the 2024 athletes. To do this:

Find the mean vector from the 2024 athletes. Call it mean_vec_24.
Find the standard deviation vector from the 2024 athletes. Call it sd_vec_24.
Standardize Warner’s 2020 results with respect to the 2024 athletes: (warner-mean_vec_24)/sd_vec_24
Find Warner’s PC coordinates using the 2024 loadings.
Add his point to the scatterplot.

Do you think his 2020 performance would have won a medal if it had happened in 2024?

mean_vec_24 <- apply(selected_data, 2, mean)
sd_vec_24 <- apply(selected_data, 2, sd)

warner_std <- (warner - mean_vec_24) / sd_vec_24

warner_pc <- as.numeric(warner_std %*% loadings[, 1:2])

warner_point <- data.frame(
  Athlete = "Damian Warner (2020)",
  PC1= warner_pc[1],
  PC2= warner_pc[2],
  Medalist = "2020 Gold"
)

ggplot(decathlon, aes(x = PC1, y = PC2, color = Medal, label = Athlete)) +
  geom_point(size = 3) +
  geom_text_repel(size = 3, show.legend = FALSE) +
  geom_point(data = warner_point, aes(x = PC1, y = PC2),
             color = "black", size = 5, shape = 8) +
  geom_text_repel(data = warner_point, aes(label = Athlete), size = 3.5,
                  color = "black") +
  theme_minimal()

Damian Warner’s 2020 performance would not have won a medal if it had happened in 2024. His performance is located in the third quadrant, not the second where the 2024 medalists are.

Question 2

Below is a screenshot of a conversation between me and chatbot Claude:

After looking at the graphs, I grew skeptical. So I said:

Behold, Claude’s three data sets which I’ve called claudeA, claudeB, and claudeC:

claudeA <- read.csv('Data/claude_dataA.csv')
claudeB <- read.csv('Data/claude_dataB.csv')
claudeC <- read.csv('Data/claude_dataC.csv')

Each data set has an X and a Y column which represent 2-dimensional variables that we need to rotate.

A)

Scale each data set and plot them side-by-side using the patchwork package. Make sure the aspect ratio of each graph is 1 (i.e., make the height and width of each graph equal). At this point, explain why you think I was skeptical. Specifically, do you think the percent variability explained by the first PC of each data set appears to exceed or fall short of the variability I asked it to?

scaled_claudeA <- scale(claudeA)
scaled_claudeB <- scale(claudeB)
scaled_claudeC <- scale(claudeC)

library(ggplot2)

plotA <- ggplot(as.data.frame(scaled_claudeA), aes(X, Y)) +
  geom_point(color = "blue") +
  coord_fixed(ratio = 1) +
  ggtitle("Claude A")

plotB <- ggplot(as.data.frame(scaled_claudeB), aes(X, Y)) +
  geom_point(color = "red") +
  coord_fixed(ratio = 1) +
  ggtitle("Claude B")

plotC <- ggplot(as.data.frame(scaled_claudeC), aes(X, Y)) +
  geom_point(color = "green") +
  coord_fixed(ratio = 1) +
  ggtitle("Claude C")

library(patchwork)

Warning: package 'patchwork' was built under R version 4.4.3

(plotA | plotB | plotC)

Data set A has a negative correlation and moderate scatter around the trend line. Data set B has positive correlation and and looks like percent variability explained by the first PC is less compared to data set A. Data et C has a strong negative correlation and the percent variability explained by the first PC is higher than the first two data sets.

B)

Use SVD to find the first PC for each data set, and find the actual percent of total variability explained by each PC using aggregation methods.

svdA <- svd(scaled_claudeA)
svdB <- svd(scaled_claudeB)
svdC <- svd(scaled_claudeC)

scoresA <- svdA$u %*% diag(svdA$d)
scoresB <- svdB$u %*% diag(svdB$d)
scoresC <- svdC$u %*% diag(svdC$d)

PC1_a <- scoresA[, 1] 
PC1_b <- scoresB[, 1] 
PC1_c <- scoresC[, 1]

(PC1_a
 %>% data.frame()
 %>% summarize(across(everything(), var))
 %>% mutate(across(everything(), \(x) x/2))
)

          .
1 0.9756201

(PC1_b
 %>% data.frame()
 %>% summarize(across(everything(), var))
 %>% mutate(across(everything(), \(x) x/2))
)

          .
1 0.9591807

(PC1_c
 %>% data.frame()
 %>% summarize(across(everything(), var))
 %>% mutate(across(everything(), \(x) x/2))
)

          .
1 0.9949096

PC_a = 97.56 %

PC_b = 95.91 %

PC_c= 99.4%

percent variability explained by the first PC is higest for data set C then A and B.