Activity 3.2 - PC concepts

SUBMISSION INSTRUCTIONS

  1. Render to html
  2. Publish your html to Posit Connect
  3. Submit a link to your published solutions

Question 1

The 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.

If not scaled, events that have higher scores will naturally have higher variance therefore carry more weight. Specifically, the 1500m race is the only one with scores in the triple digits, so that will have the most impact. The units are also different. So when I say ‘scores’ sometimes that means distance and sometimes that means time, which is a reason to scale the data.

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)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.1     ✔ readr     2.2.0
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.3     ✔ tibble    3.3.1
✔ lubridate 1.9.5     ✔ tidyr     1.3.2
✔ purrr     1.2.2     
── 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
decathlon_scaled <- (decathlon
                     %>% select(X100m:X1500m)
                     %>% scale()
)

decathlon_svd <- svd(decathlon_scaled)
U <- decathlon_svd$u
D <- diag(decathlon_svd$d)
V <- decathlon_svd$v

PCs <- U %*% D
head(PCs[, 1: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
loadings <- V 
rownames(loadings) <- colnames(decathlon_scaled)
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
(PCs
  %>% data.frame()
  %>% summarize(across(everything(), var))
  %>% mutate(across(everything(), \(x) x/10))
)
         X1        X2       X3        X4        X5         X6         X7
1 0.2913494 0.1928006 0.146543 0.1404117 0.0952215 0.05535841 0.04806271
          X8         X9         X10
1 0.01403932 0.01222556 0.003987768

The first two PCs explain about 48.41% of the total variability (PC1: 29.13%, PC2: 19.28%).

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 <- V
rownames(loadings) <- colnames(decathlon_scaled)
loadings[, 1:2] %>% round(3)
              [,1]   [,2]
X100m        0.447  0.017
LongJump    -0.493 -0.066
ShotPut     -0.309  0.531
HighJump    -0.150 -0.057
X400m        0.486  0.148
X110mHurdle  0.345  0.317
Discus      -0.139  0.596
PoleVault    0.020  0.481
Javelin      0.253  0.031
X1500m       0.006 -0.019

The medalists should be in Quadrant II (negative PC1 and positive PC2).

PC1: The biggest positive loadings are for the running events (100m, 400m, and 110m hurdles), while the long jump, shot put, high jump, and discus have negative loadings. Since lower times are better in running events, faster runners will have lower values, which pushes PC1 down. The same goes for athletes who do well in the jumping and throwing events. So overall, the better all-around athletes should have negative PC1 scores.

PC2: The biggest loadings are for discus (0.596), shot put (0.531), and pole vault (0.481). This means PC2 is mostly picking up performance in those events. Athletes who do well in these events should have higher PC2 scores. Medalists would likely have positive PC2 scores.

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.

library(ggrepel)

decathlon_pcs <- (decathlon
                  %>% bind_cols(data.frame(PCs[, 1:2]))
                  %>% rename(PC1 = X1, PC2 = X2)
                  %>% mutate(Medalist = ifelse(Medal == 'None',
                                               'No', 'Yes'))
)

ggplot(data = decathlon_pcs) +
  geom_hline(yintercept = 0, color = 'gray') +
  geom_vline(xintercept = 0, color = 'gray') +
  geom_point(aes(x = PC1, y = PC2, color = Medalist)) +
  geom_text_repel(aes(x = PC1, y = PC2,
                      label = Athlete, color = Medalist),
                  size = 3, show.legend = FALSE) +
  scale_color_manual(values = c('No' = 'black',
                                'Yes' = 'orange')) +
  labs(x = 'PC1', y = 'PC2') +
  theme_classic(base_size = 14)

Correct….medalists in Q2

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 <- (decathlon
                %>% select(X100m:X1500m)
                %>% apply(2, mean)
)

sd_vec_24 <- (decathlon
              %>% select(X100m:X1500m)
              %>% apply(2, sd)
)

warner_scaled <- (warner - mean_vec_24)/sd_vec_24

warner_pcs <- (warner_scaled %*% V[, 1:2]
               %>% data.frame()
               %>% rename(PC1 = X1, PC2 = X2)
               %>% mutate(Athlete = 'Damian Warner (2020)')
)
warner_pcs

ggplot(data = decathlon_pcs) +
  geom_hline(yintercept = 0, color = 'gray') +
  geom_vline(xintercept = 0, color = 'gray') +
  geom_point(aes(x = PC1, y = PC2, color = Medalist)) +
  geom_text_repel(aes(x = PC1, y = PC2,
                      label = Athlete, color = Medalist),
                  size = 3, show.legend = FALSE) +
  geom_point(aes(x = PC1, y = PC2),
             data = warner_pcs, color = 'red', size = 3) +
  geom_text_repel(aes(x = PC1, y = PC2, label = Athlete),
                  data = warner_pcs, color = 'red', size = 3) +
  scale_color_manual(values = c('No' = 'black',
                                'Yes' = 'orange')) +
  labs(x = 'PC1', y = 'PC2') +
  theme_classic(base_size = 14)

        PC1        PC2              Athlete
1 -3.422699 -0.5092531 Damian Warner (2020)

Based on our analysis Damian Warner’s 2020 performance would not have won a medal if it would’ve happened in 2024. His performance falls into Quadrant III while all medalists (and other non medalists) fall into Quadrant II.

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?

library(patchwork)

claudeA_scaled <- (claudeA
                   %>% scale()
                   %>% data.frame()
)
claudeB_scaled <- claudeB %>% scale() %>% data.frame()
claudeC_scaled <- claudeC %>% scale() %>% data.frame()

make_plot <- \(df, mytitle) {
  ggplot(data = df) +
    geom_point(aes(x = X, y = Y)) +
    xlim(c(-4, 4)) + ylim(c(-4, 4)) +
    ggtitle(mytitle) +
    theme_classic(base_size = 14) +
    theme(aspect.ratio = 1)
}

(make_plot(claudeA_scaled, 'A') +
    make_plot(claudeB_scaled, 'B') +
    make_plot(claudeC_scaled, 'C')
)

The 3 plots claim to be pretty different (55% variance, 75% variance, and 90% variance). But they all present as being very similar variance. Specifically plot A should look like a round cloud, when it actually it shows a pretty strong relationship. Also, the directions of the correlations don’t match Claude’s descriptions, since A is negative and B is positive.

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.

get_PCs <- \(df_scaled) {
  df_svd <- svd(as.matrix(df_scaled))
  U <- df_svd$u
  D <- diag(df_svd$d)
  U %*% D
}

PCs_A <- get_PCs(claudeA_scaled)
PCs_B <- get_PCs(claudeB_scaled)
PCs_C <- get_PCs(claudeC_scaled)

head(PCs_A[, 1])
[1] -2.2900182 -1.9202720 -1.9400319 -1.4404504 -1.4602102 -0.9606287
head(PCs_B[, 1])
[1] -2.0192225 -1.7047641 -1.8808246 -1.2802303 -1.4562907 -0.8556964
head(PCs_C[, 1])
[1] -2.270557 -1.933711 -1.914344 -1.486791 -1.467423 -1.039870
(list(A = PCs_A, B = PCs_B, C = PCs_C)
  %>% map(data.frame)
  %>% bind_rows(.id = 'Dataset')
  %>% group_by(Dataset)
  %>% summarize(across(everything(), var))
  %>% mutate(across(X1:X2, \(x) x/2))
)
# A tibble: 3 × 3
  Dataset    X1      X2
  <chr>   <dbl>   <dbl>
1 A       0.976 0.0244 
2 B       0.959 0.0408 
3 C       0.995 0.00509

Actual:

A: 97.6%

B: 95.9%

C: 99.5%