¿Cómo realizar la regresión ortogonal (mínimos cuadrados totales) a través de PCA?

29

Siempre uso lm()en R para realizar una regresión lineal de en . Esa función devuelve un coeficiente tal queyxβ

y=βx.

Hoy aprendí sobre mínimos cuadrados totales y esa princomp()función (análisis de componentes principales, PCA) se puede utilizar para realizarla. Debería ser bueno para mí (más preciso). He hecho algunas pruebas usando princomp(), como:

r <- princomp( ~ x + y)

Mi problema es: ¿cómo interpretar sus resultados? ¿Cómo puedo obtener el coeficiente de regresión? Por "coeficiente" me refiero al número que tengo que usar para multiplicar el valor de para dar un número cercano a .βxy

Dail
fuente
Un momento chicos, estoy un poco confundido. mira: zoonek2.free.fr/UNIX/48_R/09.html Esto se llama PCA (Análisis de componentes principales, también conocido como "regresión ortogonal" o "sumas perpendiculares de cuadrados" o "mínimos cuadrados totales"), así que creo que estamos hablando acerca de TLS con princomp () ¿No?
Dail
No; Esas son dos cosas diferentes, vea el artículo de Wikipedia sobre PCA. El hecho de que se use aquí es un truco (no sé cuán exacto, pero lo comprobaré); Es por eso que la extracción compleja de coeficientes.
1
Una pregunta relacionada: stats.stackexchange.com/questions/2691/… y una publicación del blog está referenciada por una de las respuestas: cerebralmastication.com/2010/09/…
Jonathan

Respuestas:

48

Mínimos cuadrados ordinarios versus mínimos cuadrados totales

Consideremos primero el caso más simple de una sola variable predictiva (independiente) . Por simplicidad, deje que x e y estén centrados, es decir, la intersección siempre es cero. La diferencia entre la regresión OLS estándar y la regresión TLS "ortogonal" se muestra claramente en esta figura (adaptada por mí) de la respuesta más popular en el hilo más popular en PCA:xxy

OLS vs TLS

MCO se ajusta a la ecuación minimizando distancias al cuadrado entre los valores observados de Y y los valores predichos Y . TLS se ajusta a la misma ecuación al minimizar las distancias al cuadrado entre los puntos ( x , y ) y su proyección en la línea. En este caso más simple, la línea TLS es simplemente el primer componente principal de los datos 2D. Para encontrar β , haga PCA en los puntos , es decir, construya la matriz de covarianza y encuentre su primer vector propio ; entoncesy=βxyy^(x,y)β2 × 2 Σ v = ( v x , v y ) β = v y / v x(x,y)2×2Σv=(vx,vy)β=vy/vx.

En Matlab:

 v = pca([x y]);    //# x and y are centered column vectors
 beta = v(2,1)/v(1,1);

En R:

 v <- prcomp(cbind(x,y))$rotation
 beta <- v[2,1]/v[1,1]

Por cierto, esto dará paso a la pendiente correcta incluso si y no se centraron (porque las funciones integradas de PCA realizar automáticamente el centrado). Para recuperar la intersección, calcule .y β 0 = ˉ y - β ˉ xxyβ0=y¯βx¯

OLS vs. TLS, regresión múltiple

Dada una variable dependiente y muchas variables independientes (nuevamente, todas centradas para simplificar), la regresión se ajusta a una ecuaciónOLS hace el ajuste al minimizar los errores al cuadrado entre los valores observados de y los valores predichos . TLS hace el ajuste minimizando las distancias al cuadrado entre los puntos observados y los puntos más cercanos en el plano de regresión / hiperplano.x i y = β 1 x 1 + + β p x p . y y ( x , y ) R p + 1yxi

y=β1x1++βpxp.
yy^(X,y)Rpags+1

¡Tenga en cuenta que ya no hay una "línea de regresión"! La ecuación anterior especifica un hiperplano : es un plano 2D si hay dos predictores, un hiperplano 3D si hay tres predictores, etc. Por lo tanto, la solución anterior no funciona: no podemos obtener la solución TLS tomando solo la primera PC (que es una linea). Aún así, la solución se puede obtener fácilmente a través de PCA.

Como antes, PCA se realiza en puntos . Esta rendimientos vectores propios en columnas de . Los primeros vectores propios definen un hiperplano dimensional que necesitamos; el último (número ) vector propio es ortogonal a él. La cuestión es cómo transformar la base de dada por el primer vectores propios en los coeficientes.p + 1 V p p H p + 1 v p + 1 H p β(X,y)pags+1VpagspagsHpags+1vpags+1Hpagsβ

Observe que si establecemos para todo y solo , entonces , es decir, el vector se encuentra en el hiperplano . Por otro lado, sabemos que es ortogonal a él. Es decir, su producto punto debe ser cero:i k x k = 1 y = β k ( 0 , ... , 1 , ... , β k ) H H v p + 1 = ( v 1 , ... , v p + 1 )Xyo=0 0yokXk=1y^=βk

(0 0,...,1,...,βk)H
Hv k + β k v p + 1 = 0 β k = - v k / v p + 1 .
vpags+1=(v1,...,vpags+1)H
vk+βkvpags+1=0 0βk=-vk/ /vpags+1.

En Matlab:

 v = pca([X y]);    //# X is a centered n-times-p matrix, y is n-times-1 column vector
 beta = -v(1:end-1,end)/v(end,end);

En R:

 v <- prcomp(cbind(X,y))$rotation
 beta <- -v[-ncol(v),ncol(v)] / v[ncol(v),ncol(v)]

Nuevamente, esto producirá pendientes correctas incluso si e no estuvieran centradas (porque las funciones PCA integradas realizan automáticamente el centrado). Para recuperar la intersección, calcule .y β 0 = ˉ y - ˉ x βXyβ0 0=y¯-X¯β

Como comprobación de cordura, observe que esta solución coincide con la anterior en caso de que solo haya un único predictor . De hecho, entonces el espacio es 2D, y por lo tanto, dado que el primer vector propio PCA es ortogonal al segundo (último), .( x , y ) v ( 1 ) y / v ( 1 ) x = - v ( 2 ) x / v ( 2 ) yX(X,y)vy(1)/ /vX(1)=-vX(2)/ /vy(2)

Solución de forma cerrada para TLS

Sorprendentemente, resulta que hay una ecuación de forma cerrada para . El siguiente argumento está tomado del libro de Sabine van Huffel "Los mínimos cuadrados totales" (sección 2.3.2).β

Sea y las matrices de datos centradas. El último vector propio de PCA es un vector propio de la matriz de covarianza de con un valor propio . Si es un vector propio, entonces también lo es . Anotando la ecuación del vector propio: Xyvpags+1[Xy]σpags+12-vpags+1/ /vpags+1=(β-1)

(XXXyyXyy)(β-1)=σpags+12(β-1),
y calculando el producto a la izquierda, inmediatamente obtenemos que recuerda fuertemente la conocida expresión OLS
βTLS=(XX-σpags+12yo)-1Xy,
βOLS=(XX)-1Xy.

Regresión múltiple multivariante

La misma fórmula puede generalizarse al caso multivariante, pero incluso para definir qué hace TLS multivariante, requeriría algo de álgebra. Ver Wikipedia en TLS . La regresión OLS multivariada es equivalente a un grupo de regresiones OLS univariadas para cada variable dependiente, pero en el caso TLS no es así.

ameba dice Reinstate Monica
fuente
1
No conozco R, pero aún quería proporcionar fragmentos de R para referencia futura. Aquí hay muchas personas competentes en R. ¡Siéntase libre de editar mis fragmentos si es necesario! Gracias.
ameba dice Reinstate Monica
Buena publicación, pero si puedo preguntar, ¿qué garantiza el hecho de que el vector encuentra en el hiperplano? (0 0,...,1,...,βk)
JohnK
@JohnK, no estoy seguro de qué es exactamente lo que no está claro. Como escribí, deje que todo sea ​​igual a cero, aparte de . Luego, si conecta esto a , obtendrá . Entonces el punto encuentra en el hiperplano definido por la ecuación . x k = 1 y = β j x j y = β k1 = β k ( 0 , , 1 , β k ) y = β j x jXyoXk=1y=βjXjy=βk1=βk(0 0,...,1,...βk)y=βjXj
ameba dice Reinstate Monica
Parece que he leído mal esa parte, pero ahora está claro. Gracias por la aclaración también.
JohnK
2
En R, puede preferir "eigen (cov (cbind (x, y))) $ vectors" sobre "prcomp (cbind (x, y)) $ rotacion" porque el primero es mucho más rápido para vectores más grandes.
Thomas Browne
9

Basado en la ingenua implementación de GNU Octave que se encuentra aquí , algo como esto podría funcionar (grano de sal, es tarde).

tls <- function(A, b){

  n <- ncol(A)
  C <- cbind(A, b)

  V <- svd(C)$v
  VAB <- V[1:n, (n+1):ncol(V)]
  VBB <- V[(n+1):nrow(V), (n+1):ncol(V)]
  return(-VAB/VBB)
}
cashoes
fuente
4

princompestá ejecutando el análisis de componentes principales en lugar de la regresión total de mínimos cuadrados. Hasta donde yo sé, no hay función R ni paquete que haga TLS; a lo sumo hay regresión de Deming en MethComp .
Sin embargo, trate esto como una sugerencia de que lo más probable es que no valga la pena.


fuente
Pensé que Deming en el paquete MethComp era TLS: ¿cuál es la diferencia?
mark999
Debe darle la razón de errores en x e y; TLS puro optimiza esto.