Régressions polynômiales avec R

Il s'agit simplement d'utiliser un modède linéaire en fournissant une formule polynômiale à R.

Pour rappel la solution à la minimisation de (ou X est une matrice et i et Formula : $\beta$beta des vecteurs colonne de dimensions appropriées)

Formula : $$ \mbox{argmin}_\beta( ||y - X\beta||_2) $$

est sous réserve d'inversibilité de Formula : $^tX X$

Formula : $$ \beta = \left(^tX X\right)^{-1} {}^t Xy$$

Il s'agit du calcul effectué par R lors de l'appel à la fonction lm. Toute l'astuce consiste à construire une matrice X de manière adaptée.

Supposons que l'on veuille approximer le jeu de données généré par la fonction suivante :

 
#Random seed initialisée avec une constante pour la reproductibilité des résultats
set.seed(0);
#echantillonnage dans -10;10
x <- runif(300,  min=-10, max=10) 

#fonction objectif (supposée inconnue et à découvrir)  bruitée
y <- 0.1*x^3 - 0.5 * x^2 - x + 10 + rnorm(length(x),0,5) 

#Affichage des données
plot(x,y,col=rgb(0.3,0.4,0.1,0.5),pch=16) 

#Affichage de la solution idéale (nécessite de trier x pour tracer une ligne) 
ix <- sort(x,index.return=T)$ix
lines(x[ix],0.1*x[ix]^3 - 0.5 * x[ix]^2 - x[ix] + 10,col=2,lwd=2)

R plot

Que l'on peut fitter par un simple :

model <- lm( y ~ x + I(x^2) + I(x^3) ) 
print(model)

Call:
lm(formula = y ~ x + I(x^2) + I(x^3))

Coefficients:
(Intercept)            x       I(x^2)       I(x^3)  
    10.3762      -1.0775      -0.5065       0.1023  

Les I() sont importants car ils servent à éviter l'interprétation de x^2. Si cela n'est pas clair, essayez la même commande sans les I().

Il est parfois préférable pour des raisons de travailler avec des polynômes orthogonaux (cela a tendance à réduire la corrélation entre les coefficients obtenus par la régression) dans ce cas on peut utiliser comme formule poly(x,3). Toutefois l'interprétation des coefficients ainsi obtenus n'est pas évidente.

Enfin, nous pouvons utiliser la fonction predict pour visualiser la courbe de la régression effectuée par notre modèle :

 
#Random seed initialisée avec une constante pour la reproductibilité des résultats
set.seed(0);
#echantillonnage dans -10;10
x <- runif(300,  min=-10, max=10) 

#fonction objectif (supposée inconnue et à découvrir)  bruitée
y <- 0.1*x^3 - 0.5 * x^2 - x + 10 + rnorm(length(x),0,5) 

plot(x,y,col=rgb(0.3,0.4,0.1,0.5),pch=16) 

myPredict <- predict( lm(y ~ x + I(x^2) + I(x^3)) ) 
ix <- sort(x,index.return=T)$ix
lines(x[ix], myPredict[ix], col=2, lwd=2 )  

R plot

Les erreurs en fonction de x (possible car on connait la fonction objectif)

R plot

On peut aussi en faisant un plot sur le résultat de lm obtenir une série de courbes de contrôle sur les résidus. Ceux-ci sont important en particulier lorsque l'on ne peut pas avoir de représentation graphique des données (à cause de la dimension par exemple). Ainsi par exemple si les résisdus présentent un biais c'est probablement que l'on a manqué une variable explicative. De même si certains points présentent une distance de Cook élevée alors il s'agit possiblement de valeurs aberrantes (en tout cas de points nécessitant une attention particulière).

R plot

Enfin on peut utiliser le résultat en prédiction pour de nouvelles données avec des intervalles de confiance (avec deux modes de calcul l'un représentant un intervalle sur la position de la courbe idéale l'autre sur la position du point que l'on va observer en tenant compte du bruit):

 
#Random seed initialisée avec une constante pour la reproductibilité des résultats
set.seed(0);
#echantillonnage dans -10;10
x <- runif(300,  min=-10, max=10) 

#fonction objectif (supposée inconnue et à découvrir)  bruitée
y <- 0.1*x^3 - 0.5 * x^2 - x + 10 + rnorm(length(x),0,5) 

# pour faire un "trou" dans l échantillonnage et se tester en prédiction 
sx <- x; sy <- y
useOnly <- x<0 | x>5 
x <- x[useOnly]
y <- y[useOnly]

# Les données à prédire
new <- data.frame(x = seq(-2, 6 , 0.5))

# prediction (rappel : les données pour x dans 0;5 sont manquantes )
predict( lm(y ~ x + I(x^2) + I(x^3)), new,  se.fit = T) 
pred.interval <- predict(lm(y ~ x + I(x^2) + I(x^3)), new, interval="prediction")
pred.confiance <- predict(lm(y ~ x + I(x^2) + I(x^3)), new, interval="confidence")
matplot(new$x,cbind(pred.confiance, pred.interval[,-1]),lwd = 2,
 lty=c(1,2,2,3,3), col=c(1,4,4,3,3), type="l", ylab="predicted y")

#la courbe idéale (en rouge)
lines(new$x,0.1*new$x^3 - 0.5 * new$x^2 - new$x + 10, col = 2 )

#avec les points présents dans les données (et non utilisés)
points(sx,sy,col=rgb(0.3,0.4,0.1,0.5),pch=16);

R plot

$fit
        1         2         3         4         5         6         7         8         9        10        11 
 9.464479 10.216388 10.484081 10.345701  9.879393  9.163298  8.275560  7.294321  6.297726  5.363917  4.571036 
       12        13        14        15        16        17 
 3.997228  3.720635  3.819400  4.371666  5.455577  7.149275 

$se.fit
        1         2         3         4         5         6         7         8         9        10        11 
0.5326349 0.5471408 0.5693318 0.5989580 0.6347678 0.6747961 0.7166866 0.7579460 0.7961075 0.8288235 0.8539198 
       12        13        14        15        16        17 
0.8694408 0.8737122 0.8654501 0.8439623 0.8095255 0.7641028 

$df
[1] 224

$residual.scale
[1] 5.17617

Et après ?

Il faut garder à l'esprit que tout ceci marche bien car un certain nombre de conditions sont réunies :

  • Nombre de points observés "suffisant"
  • L'on recherche dans la bonne classe d'hypothèse : on cherche a approximer un polynôme de degré 3 bruité uniformément avec un modèle qui correspond (essayez de faire la même chose avec un polynôme de degré 2 ou 4).
  • La dimension des données (ici 1) est petite. Avec une régression en dimension >7 tout devient beaucoup plus difficile (malédiction de la dimension) sauf si les données vivent en réalité dans un sous espace de dimension plus faible (lien avec la PCA par exemple).
Design selector