Ajuste del modelo polinomial a los datos en R

He leído las respuestas a esta pregunta y son bastante útiles, pero necesito ayuda especialmente en R.

Tengo un conjunto de datos de ejemplo en R de la siguiente manera:

x <- c(32,64,96,118,126,144,152.5,158)  
y <- c(99.5,104.8,108.5,100,86,64,35.3,15)

Quiero ajustar un modelo a estos datos para que y = f(x). Quiero que sea un modelo polinomial de tercer orden.

¿Cómo puedo hacer eso en R?

Además, ¿puede R ayudarme a encontrar el modelo que mejor se adapte?

Para obtener un polinomio de tercer orden en x (x ^ 3), puede hacer

lm(y ~ x + I(x^2) + I(x^3))

o

lm(y ~ poly(x, 3, raw=TRUE))

Podría ajustar un polinomio de décimo orden y obtener un ajuste casi perfecto, pero ¿debería hacerlo?

EDITAR: poly (x, 3) es probablemente una mejor opción (ver @hadley a continuación).

Qué modelo es el "modelo que mejor se ajusta" depende de lo que usted quiera decir con "mejor". R tiene herramientas para ayudar, pero necesita proporcionar la definición de "mejor" para elegir entre ellas. Considere los siguientes datos y código de ejemplo:

x <- 1:10
y <- x + c(-0.5,0.5)

plot(x,y, xlim=c(0,11), ylim=c(-1,12))

fit1 <- lm( y~offset(x) -1 )
fit2 <- lm( y~x )
fit3 <- lm( y~poly(x,3) )
fit4 <- lm( y~poly(x,9) )
library(splines)
fit5 <- lm( y~ns(x, 3) )
fit6 <- lm( y~ns(x, 9) )

fit7 <- lm( y ~ x + cos(x*pi) )

xx <- seq(0,11, length.out=250)
lines(xx, predict(fit1, data.frame(x=xx)), col='blue')
lines(xx, predict(fit2, data.frame(x=xx)), col='green')
lines(xx, predict(fit3, data.frame(x=xx)), col='red')
lines(xx, predict(fit4, data.frame(x=xx)), col='purple')
lines(xx, predict(fit5, data.frame(x=xx)), col='orange')
lines(xx, predict(fit6, data.frame(x=xx)), col='grey')
lines(xx, predict(fit7, data.frame(x=xx)), col='black')

¿Cuál de esos modelos es el mejor? Se podrían hacer argumentos para cualquiera de ellos (pero yo no querría usar el púrpura para la interpolación).

Con respecto a la pregunta '¿puede R ayudarme a encontrar el modelo que mejor se ajusta?', Probablemente haya una función para hacer esto, asumiendo que puede indicar el conjunto de modelos para probar, pero este sería un buen primer enfoque para el conjunto de n-1 polinomios de grado:

polyfit <- function(i) x <- AIC(lm(y~poly(x,i)))
as.integer(optimize(polyfit,interval = c(1,length(x)-1))$minimum)

Notas

• La validez de este enfoque dependerá de sus objetivos, los supuestos de optimize()y AIC()y si AIC es el criterio que desea utilizar,

• polyfit()puede que no tenga un mínimo único. verifique esto con algo como:

  for (i in 2:length(x)-1) print(polyfit(i))
  

• Usé la as.integer()función porque no me queda claro cómo interpretaría un polinomio no entero.

• para probar un conjunto arbitrario de ecuaciones matemáticas, considere el programa 'Eureqa' revisado por Andrew Gelman aquí

Actualizar

Consulte también la stepAICfunción (en el paquete MASS) para automatizar la selección del modelo.

La forma más fácil de encontrar el mejor ajuste en R es codificar el modelo como:

lm.1 <- lm(y ~ x + I(x^2) + I(x^3) + I(x^4) + ...)

Después de usar la regresión AIC descendente

lm.s <- step(lm.1)