¿Regresión polinómica en bruto u ortogonal?

Quiero hacer una regresión de una variable en . ¿Debo hacer esto usando polinomios en bruto u ortogonales? Miré la pregunta en el sitio que se ocupa de estos, pero realmente no entiendo cuál es la diferencia entre usarlos. $y$$x,x^2,\ldots,x^5$

¿Por qué no puedo simplemente hacer una regresión "normal" para obtener los coeficientes de$\beta_i$$y=\sum_{i=0}^5 \beta_i x^i$ (junto con los valores de p y todas las otras cosas buenas) y en su lugar ¿tiene que preocuparse si usa polinomios en bruto u ortogonales? Esta elección me parece estar fuera del alcance de lo que quiero hacer.

En el libro de estadísticas que estoy leyendo actualmente (ISLR por Tibshirani et al) estas cosas no fueron mencionadas. En realidad, fueron minimizados de alguna manera.
La razón es, AFAIK, que en la lm()función en R, usar y ~ poly(x, 2)cantidades para usar polinomios ortogonales y usar y ~ x + I(x^2)cantidades para usar materias primas. Pero en las págs. 116 los autores dicen que usamos la primera opción porque la segunda es "engorrosa", lo que no deja ninguna indicación de que estos comandos realmente hagan cosas completamente diferentes (y que tengan diferentes resultados como consecuencia).
(tercera pregunta) ¿Por qué los autores de ISLR confundirían así a sus lectores?

Creo que la respuesta es menos sobre la estabilidad numérica (aunque eso juega un papel) y más sobre la reducción de la correlación.

En esencia, el problema se reduce al hecho de que cuando retrocedemos contra un grupo de polinomios de alto orden, las covariables contra las que retrocedemos se correlacionan altamente. Código de ejemplo a continuación:

x = rnorm(1000)
raw.poly = poly(x,6,raw=T)
orthogonal.poly = poly(x,6)
cor(raw.poly)
cor(orthogonal.poly)

Esto es tremendamente importante. A medida que las covariables se vuelven más correlacionadas, nuestra capacidad para determinar cuáles son importantes (y cuál es el tamaño de sus efectos) se erosiona rápidamente. Esto generalmente se conoce como el problema de la multicolinealidad. En el límite, si tuviéramos dos variables que estuvieran completamente correlacionadas, cuando las retrocedemos contra algo, es imposible distinguir entre las dos; puede pensar en esto como una versión extrema del problema, pero este problema afecta nuestras estimaciones para grados menores de correlación también. Por lo tanto, en un sentido real, incluso si la inestabilidad numérica no fuera un problema, la correlación de los polinomios de orden superior hace un daño tremendo a nuestras rutinas de inferencia. Esto se manifestará como errores estándar más grandes (y, por lo tanto, estadísticas t más pequeñas) que de otro modo vería (ver ejemplo de regresión a continuación).

y = x*2 + 5*x**3 - 3*x**2 + rnorm(1000)
raw.mod = lm(y~poly(x,6,raw=T))
orthogonal.mod = lm(y~poly(x,6))
summary(raw.mod)
summary(orthogonal.mod)

Si ejecuta este código, la interpretación es un poco difícil porque todos los coeficientes cambian y, por lo tanto, las cosas son difíciles de comparar. Sin embargo, al observar las estadísticas T, podemos ver que la capacidad de determinar los coeficientes era MUCHO mayor con los polinomios ortogonales. Para los 3 coeficientes relevantes, obtuve estadísticas t de (560,21,449) para el modelo ortogonal, y solo (28, -38,121) para el modelo polinomial bruto. Esta es una gran diferencia para un modelo simple con solo unos pocos términos polinómicos de orden relativamente bajo que importaban.

Eso no quiere decir que esto viene sin costos. Hay dos costos principales a tener en cuenta. 1) perdemos cierta interpretabilidad con polinomios ortogonales. Podríamos entender qué x**3significa el coeficiente en , pero interpretar el coeficiente en x**3-3x(el tercer poli hermita, no necesariamente lo que usará) puede ser mucho más difícil. Segundo, cuando decimos que estos son polinomios son ortogonales, queremos decir que son ortogonales con respecto a alguna medida de distancia. Elegir una medida de distancia que sea relevante para su situación puede ser difícil. Sin embargo, una vez dicho esto, creo que la polyfunción está diseñada para elegir de modo que sea ortogonal con respecto a la covarianza, lo cual es útil para regresiones lineales.

¿Por qué no puedo simplemente hacer una regresión "normal" para obtener los coeficientes?

$0.4$$0.4000000059604644775390625$

El uso de polinomios en bruto causará problemas porque tendremos un gran número. Aquí hay una pequeña prueba: estamos comparando el número de condición de la matriz con el polinomio crudo y ortogonal.

> kappa(model.matrix(mpg~poly(wt,10),mtcars))
[1] 5.575962
> kappa(model.matrix(mpg~poly(wt,10, raw = T),mtcars))
[1] 2.119183e+13

También puedes consultar mi respuesta aquí para ver un ejemplo.

¿Por qué hay coeficientes grandes para polinomios de orden superior?

Siento que varias de estas respuestas pierden completamente el punto. La respuesta de Haitao aborda los problemas computacionales con el ajuste de polinomios en bruto, pero está claro que OP pregunta por las diferencias estadísticas entre los dos enfoques. Es decir, si tuviéramos una computadora perfecta que pudiera representar todos los valores exactamente, ¿por qué preferiríamos un enfoque sobre el otro?

$R^2$$X$$Y$$X=0$$X=0$$X$

data("iris")

#Raw:
fit.raw <- lm(Petal.Length ~ Petal.Width + I(Petal.Width^2) +
                  I(Petal.Width^3), data = iris)
summary(fit.raw)

#> Coefficients:
#>                  Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)        1.1034     0.1304   8.464 2.50e-14 ***
#> Petal.Width        1.1527     0.5836   1.975  0.05013 .  
#> I(Petal.Width^2)   1.7100     0.5487   3.116  0.00221 ** 
#> I(Petal.Width^3)  -0.5788     0.1408  -4.110 6.57e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.3898 on 146 degrees of freedom
#> Multiple R-squared:  0.9522, Adjusted R-squared:  0.9512 
#> F-statistic: 969.9 on 3 and 146 DF,  p-value: < 2.2e-16

#Orthogonal
fit.orth <- lm(Petal.Length ~ stats::poly(Petal.Width, 3), data = iris)

#Marginal effect of X at X=0 from orthogonal model
library(margins)
summary(margins(fit.orth, variables = "Petal.Width", 
                at = data.frame(Petal.Width = 0)))
#> Warning in check_values(data, at): A 'at' value for 'Petal.Width' is
#> outside observed data range (0.1,2.5)!
#>       factor Petal.Width    AME     SE      z      p  lower  upper
#>  Petal.Width      0.0000 1.1527 0.5836 1.9752 0.0482 0.0089 2.2965

^{Creado el 25/10/2019 por el paquete reprex (v0.3.0)}

El efecto marginal de Petal.Widtha 0 del ajuste ortogonal y su error estándar son exactamente iguales a los del ajuste polinómico bruto. El uso de polinomios ortogonales no mejora la precisión de las estimaciones de la misma cantidad entre los dos modelos.

$Y$$X$$Y$$X$

library(jtools)
data("iris")

fit.raw3 <- lm(Petal.Length ~ Petal.Width + I(Petal.Width^2) +
                  I(Petal.Width^3), data = iris)
fit.raw1 <- lm(Petal.Length ~ Petal.Width, data = iris)

round(summ(fit.raw3, part.corr = T)$coef, 3)
#>                    Est.  S.E. t val.     p partial.r part.r
#> (Intercept)       1.103 0.130  8.464 0.000        NA     NA
#> Petal.Width       1.153 0.584  1.975 0.050     0.161  0.036
#> I(Petal.Width^2)  1.710 0.549  3.116 0.002     0.250  0.056
#> I(Petal.Width^3) -0.579 0.141 -4.110 0.000    -0.322 -0.074

round(summ(fit.raw1, part.corr = T)$coef, 3)
#>              Est.  S.E. t val. p partial.r part.r
#> (Intercept) 1.084 0.073 14.850 0        NA     NA
#> Petal.Width 2.230 0.051 43.387 0     0.963  0.963

fit.orth3 <- lm(Petal.Length ~ stats::poly(Petal.Width, 3), 
               data = iris)
fit.orth1 <- lm(Petal.Length ~ stats::poly(Petal.Width, 3)[,1], 
               data = iris)

round(summ(fit.orth3, part.corr = T)$coef, 3)
#>                                Est.  S.E.  t val. p partial.r part.r
#> (Intercept)                   3.758 0.032 118.071 0        NA     NA
#> stats::poly(Petal.Width, 3)1 20.748 0.390  53.225 0     0.975  0.963
#> stats::poly(Petal.Width, 3)2 -3.015 0.390  -7.735 0    -0.539 -0.140
#> stats::poly(Petal.Width, 3)3 -1.602 0.390  -4.110 0    -0.322 -0.074

round(summ(fit.orth1, part.corr = T)$coef, 3)
#>                                    Est.  S.E. t val. p partial.r part.r
#> (Intercept)                       3.758 0.039 96.247 0        NA     NA
#> stats::poly(Petal.Width, 3)[, 1] 20.748 0.478 43.387 0     0.963  0.963

^{Creado el 25/10/2019 por el paquete reprex (v0.3.0)}

$0.001$$0.003$$0.005$$0.927$$0.927$$0.020$$0.005$$0.927$. Por el modelo polinomial ortogonal pero no por el modelo polinomial bruto, sabemos que la mayor parte de la varianza explicada en el resultado se debe al término lineal, con muy poco proveniente del término cuadrado y aún menos del término cúbico. Los valores polinómicos brutos no cuentan esa historia.

Ahora, si desea este beneficio interpretativo sobre el beneficio interpetacional de ser capaz de comprender los coeficientes del modelo, entonces debe usar polinomios ortogonales. Si prefiere mirar los coeficientes y saber exactamente lo que significan (aunque dudo que uno lo haga), entonces debería usar los polinomios en bruto. Si no le importa (es decir, solo desea controlar la confusión o generar valores pronosticados), entonces realmente no importa; ambas formas llevan la misma información con respecto a esos objetivos. También diría que los polinomios ortogonales deberían preferirse en la regularización (p. Ej., Lazo), porque eliminar los términos de orden superior no afecta los coeficientes de los términos de orden inferior, lo cual no es cierto con los polinomios en bruto,

Corroboro la excelente respuesta de @ user5957401 y agrego comentarios sobre interpolación, extrapolación e informes.

Incluso en el dominio de valores de parámetros estables, los coeficientes / parámetros modelados por los polinomios ortogonales tendrán errores estándar sustancialmente más pequeños que los coeficientes / parámetros modelados por los parámetros brutos. Esencialmente, los polinomios ortogonales son un conjunto libre de descriptores de covarianza cero. ¡Eso es PCA gratis!

El único inconveniente potencial es tener que explicar esto a alguien que no comprende la virtud de los descriptores de covarianza cero. Los coeficientes no son inmediatamente interpretables en el contexto de los efectos de primer orden (tipo velocidad) o de segundo orden (tipo aceleración). Esto puede ser bastante condenatorio en un entorno empresarial.

$10^{-d}$$\sim$$R^2$$adj ~R^2$

Por lo tanto, sería "órdenes de magnitud" más seguro al informar sobre el modelo ortogonal que el modelo en bruto. En la práctica, interpolaría con cualquiera de los modelos, pero extrapolaría solo con el ortogonal.

Hubiera comentado para mencionar esto, pero no tengo suficiente representante, así que intentaré expandirme en una respuesta. Tal vez le interese ver que en la Sección 7.8.1 del Laboratorio en "Introducción al aprendizaje estadístico" (James et. Al., 2017, 8a impresión corregida), discuten algunas diferencias entre usar polinomios ortogonales o no, que es usar el raw=TRUEo raw=FALSEen la poly()función. Por ejemplo, los coeficientes estimados cambiarán, pero los valores ajustados no:

# using the Wage dataset in the ISLR library
fit1 <- lm(wage ~ poly(age, 4, raw=FALSE), data=Wage)
fit2 <- lm(wage ~ poly(age, 4, raw=TRUE), data=Wage)
print(coef(fit1)) # coefficient estimates differ
print(coef(fit2))
all.equal(predict(fit1), predict(fit2)) #returns TRUE    

El libro también discute cómo cuando se usan polinomios ortogonales, los valores p obtenidos usando la anova()prueba F anidada (para explorar qué grado de polinomio podría estar justificado) son los mismos que los obtenidos cuando se usa la prueba t estándar, emitidos por summary(fit). Esto ilustra que el estadístico F es igual al cuadrado del estadístico t en ciertas situaciones.