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.

library(tidyverse)
library(ggplot2)
library(ggrepel)
library(dplyr)
library(patchwork)
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.

We must scale the data because the 10 decathlon events are recorded on fundamentally different scales and in different units. If we do not scale events with numerically large values and high variance such as the 1500m run would dominate the singular values and principal components, simply because of their measurement unit rather than their athletic significance.

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?

event_data <- decathlon[, 5:14]

Z <- scale(event_data)
n <- nrow(Z)


svd_decomp <- svd(Z)

U <- svd_decomp$u
D <- svd_decomp$d
V <- svd_decomp$v


pc_scores <- U[, 1:2] %*% diag(D[1:2])
colnames(pc_scores) <- c("PC1", "PC2")
 

loadings_pc1_pc2 <- V[, 1:2]
rownames(loadings_pc1_pc2) <- colnames(event_data)
colnames(loadings_pc1_pc2) <- c("PC1", "PC2")

pc_variance <- (D^2) / (n - 1)
pct_explained <- (pc_variance / sum(pc_variance)) * 100

pct_pc1 <- pct_explained[1]
pct_pc2 <- pct_explained[2]
pct_total_first_two <- pct_pc1 + pct_pc2

cat(sprintf("PC1 explains: %.2f%%\n", pct_pc1))
PC1 explains: 29.13%
cat(sprintf("PC2 explains: %.2f%%\n", pct_pc2))
PC2 explains: 19.28%
cat(sprintf("Together, PC1 and PC2 explain: %.2f%%\n", pct_total_first_two))
Together, PC1 and PC2 explain: 48.41%

  • PC1 variance explained: \(\approx \mathbf{29.13\%}\)

  • PC2 variance explained: \(\approx \mathbf{19.28\%}\)

  • Total variability explained by the first two PCs: \(\approx \mathbf{48.41\%}\)

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.

print(round(loadings_pc1_pc2, 4))
                PC1     PC2
X100m        0.4466  0.0168
LongJump    -0.4930 -0.0657
ShotPut     -0.3094  0.5313
HighJump    -0.1504 -0.0567
X400m        0.4860  0.1485
X110mHurdle  0.3453  0.3169
Discus      -0.1391  0.5957
PoleVault    0.0198  0.4808
Javelin      0.2528  0.0314
X1500m       0.0056 -0.0195

Medalists will be in Quadrant II (\(\text{PC1} < 0, \text{PC2} > 0\)):

  • PC1 (Overall Athletic Performance): Running events (100m, 400m, 110m hurdles) have positive loadings (\(+0.45, +0.49, +0.35\)), meaning faster times (negative \(z\)-scores) produce negative PC1 values. Jumping and throwing events (Long Jump, Shot Put) have negative loadings (\(-0.49, -0.31\)), meaning superior distances (positive \(z\)-scores) also produce negative PC1 values. Top decathletes consistently have \(\text{PC1} < 0\).

  • PC2 (Power & Throwing Strength): PC2 is dominated by positive loadings on heavy field events: Discus (\(+0.60\)), Shot Put (\(+0.53\)), and Pole Vault (\(+0.48\)). Elite medalists score high in these power events, producing positive \(z\)-scores and thus \(\text{PC2} > 0\).

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.

# Add PC scores to decathlon dataset
decathlon$PC1 <- pc_scores[, 1]
decathlon$PC2 <- pc_scores[, 2]
decathlon$is_medalist <- ifelse(!is.na(decathlon$Medal), "Medalist", "Non-Medalist")

# Plot
ggplot(decathlon, aes(x = PC1, y = PC2, color = is_medalist)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  geom_point(size = 3) +
  geom_text_repel(aes(label = Athlete), max.overlaps = 15) +
  scale_color_manual(values = c("Medalist" = "gold3", "Non-Medalist" = "steelblue")) +
  theme_minimal() +
  labs(
    title = "Paris 2024 Men's Decathlon - PCA via SVD",
    x = "PC1 (29.13% Variance Explained)",
    y = "PC2 (19.28% Variance Explained)",
    color = "Status"
  )

The three medalists:

  • Markus Rooth (Gold): \(\text{PC1} = -0.97, \text{PC2} = +0.91\)

  • Leo Neugebauer (Silver): \(\text{PC1} = -2.13, \text{PC2} = +2.24\)

  • Lindon Victor (Bronze): \(\text{PC1} = -0.73, \text{PC2} = +2.05\)

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)

# 1. Mean and SD vectors from the 2024 athletes
mean_vec_24 <- colMeans(event_data)
sd_vec_24 <- apply(event_data, 2, sd)

# 2. Standardize Warner's 2020 marks using 2024 scaling parameters
warner_std <- (warner - mean_vec_24) / sd_vec_24

# 3. Project Warner into the 2024 PC coordinate space using loadings (V)
warner_pcs <- warner_std %*% loadings_pc1_pc2
cat(sprintf("Damian Warner Coordinates: PC1 = %.4f, PC2 = %.4f\n", warner_pcs[1], warner_pcs[2]))
Damian Warner Coordinates: PC1 = -3.4227, PC2 = -0.5093
# 4. Add point to scatterplot
warner_df <- data.frame(
  Athlete = "Damian Warner (2020)",
  Medal = "2020 Gold",
  Nation = "Canada",
  Overall = 9018,
  PC1 = warner_pcs[1],
  PC2 = warner_pcs[2],
  is_medalist = "Warner (2020)"
)

combined_plot_data <- bind_rows(decathlon, warner_df)

ggplot(combined_plot_data, aes(x = PC1, y = PC2, color = is_medalist)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  geom_point(size = 3) +
  geom_text_repel(aes(label = Athlete), max.overlaps = 20) +
  scale_color_manual(values = c(
    "Medalist" = "gold3", 
    "Non-Medalist" = "steelblue", 
    "Warner (2020)" = "firebrick"
  )) +
  theme_minimal() +
  labs(
    title = "Damian Warner (2020) vs. Paris 2024 Decathletes",
    x = "PC1 (29.13% Variance Explained)",
    y = "PC2 (19.28% Variance Explained)"
  )

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?

Yes. Warner’s \(\text{PC1}\) score of \(-3.42\) places him further to the extreme negative side of PC1 than any athlete in the 2024 Paris field. Because negative PC1 measures elite performance across the primary decathlon events (especially his world-class \(10.12\text{s}\) in 100m, \(8.24\text{m}\) in Long Jump, and \(13.46\text{s}\) in 110m hurdles), his overall athletic dominance exceeds the 2024 medalists.

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?

# Scale each dataset
scale_A <- as.data.frame(scale(claudeA))
scale_B <- as.data.frame(scale(claudeB))
scale_C <- as.data.frame(scale(claudeC))

# Create plots with 1:1 aspect ratio
pA <- ggplot(scale_A, aes(x = X, y = Y)) +
  geom_point(color = "darkred", alpha = 0.7) +
  coord_fixed(ratio = 1) +
  theme_minimal() +
  labs(title = "Dataset A (Claude claimed 55%)", x = "Scaled X", y = "Scaled Y")

pB <- ggplot(scale_B, aes(x = X, y = Y)) +
  geom_point(color = "darkblue", alpha = 0.7) +
  coord_fixed(ratio = 1) +
  theme_minimal() +
  labs(title = "Dataset B (Claude claimed 75%)", x = "Scaled X", y = "Scaled Y")

pC <- ggplot(scale_C, aes(x = X, y = Y)) +
  geom_point(color = "darkgreen", alpha = 0.7) +
  coord_fixed(ratio = 1) +
  theme_minimal() +
  labs(title = "Dataset C (Claude claimed 90%)", x = "Scaled X", y = "Scaled Y")

# Combine
pA + pB + pC

The skepticism is justified because the actual variance explained by PC1 exceeds the claimed values (\(55\%\), \(75\%\), and \(90\%\)):

  • For the two standardized variables, a PC1 explaining only \(55\%\) variance requires an almost circular, uncorrelated scatter cloud (\(r \approx 0.10\)). Instead, Dataset A shows an extremely tight, near-perfect linear line (\(r = -0.951\), yielding \(97.56\%\) variance on PC1).

  • Datasets B and C likewise display almost zero orthogonal scatter (\(r = 0.918\) and \(r = -0.990\), yielding \(95.92\%\) and \(99.49\%\) on PC1, respectively).

  • The claim understated the linear correlation and variance explained in all three datasets.

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.

calc_svd_var <- function(df, label) {
  # Standardize
  Z <- scale(df)
  n <- nrow(Z)
  
  # SVD
  s <- svd(Z)
  

  var_pcs <- (s$d^2) / (n - 1)
  pct_var <- (var_pcs / sum(var_pcs)) * 100
  
  tibble(
    Dataset = label,
    Singular_Value_1 = s$d[1],
    Singular_Value_2 = s$d[2],
    Actual_PC1_Pct = pct_var[1],
    Claude_Claimed_Pct = case_when(
      label == "Dataset A" ~ 55,
      label == "Dataset B" ~ 75,
      label == "Dataset C" ~ 90
    )
  )
}

results <- bind_rows(
  calc_svd_var(claudeA, "Dataset A"),
  calc_svd_var(claudeB, "Dataset B"),
  calc_svd_var(claudeC, "Dataset C")
)

print(results)
# A tibble: 3 × 5
  Dataset   Singular_Value_1 Singular_Value_2 Actual_PC1_Pct Claude_Claimed_Pct
  <chr>                <dbl>            <dbl>          <dbl>              <dbl>
1 Dataset A             8.72            1.38            97.6                 55
2 Dataset B             8.65            1.78            95.9                 75
3 Dataset C             8.81            0.630           99.5                 90