Fitting modello polinomiale di dati in R

Question

ho letto le risposte a questo question e sono molto utili, ma ho bisogno di aiuto soprattutto in R.

ho un dato esempio dato in R come segue:

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) 

Voglio adattare un modello a questi dati in modo che y = f(x). Voglio che sia un modello polinomiale del 3 ° ordine.

Come posso farlo in R?

Inoltre, può R aiutarmi a trovare il modello più adatto?

Answer 1

Per ottenere un terzo polinomio di ordine in x (x^3), si può fare

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

o

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

Si potrebbe andare bene un 10 ° ordine polinomiale e ottenere una misura quasi perfetta , ma dovresti?

MODIFICA: poly (x, 3) è probabilmente una scelta migliore (vedere @hadley di seguito).

Answer 2

Per quanto riguarda la domanda "R può aiutarmi a trovare il modello più adatto", probabilmente c'è una funzione per farlo, supponendo che si possa indicare l'insieme di modelli da testare, ma questo sarebbe un buon primo approccio per il set di n-1 polinomi di grado:

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

Note

• la validità di questo approccio dipenderà dalla vostra obiettivi, le ipotesi di optimize() e AIC() e se AIC è il criterio che si desidera utilizzare ,

• polyfit() potrebbe non avere un minimo singolo. controllare questo con qualcosa di simile:

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

• ho usato la funzione as.integer() perché non è chiaro per me come vorrei interpretare un polinomio non intero.

• per testare un insieme arbitrario di equazioni matematiche, prendere in considerazione il programma 'Eureqa' recensito da Andrew Gelman here

Aggiornamento

vedere anche la funzione di stepAIC (nel pacchetto MASS) per automatizzare selezione del modello.

Answer 3

Quale modello è il "miglior modello di adattamento" dipende da cosa intendi per "migliore". R ha strumenti per aiutare, ma è necessario fornire la definizione di "migliore" tra cui scegliere. Considera il seguente esempio di dati e codice:

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') 

Quale di questi modelli è il migliore?argomenti potrebbero essere fatti per nessuno di essi (ma io per primo non vorrei usare quello viola per l'interpolazione).

Answer 4

Il modo più semplice per trovare la soluzione migliore in R è quello di codificare il modello come:

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

Dopo aver utilizzato passo verso il basso AIC regressione

lm.s <- step(lm.1)