¿Qué hace realmente la función R "poly"?

Question

He leído la página del manual ?poly (que admito que no compliqué por completo) y también leí la descripción de la función en el libro Introducción al aprendizaje estadístico .

Mi comprensión actual es que una llamada a poly(horsepower, 2) debería ser equivalente a escribir horsepower + I(horsepower^2) . Sin embargo, esto parece estar en contradicción con el resultado del siguiente código:

library(ISLR)
summary(lm(mpg~poly(horsepower,2), data=Auto))$coef
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) 23.44592 0.2209163 106.13030 2.752212e-289
#poly(horsepower, 2)1 -120.13774 4.3739206 -27.46683 4.169400e-93
#poly(horsepower, 2)2 44.08953 4.3739206 10.08009 2.196340e-21
summary(lm(mpg~horsepower+I(horsepower^2), data=Auto))$coef
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) 56.900099702 1.8004268063 31.60367 1.740911e-109
#horsepower -0.466189630 0.0311246171 -14.97816 2.289429e-40
#I(horsepower^2) 0.001230536 0.0001220759 10.08009 2.196340e-21

Mi pregunta es, ¿por qué las salidas no coinciden, y qué está realmente haciendo poly ?

Answer 1

Cuando se introducen términos polinomiales en un modelo estadístico, la motivación habitual es determinar si la respuesta es "curva" y si la curvatura es "significativa" cuando se agrega ese término. El resultado de lanzar en un +I(x^2) términos es que las desviaciones menores pueden ser "magnificadas" por el proceso de ajuste dependiendo de su ubicación, y malinterpretadas como debidas al término de curvatura cuando solo eran fluctuaciones en un extremo u otro del rango de datos. Esto resulta en la asignación inapropiada de declaraciones de "significado".

Si solo agregas un término al cuadrado con I(x^2) , necesariamente también estará altamente correlacionado con x al menos en el dominio donde x > 0 . En su lugar, utilizando: poly(x,2) crea un conjunto "curvado" de variables donde el término lineal no está tan correlacionado con x, y donde la curvatura es aproximadamente la misma en todo el rango de datos. (Si desea leer sobre la teoría estadística, busque "polinomios ortogonales".) Simplemente escriba poly(1:10, 2) y observe las dos columnas.

poly(1:10, 2)
1 2
[1,] -0.49543369 0.52223297
[2,] -0.38533732 0.17407766
[3,] -0.27524094 -0.08703883
[4,] -0.16514456 -0.26111648
[5,] -0.05504819 -0.34815531
[6,] 0.05504819 -0.34815531
[7,] 0.16514456 -0.26111648
[8,] 0.27524094 -0.08703883
[9,] 0.38533732 0.17407766
[10,] 0.49543369 0.52223297
attr(,"degree")
[1] 1 2
attr(,"coefs")
attr(,"coefs")$alpha
[1] 5.5 5.5
attr(,"coefs")$norm2
[1] 1.0 10.0 82.5 528.0
attr(,"class")
[1] "poly" "matrix"

El término "cuadrático" se centra en 5.5 y el término lineal se ha desplazado hacia abajo, por lo que es 0 en el mismo punto x (con el término implícito (Intercept) en el modelo del que se depende para desplazar todo de vuelta en el momento en que se predicen pedido.)

Answer 2

Para obtener polinomios comunes como en la pregunta use raw = TRUE . Desafortunadamente, hay un aspecto indeseable con los polinomios ordinarios en la regresión. Si ajustamos un cuadrático, por ejemplo, y luego un cúbico, los coeficientes de orden inferior del cubo son todos diferentes que para el cuadrático, es decir, 56.900099702, -0.466189630, 0.001230536 para el cuadrático frente a 6.068478e + 01, -5.688501e-01, 2.079011e-03 después de volver a montar con un cúbico a continuación.

library(ISLR)
fm2raw <- lm(mpg ~ poly(horsepower, 2, raw = TRUE), Auto)
cbind(coef(fm2raw))
## [,1]
## (Intercept) 56.900099702
## poly(horsepower, 2, raw = TRUE)1 -0.466189630
## poly(horsepower, 2, raw = TRUE)2 0.001230536
fm3raw <- lm(mpg ~ poly(horsepower, 3, raw = TRUE), Auto)
cbind(coef(fm3raw))
## [,1]
## (Intercept) 6.068478e+01
## poly(horsepower, 3, raw = TRUE)1 -5.688501e-01
## poly(horsepower, 3, raw = TRUE)2 2.079011e-03
## poly(horsepower, 3, raw = TRUE)3 -2.146626e-06

Lo que realmente nos gustaría es agregar el término cúbico de tal forma que los coeficientes de orden inferior que se ajusten usando el cuadrático permanezcan iguales después de volver a colocarlos con un ajuste cúbico. Para ello, tome combinaciones lineales de las columnas de poly(horsepower, 2, raw = TRUE) y haga lo mismo con poly(horsepower, 3, raw = TRUE) modo que las columnas en el ajuste cuadrático sean ortogonales entre sí y de manera similar para el ajuste cúbico. Eso es suficiente para garantizar que los coeficientes de orden inferior no cambien cuando agreguemos coeficientes de orden superior. Observe cómo los primeros tres coeficientes son ahora los mismos en los dos conjuntos a continuación (mientras que arriba difieren). Es decir, en ambos casos, los 3 coeficientes de orden inferior son 23.44592, -120.13774 y 44.08953.

fm2 <- lm(mpg ~ poly(horsepower, 2), Auto)
cbind(coef(fm2))
## [,1]
## (Intercept) 23.44592
## poly(horsepower, 2)1 -120.13774
## poly(horsepower, 2)2 44.08953
fm3 <- lm(mpg ~ poly(horsepower, 3), Auto)
cbind(coef(fm3))
## [,1]
## (Intercept) 23.445918
## poly(horsepower, 3)1 -120.137744
## poly(horsepower, 3)2 44.089528
## poly(horsepower, 3)3 -3.948849

Es importante destacar que, dado que las columnas de poly(horsepwer, 2) son simplemente combinaciones lineales de las columnas de poly(horsepower, 2, raw = TRUE) los dos modelos cuadráticos (ortogonales y brutos) representan los mismos modelos (es decir, dan las mismas predicciones ) y solo difieren en la parametrización. Por ejemplo, los valores ajustados son los mismos:

all.equal(fitted(fm2), fitted(fm2raw))
## [1] TRUE

También podemos verificar que los polinomios tienen columnas ortogonales que también son ortogonales a la intersección:

nr <- nrow(Auto)
e <- rep(1, nr) / sqrt(nr) # constant vector of unit length
p <- cbind(e, poly(Auto$horsepower, 2))
zapsmall(crossprod(p))
## e 1 2
## e 1 0 0
## 1 0 1 0
## 2 0 0 1

Answer 3

Una respuesta rápida es que el poly de un vector es x esencialmente equivalente a la descomposición QR de la matriz cuyas columnas son potencias de x (después del centrado). Por ejemplo:

> x<-rnorm(50)
> x0<-sapply(1:5,function(z) x^z)
> x0<-apply(x0,2,function(z) z-mean(z))
> x0<-qr.Q(qr(x0))
> cor(x0,poly(x,5))
1 2 3 4 5
[1,] -1.000000e+00 -1.113975e-16 -3.666033e-17 7.605615e-17 -1.395624e-17
[2,] -3.812474e-17 1.000000e+00 1.173755e-16 -1.262333e-17 -3.988085e-17
[3,] -7.543077e-17 -7.778452e-17 1.000000e+00 3.104693e-16 -8.472204e-17
[4,] 1.722929e-17 -1.952572e-16 1.013803e-16 -1.000000e+00 -1.611815e-16
[5,] -5.973583e-17 -1.623762e-18 9.163891e-17 -3.037121e-16 1.000000e+00