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

36

Estoy tratando de crear un ajuste polinómico de segundo orden para algunos datos que tengo. Digamos que trazo esto en forma con ggplot():

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

Yo obtengo:

trama de ajuste parabólico con banda de confianza en diagrama de dispersión

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

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

Y obtengo:

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, supongo que la fórmula para mi ajuste es:

bar=3.2680.122foo+1.575foo2

Pero eso solo me da los valores incorrectos. Por ejemplo, con siendo 3, esperaría que convierta en algo alrededor de 3.15. Sin embargo, insertando en la fórmula anterior obtengo: foobar

bar=3.2680.1223+1.57532=17.077

¿Lo que da? ¿Estoy interpretando incorrectamente los coeficientes del modelo?

usuario13907
fuente
2
Esta pregunta se responde en varios hilos que se pueden encontrar buscando en nuestro sitio polinomio ortogonal
whuber
66
@whuber Si hubiera sabido que el problema era con los "polinomios ortogonales", probablemente habría encontrado una respuesta. Pero si no sabes qué buscar, es un poco difícil.
user13907
2
También puede encontrar respuestas al buscar en poli , que aparece prominentemente en su código. Pongo dicha información en los comentarios por dos razones: (1) los enlaces pueden ayudar a los futuros lectores, así como a usted mismo y (2) pueden ayudar a mostrarle cómo explotar nuestro sistema de búsqueda (algo idiosincrásico).
whuber
77
¿Ha publicado una pregunta relacionada con su uso de polysin escribir ?polyprimero en R? Eso dice ' Calcular polinomios ortogonales ' en la parte superior en letras grandes y amigables.
Glen_b -Reinstate a Monica
44
@Glen_b Sí, bueno, lo hice tipo de ?polyentender la sintaxis. Es cierto que tengo muy poco conocimiento de los conceptos detrás de esto. No sabía que había algo más (o una diferencia tan grande entre los polinomios "normales" y los polinomios ortogonales), y los ejemplos que vi en línea todos se usaban poly()para el ajuste, especialmente con ggplot, entonces, ¿por qué no usaría eso y confundirse si el resultado fue "incorrecto"? Eso sí, no soy experto en matemáticas, simplemente estoy aplicando lo que he visto hacer a otros y tratando de entenderlo.
user13907

Respuestas:

55

Mi respuesta detallada está a continuación, pero la respuesta general (es decir, real) a este tipo de pregunta es: 1) experimentar, atornillar, mirar los datos, no puedes romper la computadora sin importar lo que hagas, entonces. . . experimentar; o 2) RTFM .

Aquí hay un Rcó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))

El primero lmdevuelve 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 lmdevuelve algo extraño:

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

Dado que lmes lo mismo en las dos llamadas, tiene que ser los argumentos de los lmcuales son diferentes. Entonces, veamos los argumentos. Obviamente, yes lo mismo. Son las otras partes. Veamos las primeras observaciones sobre 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 como se esperaba. La primera columna es xy la segunda columna es x^2. ¿Qué tal la segunda llamada de lm, la 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 muy diferente. La primera columna no es x, y la segunda columna no lo es x^2. Entonces, poly(x,2)haga lo que haga, no regresa xy x^2. Si queremos saber qué polyhace, podríamos comenzar leyendo su archivo de ayuda. Eso decimos nosotros help(poly). La descripción dice:

Devuelve o evalúa polinomios ortogonales de grado 1 a grado sobre el conjunto especificado de puntos x. Todos estos son ortogonales al polinomio constante de grado 0. Alternativamente, evalúe los polinomios en bruto.

Ahora, o sabes qué son los "polinomios ortogonales" o no. Si no lo hace, use Wikipedia o Bing (no Google, por supuesto, porque Google es malvado, no tan malo como Apple, naturalmente, pero sigue siendo malo). O bien, puede decidir que no le importa qué son los polinomios ortogonales. Puede notar la frase "polinomios en bruto" y puede notar un poco más abajo en el archivo de ayuda que polytiene una opción rawque, por defecto, es igual a FALSE. Esas dos consideraciones pueden inspirarlo a probar head(poly(x, 2, raw=TRUE))qué retornos:

            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

Emocionado por este descubrimiento (parece correcto, ahora, ¿sí?), Puede continuar. summary(lm(y ~ poly(x, 2, raw=TRUE))) Esto devuelve:

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

Hay al menos dos niveles para la respuesta anterior. Primero, respondí tu pregunta. Segundo, y mucho más importante, ilustré cómo se supone que debes responder preguntas como esta tú mismo. Cada persona que "sabe programar" ha pasado por una secuencia como la que se menciona más de sesenta millones de veces. Incluso las personas tan deprimentemente malas en la programación como yo pasan por esta secuencia todo el tiempo. Es normal que el código no funcione. Es normal entender mal las funciones que hacen. La forma de lidiar con esto es atornillar, experimentar, mirar los datos y RTFM. Sal del modo de "seguir una receta sin pensar" y entra en el modo "detective".

Cuenta
fuente
77
Creo que esto merece un +6. Intentaré recordar en un par de días cuando eso sea posible. FTR, creo que no tiene por qué ser tan sarcástico, pero hace un buen trabajo al mostrar qué son los polinomios ortogonales / cómo funcionan, y al mostrar el proceso que usas para resolver tales cosas.
gung - Restablece a Monica
13
Gran respuesta, gracias. Aunque estoy un poco ofendido por un "RTFM" (pero tal vez solo soy yo): el problema es que, en todo lo que he leído, al menos con respecto a la regresión lineal en R, la gente a veces hace esto, otros lo hacen. Francamente, no entiendo la entrada de Wikipedia sobre polinomios ortogonales. No se me ocurre por qué uno usaría esto para la regresión si los coeficientes que obtiene son "incorrectos". No soy matemático, trato de seguir las recetas porque no soy un cocinero erudito, pero de todos modos necesito comer algo.
user13907
12
@ user13907, no eres solo tú. De hecho, esta es una buena respuesta que merece ser votada, pero se beneficiaría de tener un tono más agradable.
Waldir Leoncio
8
Realmente no necesita comprender qué polinomios ortogonales hay aquí, solo necesita comprender que no son lo que desea. ¿Por qué alguien podría querer polinomios ortogonales? Envíe cov (poly (x, 2)) para encontrar que la covarianza entre los dos términos en el polinomio es cero (hasta el error de redondeo). Esta es la propiedad clave de los polinomios ortogonales: sus términos tienen cero covarianza entre sí. A veces es conveniente que sus variables de RHS tengan una correlación cero entre sí. Sus coeficientes no son incorrectos, realmente, solo tienen que ser interpretados de manera diferente.
Bill
2
Oh, está bien, esa explicación en inglés simple ahora tiene sentido. Gracias.
user13907
5

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

Y=β0 0+β1X+β2X2+tu

como

Y=metro+β2(F-X)2+tu

metro=β0 0-β12/ /4 4β2β2F=-β1/ /2β2

Durden
fuente
2
+1 Para análisis relacionados, consulte stats.stackexchange.com/questions/28730 y stats.stackexchange.com/questions/157629 .
whuber
4

Si solo desea un empujón en la dirección correcta sin tanto juicio: poly()crea polinomios ortogonales (no correlacionados), en lugar de I(), que ignora por completo la correlación entre los polinomios resultantes. La correlación entre las variables predictoras puede ser un problema en los modelos lineales (consulte aquí para obtener más información sobre por qué la correlación puede ser problemática), por lo que probablemente sea mejor (en general) usar en poly()lugar de I(). Ahora, ¿por qué los resultados se ven tan diferentes? Bueno, ambos poly()y I()tome x y conviértalo en una nueva x (en el caso de I(), la nueva x es solo x ^ 1 o x ^ 2, en el caso de poly(), las nuevas x son mucho más complicadas (si quiere saber de dónde vienen (y probablemente no), puede comenzaraquí o la página de Wikipedia antes mencionada o un libro de texto). El punto es que, cuando estás calculando (prediciendo) y en función de un conjunto particular de valores de x, debes usar los valores de x convertidos producidos por cualquiera poly()o I()(dependiendo de cuál estaba en tu modelo lineal). Asi que:

library(ggplot2)    

set.seed(3)
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))

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

# predict y using modI
coef(modI)[1] + coef(modI)[2] * 3^1 + coef(modI)[3] * 3^2

# predict y using modp
# calculate the new x values using predict.poly()
x_poly <- stats:::predict.poly(object = poly(x,2), newdata = 3)
coef(modp)[1] + coef(modp)[2] * x_poly[1] + coef(modp)[3] * x_poly[2]

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

filups21
fuente
1

'poly' realiza la orto-normalizació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 estaba interesado en la forma funcional. Entonces, ¿cómo expresamos el resultado de 'poli' como expresión? Simplemente invierta el procedimiento de Graham-Schmidt. ¡Terminarás con un desastre!

izmirlig
fuente