2010-09-29 12 views
67

he leído las respuestas a esta question y son bastante útiles, pero necesito ayuda sobre todo en R.modelo polinomial de montaje a los datos en R

tengo un ejemplo de datos establecidos en I 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 ajusta?

Respuesta

71

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

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

o

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

Se podría colocar una orden 10 de polinomio y obtener un ajuste casi perfecto , pero deberias?

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

+6

en el clavo en pedir "en caso de que". Los datos de muestra solo tienen 8 puntos. Los grados de libertad son bastante bajos aquí. Los datos de la vida real pueden tener mucho más, por supuesto. –

+1

Gracias por su respuesta. ¿Qué hay de hacer que R encuentre el modelo que mejor se adapte? ¿Hay alguna función para esto? –

+4

Depende de su definición de "mejor modelo". El modelo que le proporcione el mayor R^2 (que sería un polinomio de décimo orden) no es necesariamente el "mejor" modelo. Los términos en su modelo deben ser razonablemente elegidos. Puede obtener un ajuste casi perfecto con una gran cantidad de parámetros, pero el modelo no tendrá poder predictivo y será inútil para otra cosa que no sea dibujar una línea que mejor se ajuste a los puntos. – Greg

12

En cuanto a la pregunta "¿puede ayudarme R a encontrar el mejor modelo de ajuste?", Probablemente exista una función para hacerlo, suponiendo que pueda establecer el conjunto de modelos a probar, pero este sería un buen primer acercamiento 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 no tener un mínimo. comprobar esto con algo como:

    for (i in 2:length(x)-1) print(polyfit(i)) 
    
  • que utiliza la función as.integer() porque no me queda claro cómo iba a interpretar un polinomio no entero.

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

actualización

Véase también la función stepAIC (en el paquete MASS) para automatizar selección de modelos.

+0

¿Cómo puedo conectar Eurequa con R? –

+0

@ adam.888 gran pregunta: no sé la respuesta, pero podría publicarla por separado. Ese último punto fue un poco de una digresión. –

+0

Nota: AIC es el _Akaike Information Criterion_, que premia un ajuste perfecto y penaliza un mayor número de parámetros de un modelo, de una manera que ha demostrado ser óptima en varios sentidos. http://en.wikipedia.org/wiki/Akaike_information_criterion –

37

¿Qué modelo es el "mejor modelo de ajuste" depende de lo que quiere decir con "mejor". R tiene herramientas para ayudar, pero debe proporcionar la definición de "mejor" para elegir entre ellas. Tenga en cuenta los siguientes datos de ejemplo y código:

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 los modelos es el mejor?se podrían hacer argumentos para cualquiera de ellos (pero yo, por mi parte, no querría usar el morado para la interpolación).

5

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

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

Después de usar paso hacia abajo AIC regresión

lm.s <- step(lm.1) 
+2

El uso de 'I (x^2)', etc. no proporciona los polinomios ortogonales apropiados para el ajuste. –

Cuestiones relacionadas