¿Cómo interpretar los coeficientes del ajuste de un modelo polinómico?

Question

Estoy intentando crear un ajuste polinómico de segundo orden a unos datos que tengo. Digamos que trazo este ajuste con ggplot() :

ggplot(data, aes(foo, bar)) + geom_point() + 
       geom_smooth(method="lm", formula=y~poly(x, 2))

Lo entiendo:

Por lo tanto, un ajuste de segundo orden funciona bastante bien. Lo calculo con R:

summary(lm(data$bar ~ poly(data$foo, 2)))

Y lo consigo:

lm(formula = data$bar ~ poly(data$foo, 2))
# ...
# Coefficients:
#                     Estimate Std. Error t value Pr(>|t|)    
# (Intercept)         3.268162   0.008282 394.623   <2e-16 ***
# poly(data$foo, 2)1 -0.122391   0.096225  -1.272    0.206
# poly(data$foo, 2)2  1.575391   0.096225  16.372   <2e-16 ***
# ....

Ahora, yo asumiría que la fórmula para mi ajuste es:

$$ \text{bar} = 3.268 - 0.122 \cdot \text{foo} + 1.575 \cdot \text{foo}^2 $$

Pero eso me da los valores equivocados. Por ejemplo, con $\text{foo}$ siendo 3 yo esperaría $\text{bar}$ para convertirse en algo parecido a 3,15. Sin embargo, insertando en la fórmula anterior obtengo:

$$ \text{bar} = 3.268 - 0.122 \cdot 3 + 1.575 \cdot 3^2 = 17.077 $$

¿Qué pasa? ¿Estoy interpretando mal los coeficientes del modelo?

Answer 1

Mi respuesta detallada está más abajo, pero la respuesta general (es decir, real) a este tipo de preguntas es: 1) experimentar, trastear, mirar los datos, no puedes romper el ordenador hagas lo que hagas, así que... experimenta; o 2) lee la documentación.

Aquí hay algunos R código que replica el problema identificado en esta pregunta, más o menos:

# This program written in response to a Cross Validated question
# http://stats.stackexchange.com/questions/95939/
# 
# It is an exploration of why the result from lm(y_x+I(x^2))
# looks so different from the result from lm(y~poly(x,2))

library(ggplot2)

epsilon <- 0.25*rnorm(100)
x       <- seq(from=1, to=5, length.out=100)
y       <- 4 - 0.6*x + 0.1*x^2 + epsilon

# Minimum is at x=3, the expected y value there is
4 - 0.6*3 + 0.1*3^2

ggplot(data=NULL,aes(x, y)) + geom_point() + 
       geom_smooth(method = "lm", formula = y ~ poly(x, 2))

summary(lm(y~x+I(x^2)))       # Looks right
summary(lm(y ~ poly(x, 2)))   # Looks like garbage

# What happened?
# What do x and x^2 look like:
head(cbind(x,x^2))

#What does poly(x,2) look like:
head(poly(x,2))

La primera lm devuelve la respuesta esperada:

Call:
lm(formula = y ~ x + I(x^2))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.53815 -0.13465 -0.01262  0.15369  0.61645 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.92734    0.15376  25.542  < 2e-16 ***
x           -0.53929    0.11221  -4.806 5.62e-06 ***
I(x^2)       0.09029    0.01843   4.900 3.84e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2241 on 97 degrees of freedom
Multiple R-squared:  0.1985,    Adjusted R-squared:  0.182 
F-statistic: 12.01 on 2 and 97 DF,  p-value: 2.181e-05

El segundo lm devuelve algo impar:

Call:
lm(formula = y ~ poly(x, 2))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.53815 -0.13465 -0.01262  0.15369  0.61645 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.24489    0.02241 144.765  < 2e-16 ***
poly(x, 2)1  0.02853    0.22415   0.127    0.899    
poly(x, 2)2  1.09835    0.22415   4.900 3.84e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2241 on 97 degrees of freedom
Multiple R-squared:  0.1985,    Adjusted R-squared:  0.182 
F-statistic: 12.01 on 2 and 97 DF,  p-value: 2.181e-05

Desde lm es el mismo en las dos llamadas, tienen que ser los argumentos de lm que son diferentes. Así pues, veamos los argumentos. Obviamente, y es el mismo. Son las otras partes. Veamos las primeras observaciones de las variables del lado derecho en la primera llamada de lm . El regreso de head(cbind(x,x^2)) parece:

            x         
[1,] 1.000000 1.000000
[2,] 1.040404 1.082441
[3,] 1.080808 1.168146
[4,] 1.121212 1.257117
[5,] 1.161616 1.349352
[6,] 1.202020 1.444853

Esto es lo que se esperaba. La primera columna es x y la segunda columna es x^2 . ¿Qué tal la segunda llamada de lm ¿el que tiene poli? El regreso de head(poly(x,2)) parece:

              1         2
[1,] -0.1714816 0.2169976
[2,] -0.1680173 0.2038462
[3,] -0.1645531 0.1909632
[4,] -0.1610888 0.1783486
[5,] -0.1576245 0.1660025
[6,] -0.1541602 0.1539247

OK, eso es realmente diferente. La primera columna no es x y la segunda columna no es x^2 . Por lo tanto, lo que poly(x,2) no devuelve x y x^2 . Si queremos saber qué poly podemos empezar por leer su archivo de ayuda. Así que decimos help(poly) . La descripción dice:

Devuelve o evalúa polinomios ortogonales de grado 1 a grado sobre el conjunto especificado de puntos x. Estos son todos ortogonales al polinomio constante de grado 0. Alternativamente, evalúa polinomios crudos.

Ahora bien, o sabes lo que son los "polinomios ortogonales" o no lo sabes. Si no lo sabes, entonces usa Wikipedia o Bing (no Google, por supuesto, porque Google es malo - no tan malo como Apple, naturalmente, pero sigue siendo malo). O puede que decidas que no te importa lo que son los polinomios ortogonales. Puede que te fijes en la frase "polinomios brutos" y que notes un poco más abajo en el archivo de ayuda que poly tiene una opción raw que es, por defecto, igual a FALSE . Estas dos consideraciones podrían inspirarle a probar head(poly(x, 2, raw=TRUE)) que devuelve:

            1        2
[1,] 1.000000 1.000000
[2,] 1.040404 1.082441
[3,] 1.080808 1.168146
[4,] 1.121212 1.257117
[5,] 1.161616 1.349352
[6,] 1.202020 1.444853

Entusiasmado por este descubrimiento (parece correcto, ahora, ¿sí?), podría pasar a probar summary(lm(y ~ poly(x, 2, raw=TRUE))) Esto vuelve:

Call:
lm(formula = y ~ poly(x, 2, raw = TRUE))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.53815 -0.13465 -0.01262  0.15369  0.61645 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)              3.92734    0.15376  25.542  < 2e-16 ***
poly(x, 2, raw = TRUE)1 -0.53929    0.11221  -4.806 5.62e-06 ***
poly(x, 2, raw = TRUE)2  0.09029    0.01843   4.900 3.84e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2241 on 97 degrees of freedom
Multiple R-squared:  0.1985,    Adjusted R-squared:  0.182 
F-statistic: 12.01 on 2 and 97 DF,  p-value: 2.181e-05

La respuesta anterior tiene al menos dos niveles. En primer lugar, he respondido a su pregunta. En segundo lugar, y mucho más importante, he ilustrado cómo se supone que debes responder tú mismo a preguntas como ésta. Todas las personas que "saben programar" han pasado por una secuencia como la anterior sesenta millones de veces. Incluso personas tan deprimentemente malas en programación como yo pasan por esta secuencia todo el tiempo. Es normal que el código no funcione. Es normal que no se entienda lo que hacen las funciones. La manera de lidiar con ello es atornillar, experimentar, mirar los datos, y RTFM. Salir del modo "seguir una receta sin sentido" y entrar en el modo "detective".

Answer 2

Si sólo quieres un empujón en la dirección correcta sin tanto juicio: poly() crea polinomios ortogonales (no correlacionados), a diferencia de I() que ignora completamente la correlación entre los polinomios resultantes. La correlación entre las variables predictoras puede ser un problema en los modelos lineales (véase ici para más información sobre por qué la correlación puede ser problemática), por lo que probablemente sea mejor (en general) utilizar poly() en lugar de I() .

Ahora bien, ¿por qué los resultados son tan diferentes? Bueno, ambos poly() y I() tomar x y convertirla en una nueva x. En el caso de I() la nueva x es simplemente x^1 o x^2. En el caso de poly() Los nuevos X son mucho más complicados. Si quieres saber de dónde vienen (y probablemente no lo sepas), puedes empezar ici o el mencionado Página de Wikipedia o un libro de texto.

La cuestión es que, cuando se está calculando (prediciendo) y en base a un conjunto particular de valores x, es necesario utilizar los valores x convertidos producidos por poly() o I() (dependiendo de cuál fuera su modelo lineal). Así que:

library(ggplot2)    

# set the seed to make the results reproducible.
set.seed(3)

#### simulate some data ####
# epsilon = random error term
epsilon <- 0.25*rnorm(100)
# x values are just a sequence from 1 to 5
x       <- seq(from=1, to=5, length.out=100)
# y is a polynomial function of x (plus some error)
y       <- 4 - 0.6*x + 0.1*x^2 + epsilon

# Minimum is at x=3, the expected y value there is
4 - 0.6*3 + 0.1*3^2

# visualize the data (with a polynomial best-fit line)
ggplot(data=NULL,aes(x, y)) + geom_point() + 
   geom_smooth(method = "lm", formula = y ~ poly(x, 2))

#### Model the data ####
# first we'll try to model the data with just I()
modI <- lm(y~x+I(x^2)) 
# the model summary looks right
summary(modI)

# next we'll try it with poly()
modp <- lm(y ~ poly(x, 2))
# the model summary looks weird
summary(modp)

#### make predictions at x = 3 based on each model ####
# predict y using modI
# all we need to do is take the model coefficients and plug them into the formula: intercept + beta1 * x^1 + beta2 * x^2
coef(modI)[1] + coef(modI)[2] * 3^1 + coef(modI)[3] * 3^2

(Intercepción)
3.122988

# predict y using modp
# this takes an extra step.
# first, calculate the new x values using predict.poly()
x_poly <- stats:::predict.poly(object = poly(x,2), newdata = 3)
# then use the same formula as above, but this time instead of the raw x value (3), use the polynomial x-value (x_poly)
coef(modp)[1] + coef(modp)[2] * x_poly[1] + coef(modp)[3] * x_poly[2]

(Intercepción)
3.122988

En este caso, ambos modelos devuelven la misma respuesta, lo que sugiere que la correlación entre las variables predictoras no está influyendo en sus resultados. Si la correlación fuera un problema, los dos métodos predecirían valores diferentes.

Answer 3

Hay un enfoque interesante para la interpretación de la regresión polinómica por Stimson et al. (1978) . Se trata de reescribir

$Y = \beta_{0} + \beta_{1} X + \beta_{2} X^{2} + u$

como

$Y = m + \beta_{2} \left( f - X \right)^{2} + u$

donde $m = \beta_{0} - \left. \beta_{1}^{2} \right/ 4 \beta_{2}$ es el mínimo o el máximo (según el signo de $\beta_{2}$ ) y $f = \left. -\beta_{1} \right/ 2 \beta_{2}$ es el valor focal. Básicamente, transforma la combinación tridimensional de pendientes en una parábola en dos dimensiones. Su documento ofrece un ejemplo de ciencias políticas.

Answer 4

'poly' realiza la ortonormalización de Graham-Schmidt en los polinomios 1, x, x^2, ..., x^deg Por ejemplo, esta función hace lo mismo que 'poly' sin devolver los atributos 'coef', por supuesto.

MyPoly <- 
function(x, deg)
{
    n <- length(x)
    ans <- NULL
    for(k in 1:deg)
    {
        v <- x^k
        cmps <- rep(0, n)
        if(k>0) for(j in 0:(k-1)) cmps <- cmps + c(v%*%ans[,j+1])*ans[,j+1]
        p <- v - cmps
        p <- p/sum(p^2)^0.5
        ans <- cbind(ans, p)
    }
    ans[,-1]
}

Llegué a este hilo porque me interesaba la forma funcional. Entonces, ¿cómo expresamos el resultado de 'poly' como una expresión? Simplemente invierte el procedimiento Graham-Schmidt. ¡Acabarás con un lío!