¿Por qué obtengo diferentes predicciones para la expansión polinómica manual y uso la función R `poly`?

¿Por qué obtengo diferentes predicciones para la expansión polinómica manual y uso de la polyfunción R ?

set.seed(0)
x <- rnorm(10)
y <- runif(10)
plot(x,y,ylim=c(-0.5,1.5))
grid()

# xp is a grid variable for ploting
xp <- seq(-3,3,by=0.01)
x_exp <- data.frame(f1=x,f2=x^2)
fit <- lm(y~.-1,data=x_exp)
xp_exp <- data.frame(f1=xp,f2=xp^2)
yp <- predict(fit,xp_exp)
lines(xp,yp)

# using poly function
fit2 <- lm(y~ poly(x,degree=2) -1)
yp <- predict(fit2,data.frame(x=xp))
lines(xp,yp,col=2)

Mi intento:

• Parece ser un problema con la intercepción, cuando ajusto el modelo con intercepción, es decir, no -1en el modelo formula, las dos líneas son las mismas. Pero ¿por qué sin la intersección las dos líneas son diferentes?

• Otra "solución" es usar rawexpansión polinómica en lugar de polinomio ortogonal. Si cambiamos el código en fit2 = lm(y~ poly(x,degree=2, raw=T) -1), haremos 2 líneas iguales. ¿Pero por qué?

Como notó correctamente, la diferencia original se debe a que en el primer caso usa los polinomios "en bruto" mientras que en el segundo caso usa los polinomios ortogonales. Por lo tanto, si la lmllamada posterior se modificara en: fit3<-lm(y~ poly(x,degree=2, raw = TRUE) -1)obtendríamos los mismos resultados entre fity fit3. La razón por la que obtenemos los mismos resultados en este caso es "trivial"; Nos ajustamos exactamente al mismo modelo que teníamos fit<-lm(y~.-1,data=x_exp), sin sorpresas.

Uno puede verificar fácilmente que las matrices de modelos de los dos modelos sean las mismas all.equal( model.matrix(fit), model.matrix(fit3) , check.attributes= FALSE) # TRUE).

Lo que es más interesante es por qué obtendrá los mismos gráficos al usar una intersección. Lo primero que debe notar es que, al ajustar un modelo con una intersección

• En el caso de fit2, simplemente movemos las predicciones del modelo verticalmente; La forma real de la curva es la misma.

• Por otro lado, incluir una intersección en el caso de los fitresultados en no solo una línea diferente en términos de ubicación vertical sino con una forma totalmente diferente en general.

Podemos ver eso fácilmente simplemente agregando los siguientes ajustes en la trama existente.

fit_b<-lm(y~. ,data=x_exp)
yp=predict(fit_b,xp_exp)
lines(xp,yp, col='green', lwd = 2)

fit2_b<-lm(y~ poly(x,degree=2, raw = FALSE) )
yp=predict(fit2_b,data.frame(x=xp))
lines(xp,yp,col='blue')

OK ... ¿Por qué los ajustes sin intersección son diferentes mientras que los ajustes con intersección son los mismos? La captura vuelve a estar en la condición de ortogonalidad.

En el caso de que fit_bla matriz modelo utilizada contenga elementos no ortogonales, la matriz de Gram crossprod( model.matrix(fit_b) )está lejos de ser diagonal; en el caso de fit2_blos elementos son ortogonales ( crossprod( model.matrix(fit2_b) )es efectivamente diagonal).

Como tal, en el caso de fitcuando lo expandimos para incluir una intersección fit_b, cambiamos las entradas fuera de la diagonal de la matriz de Gram y, por lo tanto, el ajuste resultante es diferente en su conjunto (curvatura, intersección, etc.) diferentes en comparación con el ajuste proporcionado por . Sin embargo, en el caso de que cuando lo expandimos para incluir una intersección ya que solo agregamos una columna que ya es ortogonal a las columnas que teníamos, la ortogonalidad está en contra del polinomio constante de grado 0 . Esto simplemente resulta en mover verticalmente nuestra línea ajustada por la intersección. Por eso las tramas son diferentes.$X^TX$fitfit2fit2_b

La interesante pregunta secundaria es por qué fit_by fit2_bson iguales; después de todo, las matrices modelo de fit_by fit2_bno son iguales en valor nominal . Aquí solo necesitamos recordar eso en última instancia fit_by fit2_btener la misma información. fit2_bes solo una combinación lineal de fit_blo que esencialmente sus ajustes resultantes serán los mismos. Las diferencias observadas en el coeficiente ajustado reflejan la recombinación lineal de los valores de fit_bpara obtenerlos ortogonales. (Vea la respuesta de G. Grothendieck aquí también para un ejemplo diferente).