library(MASS) library(dplyr) library(plotly) normv <- function(v) { sqrt(sum(v^2)) } # building intuition ---- # unit-vectors in 2D R <- 10 plot(0, 0, type = "n", asp = 1, xlab = "x", ylab = "y") theta <- seq(0, 2*pi, length.out = 500) lines(cos(theta), sin(theta)) for (r in 1:R) { w <- rnorm(2) w <- w / normv(w) arrows(0, 0, w[1], w[2], lwd = 2) } # intuition with data --- df <- read.csv("https://raw.githubusercontent.com/roualdes/data/refs/heads/master/penguins.csv") %>% select(body_mass_g, flipper_length_mm, bill_depth_mm) %>% na.omit() %>% mutate(c_bm = body_mass_g - mean(body_mass_g), c_fl = flipper_length_mm - mean(flipper_length_mm), c_bd = bill_depth_mm - mean(bill_depth_mm)) # multidimensional data plot(df$body_mass_g, df$flipper_length_mm, xlab = "body mass (g)", ylab = "flipper length (mm)") # centered data plot(df$c_bm, df$c_fl, xlab = "body mass (g, centered)", ylab = "flipper length (mm, centered)") X <- df %>% select(c_bm, c_fl) %>% as.matrix eig <- eigen(cov(X)) v <- sqrt(eig$values) w1 <- eig$vectors[,1] * v[1] w2 <- eig$vectors[,2] * v[2] plot(df$c_bm, df$c_fl, xlab = "body mass (g, centered)", ylab = "flipper length (mm, centered)", col = "grey70") arrows(0, 0, w1[1], w1[2], lwd = 2) arrows(0, 0, w2[1], w2[2], lwd = 2) # remove portion of data that corresponds to first eigenvector Xk <- X - X %*% eig$vectors[,1] %*% t(eig$vectors[,1]) # plot what remains and overlay second eigenvector plot(Xk[,1], Xk[,2], col = "grey70") arrows(0, 0, w2[1], w2[2], lwd = 2) # maybe even an attempt at intuition in 3d will be helpful? ---- df <- read.csv("https://roualdes.sfo3.digitaloceanspaces.com/data/abalone.csv") %>% mutate(c_d = diameter - mean(diameter), c_l = length - mean(length), c_h = height - mean(height)) X <- df %>% select(c_d, c_l, c_h) %>% as.matrix eig <- eigen(cov(X)) v <- sqrt(eig$values) w1 <- eig$vectors[,1] * v[1] w2 <- eig$vectors[,2] * v[2] w3 <- eig$vectors[,3] * v[3] plot_ly() %>% # points add_trace( type = "scatter3d", mode = "markers", x = df$c_d, y = df$c_l, z = df$c_h, marker = list(size = 5, color = "grey", opacity = 0.005), showlegend = FALSE ) %>% # vector add_trace( type = "scatter3d", mode = "lines", x = c(0, w1[1]), y = c(0, w1[2]), z = c(0, w1[3]), line = list(color = "black", width = 6), showlegend = FALSE ) %>% add_trace( type = "scatter3d", mode = "lines", x = c(0, w2[1]), y = c(0, w2[2]), z = c(0, w2[3]), line = list(color = "black", width = 6), showlegend = FALSE ) %>% add_trace( type = "scatter3d", mode = "lines", x = c(0, w3[1]), y = c(0, w3[2]), z = c(0, w3[3]), line = list(color = "black", width = 6), showlegend = FALSE ) %>% layout( scene = list( xaxis = list(title = "diameter (mm, centered)"), yaxis = list(title = "length (mm, centered)"), zaxis = list(title = "height (mm, centered)") ) ) # example ---- df <- read.csv("https://roualdes.sfo3.digitaloceanspaces.com/data/abalone.csv") %>% mutate(c_d = diameter - mean(diameter), c_l = length - mean(length), c_h = height - mean(height)) pca <- df %>% select(!rings & !age & where(is.numeric)) %>% as.matrix %>% prcomp dfpca <- data.frame(pca$x) dfpca$age <- df$age summary(lm(age ~ length + diameter + height + whole_weight + shucked_weight + viscera_weight + shell_weight, data = df)) summary(lm(age ~ PC1, data = dfpca)) # References ---- # 1. https://en.wikipedia.org/wiki/Principal_component_analysis # 2. GPT for plotly plot