Comparación visual de regresión y PCA

Question

Estoy tratando de perfeccionar un método para comparar regresión y PCA, inspirado en el blog Cerebral Mastication que también se ha discutido desde un ángulo diferente en SO. Antes de que me olvide, muchas gracias a JD Long y Josh Ulrich por la mayor parte de esto. Voy a usar esto en un curso el próximo semestre. Lo siento, esto es largo!

ACTUALIZACIÓN: Encontré un enfoque diferente que casi funciona (¡por favor, arréglenlo si pueden!). Lo publiqué en la parte inferior. ¡Un enfoque mucho más inteligente y más corto de lo que fui capaz de hacer!

Básicamente seguí los esquemas anteriores hasta cierto punto: Genere datos aleatorios, descubra la línea de mejor ajuste, dibuje los residuos. Esto se muestra en el segundo fragmento de código a continuación. Pero también busqué y escribí algunas funciones para dibujar líneas normales a una línea a través de un punto aleatorio (los puntos de datos en este caso). Creo que funcionan bien, y se muestran en First Code Chunk junto con la prueba de que funcionan.

Ahora, el segundo fragmento de código muestra todo en acción usando el mismo flujo que @JDLong y estoy agregando una imagen de la gráfica resultante. Los datos en negro, rojo es la regresión con residuos rosa, azul es la primera PC y el azul claro debería ser las normales, pero obviamente no lo son. Las funciones en First Code Chunk que dibujan estas normales parecen estar bien, pero algo no está bien con la demostración: creo que debo estar malinterpretando algo o transmitiendo los valores incorrectos. Mis normales vienen en horizontal, lo cual parece una pista útil (pero hasta ahora, no para mí). ¿Alguien puede ver lo que está mal aquí?

Gracias, esto me ha sido inquietado por un tiempo ...

Primer Código Chunk (Funciones para dibujar normales y prueba de que el trabajo):

##### The functions below are based very loosely on the citation at the end 

pointOnLineNearPoint <- function(Px, Py, slope, intercept) { 
    # Px, Py is the point to test, can be a vector. 
    # slope, intercept is the line to check distance. 

    Ax <- Px-10*diff(range(Px)) 
    Bx <- Px+10*diff(range(Px)) 
    Ay <- Ax * slope + intercept 
    By <- Bx * slope + intercept 
    pointOnLine(Px, Py, Ax, Ay, Bx, By) 
    } 

pointOnLine <- function(Px, Py, Ax, Ay, Bx, By) { 

    # This approach based upon comingstorm's answer on 
    # stackoverflow.com/questions/3120357/get-closest-point-to-a-line 
    # Vectorized by Bryan 

    PB <- data.frame(x = Px - Bx, y = Py - By) 
    AB <- data.frame(x = Ax - Bx, y = Ay - By) 
    PB <- as.matrix(PB) 
    AB <- as.matrix(AB) 
    k_raw <- k <- c() 
    for (n in 1:nrow(PB)) { 
     k_raw[n] <- (PB[n,] %*% AB[n,])/(AB[n,] %*% AB[n,]) 
     if (k_raw[n] < 0) { k[n] <- 0 
      } else { if (k_raw[n] > 1) k[n] <- 1 
       else k[n] <- k_raw[n] } 
     } 
    x = (k * Ax + (1 - k)* Bx) 
    y = (k * Ay + (1 - k)* By) 
    ans <- data.frame(x, y) 
    ans 
    } 

# The following proves that pointOnLineNearPoint 
# and pointOnLine work properly and accept vectors 

par(mar = c(4, 4, 4, 4)) # otherwise the plot is slightly distorted 
# and right angles don't appear as right angles 

m <- runif(1, -5, 5) 
b <- runif(1, -20, 20) 
plot(-20:20, -20:20, type = "n", xlab = "x values", ylab = "y values") 
abline(b, m) 

Px <- rnorm(10, 0, 4) 
Py <- rnorm(10, 0, 4) 

res <- pointOnLineNearPoint(Px, Py, m, b) 
points(Px, Py, col = "red") 
segments(Px, Py, res[,1], res[,2], col = "blue") 

##======================================================== 
## 
## Credits: 
## Theory by Paul Bourke http://local.wasp.uwa.edu.au/~pbourke/geometry/pointline/ 
## Based in part on C code by Damian Coventry Tuesday, 16 July 2002 
## Based on VBA code by Brandon Crosby 9-6-05 (2 dimensions) 
## With grateful thanks for answering our needs! 
## This is an R (http://www.r-project.org) implementation by Gregoire Thomas 7/11/08 
## 
##======================================================== 

segundo código Chunk (Plots las demostración):

set.seed(55) 
np <- 10 # number of data points 
x <- 1:np 
e <- rnorm(np, 0, 60) 
y <- 12 + 5 * x + e 

par(mar = c(4, 4, 4, 4)) # otherwise the plot is slightly distorted 

plot(x, y, main = "Regression minimizes the y-residuals & PCA the normals") 
yx.lm <- lm(y ~ x) 
lines(x, predict(yx.lm), col = "red", lwd = 2) 
segments(x, y, x, fitted(yx.lm), col = "pink") 

# pca "by hand" 
xyNorm <- cbind(x = x - mean(x), y = y - mean(y)) # mean centers 
xyCov <- cov(xyNorm) 
eigenValues <- eigen(xyCov)$values 
eigenVectors <- eigen(xyCov)$vectors 

# Add the first PC by denormalizing back to original coords: 

new.y <- (eigenVectors[2,1]/eigenVectors[1,1] * xyNorm[x]) + mean(y) 
lines(x, new.y, col = "blue", lwd = 2) 

# Now add the normals 

yx2.lm <- lm(new.y ~ x) # zero residuals: already a line 
res <- pointOnLineNearPoint(x, y, yx2.lm$coef[2], yx2.lm$coef[1]) 
points(res[,1], res[,2], col = "blue", pch = 20) # segments should end here 
segments(x, y, res[,1], res[,2], col = "lightblue1") # the normals 

############ ACTUALIZACIÓN

Más en 0.123.Encontré casi exactamente lo que quería. Pero, no funciona del todo (obviamente, solía funcionar). Aquí es un extracto de código de ese sitio que traza normales a la primera PC reflejan a través de un eje vertical:

set.seed(1) 
x <- rnorm(20) 
y <- x + rnorm(20) 
plot(y~x, asp = 1) 
r <- lm(y~x) 
abline(r, col='red') 

r <- princomp(cbind(x,y)) 
b <- r$loadings[2,1]/r$loadings[1,1] 
a <- r$center[2] - b * r$center[1] 
abline(a, b, col = "blue") 
title(main='Appears to use the reflection of PC1') 

u <- r$loadings 
# Projection onto the first axis 
p <- matrix(c(1,0,0,0), nrow=2) 
X <- rbind(x,y) 
X <- r$center + solve(u, p %*% u %*% (X - r$center)) 
segments(x, y, X[1,], X[2,] , col = "lightblue1") 

Y aquí está el resultado:

Answer 1

¡De acuerdo, tendré que responder mi propia pregunta! Después de leer más y comparar los métodos que la gente ha puesto en Internet, he resuelto el problema. No estoy seguro de poder decir claramente lo que "arreglé" porque pasé bastantes iteraciones. De todos modos, aquí está la trama y el código (MWE). Las funciones auxiliares están al final para mayor claridad.

# Comparison of Linear Regression & PCA 
# Generate sample data 

set.seed(39) # gives a decent-looking example 
np <- 10 # number of data points 
x <- -np:np 
e <- rnorm(length(x), 0, 10) 
y <- rnorm(1, 0, 2) * x + 3*rnorm(1, 0, 2) + e 

# Plot the main data & residuals 

plot(x, y, main = "Regression minimizes the y-residuals & PCA the normals", asp = 1) 
yx.lm <- lm(y ~ x) 
lines(x, predict(yx.lm), col = "red", lwd = 2) 
segments(x, y, x, fitted(yx.lm), col = "pink") 

# Now the PCA using built-in functions 
# rotation = loadings = eigenvectors 

r <- prcomp(cbind(x,y), retx = TRUE) 
b <- r$rotation[2,1]/r$rotation[1,1] # gets slope of loading/eigenvector 1 
a <- r$center[2] - b * r$center[1] 
abline(a, b, col = "blue") # Plot 1st PC 

# Plot normals to 1st PC 

X <- pointOnLineNearPoint(x, y, b, a) 
segments(x, y, X[,1], X[,2], col = "lightblue1") 

###### Needed Functions 

pointOnLineNearPoint <- function(Px, Py, slope, intercept) { 
    # Px, Py is the point to test, can be a vector. 
    # slope, intercept is the line to check distance. 

    Ax <- Px-10*diff(range(Px)) 
    Bx <- Px+10*diff(range(Px)) 
    Ay <- Ax * slope + intercept 
    By <- Bx * slope + intercept 
    pointOnLine(Px, Py, Ax, Ay, Bx, By) 
    } 

pointOnLine <- function(Px, Py, Ax, Ay, Bx, By) { 

    # This approach based upon comingstorm's answer on 
    # stackoverflow.com/questions/3120357/get-closest-point-to-a-line 
    # Vectorized by Bryan 

    PB <- data.frame(x = Px - Bx, y = Py - By) 
    AB <- data.frame(x = Ax - Bx, y = Ay - By) 
    PB <- as.matrix(PB) 
    AB <- as.matrix(AB) 
    k_raw <- k <- c() 
    for (n in 1:nrow(PB)) { 
     k_raw[n] <- (PB[n,] %*% AB[n,])/(AB[n,] %*% AB[n,]) 
     if (k_raw[n] < 0) { k[n] <- 0 
      } else { if (k_raw[n] > 1) k[n] <- 1 
       else k[n] <- k_raw[n] } 
     } 
    x = (k * Ax + (1 - k)* Bx) 
    y = (k * Ay + (1 - k)* By) 
    ans <- data.frame(x, y) 
    ans 
    } 

Answer 2

trate de cambiar esta línea de su código :

res <- pointOnLineNearPoint(x, y, yx2.lm$coef[2], yx2.lm$coef[1]) 

a

res <- pointOnLineNearPoint(x, new.y, yx2.lm$coef[2], yx2.lm$coef[1]) 

Así que está llamando a los valores y correctos.

Answer 3

En Vincent Zoonekynd's code, cambie la línea u <- r$loadings a u <- solve(r$loadings). En la segunda instancia de solve(), las puntuaciones del componente predicho a lo largo del primer eje principal (es decir, la matriz de puntuaciones pronosticadas con las puntuaciones de los componentes pronosticados establecidos en cero) deben multiplicarse por inversa de las cargas/vectores propios. Al multiplicar los datos por las cargas se obtienen puntajes pronosticados; dividir los puntajes pronosticados por las cargas dan datos. Espero que ayude.