Si no puedes hacerlo ortogonalmente, hazlo en bruto (regresión polinómica)

Question

Al realizar la regresión polinómica para $Y$ en $X$ A veces se utilizan polinomios brutos, a veces polinomios ortogonales. Pero cuando usan lo que parece completamente arbitrario.

Aquí y aquí se utilizan polinomios brutos. Pero aquí y aquí Los polinomios ortogonales parecen dar los resultados correctos. ¿Qué, cómo, por qué?

Por el contrario, cuando se aprende sobre regresión polinómica en un libro de texto (p. ej. ISLR ), que ni siquiera menciona los polinomios brutos u ortogonales - sólo se da el modelo a ajustar.

Entonces, ¿cuándo tenemos que usar qué?
¿Y por qué los individual valores p para $X$ , $X^2$ etc. ¿difieren mucho entre estos dos valores?

Answer 1

Las variables $X$ y $X^2$ no son linealmente independientes. Por lo tanto, incluso si no hay efecto cuadrático, la adición de $X^2$ al modelo modificará el efecto estimado de $X$ .

Veámoslo con una simulación muy sencilla.

> x <- runif(1e3)
> y <- x + rnorm(length(x))
> summary(lm(y~x))

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.03486    0.06233  -0.559    0.576    
x            1.05843    0.10755   9.841   <2e-16 ***

Ahora con un término cuadrático en el modelo a ajustar.

> summary(lm(y~x+I(x^2)))

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.03275    0.09528   0.344    0.731
x            0.65742    0.44068   1.492    0.136
I(x^2)       0.39914    0.42537   0.938    0.348

Por supuesto, la prueba ómnibus sigue siendo significativa, pero creo que el resultado que buscamos no es éste. La solución es utilizar polinomios ortogonales.

 > summary(lm(y~poly(x,2)))

 Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.49744    0.03098  16.059   <2e-16 ***
poly(x, 2)1  9.63943    0.97954   9.841   <2e-16 ***
poly(x, 2)2  0.91916    0.97954   0.938    0.348    

Obsérvese que los coeficientes de x en el primer modelo y de poly(x,2)1 en el segundo modelo no son iguales, e incluso los interceptos son diferentes. Esto se debe a que poly entrega vectores ortonormales, que también son ortogonales al vector rep(1, length(x)) . Así que poly(x,2)1 no es x sino (x -mean(x))/sqrt(sum((x-mean(x))**2)) ...

Un punto importante es que las pruebas de Wald, en este último modelo, son independientes. Se pueden utilizar polinomios ortogonales para decidir hasta qué grado se quiere llegar, simplemente observando la prueba de Wald: aquí se decide mantener $X$ pero no $X^2$ . Por supuesto, se encontraría el mismo modelo comparando los dos primeros modelos ajustados, pero es más sencillo de este modo; si se considera la posibilidad de subir a grados superiores, es realmente mucho más sencillo.

Una vez que hayas decidido qué términos conservar, puedes volver a los polinomios brutos $X$ y $X^2$ para la interpretabilidad o para la predicción.

Answer 2

Hacer una valoración ingenua de la situación:

En general: supongamos que tenemos dos sistemas diferentes de funciones de base $\{p_n\}_{n=1}^\infty$ así como $\{\tilde{p}\}_{n=1}^\infty$ para algún espacio de funciones (hilbert-), usual $L_2([a,b])$ es decir, el espacio de todas las funciones cuadradas integrables.

Esto significa que cada una de las dos bases puede utilizarse para explicar cada elemento de $L_2([a,b])$ es decir, para $y \in L_2([a,b])$ se tiene para algunos coeficientes $\theta_n$ y $\tilde{\theta}_n \in \mathbb{R}$ , $n=1,2,\dots$ (en el $L_2$ -sentido): $$ \sum_{n=1}^\infty \tilde{\theta}_n \tilde{p}_n = y= \sum_{n=1}^\infty \theta_n p_n.$$

Sin embargo, por otro lado, si se truncan ambos conjuntos de funciones de base en algún número $k<\infty$ es decir, se toma $$\{p_n\}_{n=1}^k$$ así como $$\{\tilde{p}\}_{n=1}^k,$$ es muy probable que estos conjuntos truncados de funciones de base describan "partes diferentes" de $L_2([a,b])$ .

Sin embargo, aquí en el caso especial en el que una base, $\{\tilde{p}\}_{n=1}^\infty$ no es más que una ortogonalización de la otra base, $\{p_n\}_{n=1}^\infty$ la predicción global de $y$ será el mismo para cada modelo truncado ( $\{p\}_{n=1}^k$ y su contraparte ortogonalizada describirán el mismo $k$ -subespacio dimensional de $L_2([a,b])$ ).

Pero cada función de base individual de las dos bases "diferentes" producirá una contribución diferente a esta predicción (¡obviamente como las funciones/predictores son diferentes!) resultando en diferentes $p$ -valores y coeficientes.

Por lo tanto, en términos de predicción no hay (en este caso) ninguna diferencia.

Desde un punto de vista computacional, una matriz modelo formada por funciones de base ortogonales tiene buenas propiedades numéricas/computacionales para el estimador de mínimos cuadrados. Al mismo tiempo, desde el punto de vista estadístico, la ortogonalización da lugar a estimaciones no correlacionadas, ya que $var(\hat{\tilde{\theta}}) = I \sigma²$ bajo los supuestos estándar.

---

La pregunta natural que surge es si existe el mejor sistema de bases truncadas. Sin embargo, la respuesta a la pregunta no es sencilla ni única y depende, por ejemplo, de la definición de la palabra "mejor", es decir, de lo que se intente archivar.