Ajustement du modèle polynomial aux données dans R

Question

J'ai lu les réponses à cette question question et ils sont assez utiles, mais j'ai besoin d'aide surtout en R.

J'ai un exemple d'ensemble de données dans R comme suit :

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)

Je veux ajuster un modèle à ces données pour que y = f(x) . Je veux que ce soit un modèle polynomial d'ordre 3.

Comment puis-je faire cela en R ?

En outre, R peut-il m'aider à trouver le modèle le mieux adapté ?

Answer 1

Pour obtenir un polynôme d'ordre 3 en x (x^3), vous pouvez effectuer les opérations suivantes

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

ou

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

Vous pourriez ajuster un polynôme d'ordre 10 et obtenir un ajustement presque parfait, mais le devriez-vous ?

EDIT : poly(x, 3) est probablement un meilleur choix (voir @hadley ci-dessous).

Answer 2

Le modèle qui est le "modèle le mieux adapté" dépend de ce que vous entendez par "le mieux". R dispose d'outils pour vous aider, mais vous devez fournir la définition de "meilleur" pour choisir entre eux. Considérez l'exemple de données et de code suivant :

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

Lequel de ces modèles est le meilleur ? On pourrait trouver des arguments pour chacun d'entre eux (mais pour ma part, je ne voudrais pas utiliser le modèle violet pour l'interpolation).

Answer 3

En ce qui concerne la question "R peut-il m'aider à trouver le modèle le mieux adapté ?", il existe probablement une fonction pour le faire, en supposant que vous puissiez indiquer l'ensemble des modèles à tester, mais ce serait une bonne première approche pour l'ensemble des polynômes de degré n-1 :

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

Notes

• La validité de cette approche dépendra de vos objectifs, des hypothèses de optimize() y AIC() et si l'AIC est le critère que vous voulez utiliser,

• polyfit() peut ne pas avoir un seul minimum. vérifiez cela avec quelque chose comme :

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

• J'ai utilisé le as.integer() car je ne sais pas comment interpréter un polynôme non entier.

• pour tester un ensemble arbitraire d'équations mathématiques, considérons la Eureqa programme revu par Andrew Gelman aquí

Mise à jour

Voir également le stepAIC (dans le paquet MASS) pour automatiser la sélection des modèles.

Answer 4

La manière la plus simple de trouver le meilleur ajustement dans R est de coder le modèle comme :

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

Après avoir utilisé la régression AIC descendante

lm.s <- step(lm.1)

Answer 5

Par exemple, si nous voulons ajuster un polynôme de degré 2, nous pouvons le faire directement en résolvant un système d'équations linéaires de la manière suivante :

L'exemple suivant montre comment ajuster une parabole y = ax^2 + bx + c en utilisant les équations ci-dessus et le compare avec lm() solution de régression polynomiale. J'espère que cela aidera quelqu'un à comprendre,

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)

x4 <- sum(x^4)
x3 <- sum(x^3)
x2 <- sum(x^2)
x1 <- sum(x)
yx1 <- sum(y*x)
yx2 <- sum(y*x^2)
y1 <- sum(y)

A <- matrix(c(x4, x3, x2,
              x3, x2, x1,
              x2, x1, length(x)), nrow=3, byrow=TRUE)
B <- c(yx2,
       yx1,
       y1)

coef <- solve(A, B)   # solve the linear system of equations, assuming A is not singular
coef1 <- lm(y ~ x + I(x^2))$coef  # solution with lm

coef
# [1] -0.01345808  2.01570523 42.51491582
rev(coef1)
#     I(x^2)       x         (Intercept) 
#     -0.01345808  2.01570523 42.51491582 

plot(x, y, xlim=c(min(x), max(x)), ylim=c(min(y), max(y)+10), pch=19)
xx <- seq(min(x), max(x), 0.01)
lines(xx, coef[1]*xx^2+coef[2]*xx+coef[3], col='red', lwd=3, lty=5)
lines(xx,  coef1[3]*xx^2+ coef1[2]*xx+ coef1[1], col='blue')
legend('topright', legend=c("solve", "lm"),
       col=c("red", "blue"), lty=c(5,1), lwd=c(3,1), cex=0.8,
       title="quadratic fit", text.font=4)