¿Regresión polinómica bruta u ortogonal?

Question

Quiero hacer una regresión de una variable $y$ en $x,x^2,\ldots,x^5$ . ¿Debo hacerlo con polinomios brutos u ortogonales? Miré la pregunta en el sitio que trata de estos, pero no entiendo realmente cuál es la diferencia entre el uso de ellos.

¿Por qué no puedo hacer una regresión "normal" para obtener los coeficientes $\beta_i$ de $y=\sum_{i=0}^5 \beta_i x^i$ (junto con los valores p y todas las demás cosas bonitas) y en su lugar tener que preocuparse de si se utilizan polinomios brutos u ortogonales? Esta elección me parece que está fuera del alcance de lo que quiero hacer.

En el libro de estadísticas que estoy leyendo actualmente (ISLR de Tibshirani et al) no se mencionan estas cosas. En realidad, se les restó importancia en cierto modo.
La razón es que, según la información disponible, en el lm() en R, utilizando y ~ poly(x, 2) equivale a utilizar polinomios ortogonales y usar y ~ x + I(x^2) equivale a utilizar las crudas. Pero en la página 116 los autores dicen que usamos la primera opción porque la segunda es "engorrosa", lo que no indica que estos comandos realmente hacen dos cosas completamente diferentes (y tienen resultados diferentes como consecuencia).
(tercera pregunta) ¿Por qué los autores de ISLR confunden así a sus lectores?

Answer 1

Me parece que varias de estas respuestas no entienden nada. La respuesta de Haitao aborda la computacional problemas con el ajuste de polinomios crudos, pero está claro que el OP está preguntando por el estadística diferencias entre los dos enfoques. Es decir, si tuviéramos un ordenador perfecto que pudiera representar todos los valores con exactitud, ¿por qué preferiríamos un enfoque sobre el otro?

El usuario5957401 argumenta que los polinomios ortogonales reducen la colinealidad entre las funciones polinómicas, lo que hace que su estimación sea más estable. Estoy de acuerdo con la crítica de Jake Westfall; los coeficientes de los polinomios ortogonales representan cantidades completamente diferentes de los coeficientes de los polinomios brutos. La función dosis-respuesta implícita en el modelo, $R^2$ El valor de la predicción, el MSE, los valores predichos y los errores estándar de los valores predichos serán idénticos independientemente de si se utilizan polinomios ortogonales o crudos. El coeficiente de $X$ en una regresión polinómica bruta de orden 2 tiene la interpretación de "el cambio instantáneo en $Y$ cuando $X=0$ ." Si se realizó un procedimiento de efectos marginales en el polinomio ortogonal donde $X=0$ se obtendría exactamente la misma pendiente y el mismo error estándar, aunque el coeficiente y el error estándar del término de primer orden en la regresión polinómica ortogonal sean completamente diferentes de su valor en la regresión polinómica bruta. Es decir, cuando se trata de obtener las mismas cantidades de ambas regresiones (es decir, cantidades que pueden interpretarse de la misma manera), las estimaciones y los errores estándar serán idénticos. Utilizar polinomios ortogonales no significa que mágicamente se tenga más certeza de la pendiente de $X$ en un punto determinado. La estabilidad de los modelos es idéntica. Véase a continuación:

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 2019-10-25 por el paquete reprex (v0.3.0)

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

La clave es la siguiente: el uso de polinomios ortogonales permite aislar la contribución de cada término a la explicación de la varianza en el resultado, por ejemplo, medida por la correlación semiparcial al cuadrado. Si se ajusta un polinomio ortogonal de orden 3, la correlación semiparcial al cuadrado de cada término representa la varianza del resultado explicada por ese término en el modelo. Por lo tanto, si quiere responder a la pregunta "¿Qué parte de la varianza en $Y$ se explica la componente lineal de $X$ ?" se podría ajustar una regresión polinómica ortogonal, y la correlación semiparcial al cuadrado en el término lineal representaría esta cantidad. Esto no es así con los polinomios brutos. Si se ajusta un modelo polinómico bruto del mismo orden, la correlación parcial al cuadrado en el término lineal no representa la proporción de la varianza en $Y$ explicada por el componente lineal de $X$ . Véase más abajo.

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 2019-10-25 por el paquete reprex (v0.3.0)

Las correlaciones semiparciales al cuadrado para los polinomios brutos cuando se ajusta el polinomio de orden 3 son $0.001$ , $0.003$ y $0.005$ . Cuando sólo se ajusta el término lineal, la correlación semiparcial al cuadrado es $0.927$ . Las correlaciones semiparciales al cuadrado para los polinomios ortogonales cuando se ajusta el polinomio de orden 3 son $0.927$ , $0.020$ y $0.005$ . Cuando sólo se ajusta el término lineal, la correlación semiparcial al cuadrado sigue siendo $0.927$ . A partir del modelo polinómico ortogonal, pero no del modelo polinómico en bruto, sabemos que la mayor parte de la varianza explicada en el resultado se debe al término lineal, con muy poco procedente del término cuadrado y aún menos del término cúbico. Los valores polinómicos brutos no cuentan esa historia.

Ahora bien, si quiere este beneficio interpretativo sobre el beneficio interpretativo de poder entender realmente los coeficientes del modelo, entonces debe utilizar polinomios ortogonales. Si prefiere mirar los coeficientes y saber exactamente lo que significan (aunque dudo que uno lo haga típicamente), entonces debería usar los polinomios brutos. Si no le importa (es decir, sólo quiere controlar la confusión o generar valores predichos), 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 ser preferidos en la regularización (por ejemplo, el lazo), porque la eliminación de los términos de orden superior no afecta a los coeficientes de los términos de orden inferior, lo que no es cierto con los polinomios crudos, y las técnicas de regularización a menudo se preocupan por el tamaño de cada coeficiente.

Answer 2

Creo que la respuesta no tiene tanto que ver con la estabilidad numérica (aunque ésta juega un papel importante) como con la reducción de la correlación.

En esencia, el problema se reduce al hecho de que cuando hacemos una regresión contra un grupo de polinomios de alto orden, las covariables contra las que hacemos la regresión se vuelven altamente correlacionadas. 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 están más correlacionadas, nuestra capacidad para determinar cuáles son importantes (y cuál es la magnitud de sus efectos) se erosiona rápidamente. Esto se conoce como el problema de la multicolinealidad. En el límite, si tuviéramos dos variables totalmente correlacionadas, cuando las regresamos contra algo, es imposible distinguir entre las dos -- se puede pensar en esto como una versión extrema del problema, pero este problema afecta a 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 tanto, estadísticas t más pequeñas) que de otro modo se verían (véase el ejemplo de regresión más abajo). Por esta razón, podríamos optar por ortogonalizar nuestros polinomios antes de hacer la regresió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 se 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, si miramos 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 sólo (28,-38,121) para el modelo polinómico crudo. Esta es una gran diferencia para un modelo simple con sólo unos pocos términos polinómicos de orden relativamente bajo que importaban.

Eso no quiere decir que esto no tenga costes. Hay que tener en cuenta dos costes principales 1) Perdemos algo de interpretabilidad con los polinomios ortogonales. Podemos entender lo que el coeficiente de x**3 significa, pero la interpretación del coeficiente de x**3-3x (el tercer polígono hermético - no necesariamente el que se utilizará) puede ser mucho más difícil. En segundo lugar, cuando decimos que estos polinomios son ortogonales, queremos decir que son ortogonales con respecto a alguna medida de distancia. Escoger una medida de distancia que sea relevante para tu situación puede ser difícil. Sin embargo, una vez dicho esto, creo que el poly está diseñada para elegir de manera que sea ortogonal con respecto a la covarianza, lo cual es útil para las regresiones lineales.

Answer 3

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

Porque no es numéricamente estable. Recuerda que el ordenador utiliza un número fijo de bits para representar un número flotante. Compruebe IEEE754 para los detalles, le sorprenderá que incluso un simple número $0.4$ El ordenador tiene que almacenarlo como $0.4000000059604644775390625$ . Puedes probar con otros números aquí

Usar el polinomio crudo causará problemas porque tendremos un número enorme. Aquí hay una pequeña prueba: estamos comparando número de condición de la matriz con polinomio bruto 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 grandes coeficientes para el polinomio de orden superior

Answer 4

Sólo habría comentado para mencionar esto, pero no tengo suficiente rep, así que intentaré ampliar en una respuesta. Quizá te interese ver que en la sección 7.8.1 del laboratorio de "Introducción al aprendizaje estadístico" (James et. al., 2017, octava edición corregida), sí que se comentan algunas diferencias entre usar polinomios ortogonales o no, que es usar el raw=TRUE o raw=FALSE en el poly() función. Por ejemplo, las estimaciones de los coeficientes 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 analiza cómo cuando se utilizan polinomios ortogonales, los valores p obtenidos mediante el anova() La prueba F anidada (para explorar qué grado de polinomio podría estar justificado) son las mismas que las obtenidas al utilizar la prueba t estándar, producida por summary(fit) . Esto ilustra que el estadístico F es igual al cuadrado del estadístico t en determinadas situaciones.