La régression loess (ou "lowess") est une méthode de régression non-paramétrique (c'est-à-dire qu'elle n'est pas associée à une équation, comme par exemple une régression linéaire ou polynomiale classique). Elle permet de produire des courbes lissées, ajustées à un nuage de point, un peu comme ici, par exemple:

Mise en oeuvre

Considérons le jeu de données data.csv:

data=read.table(paste(dat.path,"data.csv",sep=""),sep=";",header=T)
attach(data)

Pour ajuster une régression loess, rien de plus simple... Il suffit de faire appel à la fonction ... loess.

mod=loess(y~x)
print(mod)

## Call:
## loess(formula = y ~ x)
## 
## Number of Observations: 40 
## Equivalent Number of Parameters: 4.4 
## Residual Standard Error: 9.486

On récupère les coordonnées des points de la courbe de la manière suivante:

yfit=predict(mod, newdata=x)

On peut alors représenter le nuage de points et la courbe loess:

plot(x,y)
lines(x,yfit, col="blue",lty=3, lwd=2)

En modifiant le paramètre span de la fonction loess, on peut produire des courbes plus ou moins lissées:

mod1=loess(y~x,span=0.3)
mod2=loess(y~x,span=0.75)
mod3=loess(y~x,span=2)

xfit=seq(from=min(x),to=max(x),length.out=100)
yfit1=predict(mod1,newdata=xfit)
yfit2=predict(mod2,newdata=xfit)
yfit3=predict(mod3,newdata=xfit)
plot(x,y)
points(xfit,yfit1,type="l",lwd=2,col="red")
points(xfit,yfit2,type="l",lwd=2,col="blue")
points(xfit,yfit3,type="l",lwd=2,col="forestgreen")
legend("topleft",c(paste("span=",c(0.3,0.75,2))), lwd=2,lty=1, col=c("red","blue","forestgreen"))

Principe

Dans une régression loess, l'ajustement de la courbe se fait localement. Pour déterminer la valeur y que prend la courbe au point d'abscisse xi, on ajuste un polynôme de degré 1 ou 2 aux points au voisinage de xi. Cet ajustement se fait avec pondération: les points les plus proches de xi ont davantage de poids dans l'ajustement.

Voisinage:

La taille du voisinage est déterminée par un paramètre α. Quand \alpha < 1, le voisinage inclut une proportion α des points du jeu de données.

Les points dans le voisinage d'un xi donné sont représentés en rouge dans l'animation ci-dessus. La fonction loess appelle ce paramètre α span, et par défaut span=0.75.

Quand α ≥ 1, tous les points sont utilisés. L'ajustement ne dépend alors plus de α à travers la détermination du voisinage. En revanche α joue toujours sur le calcul des poids (cf ci-après).

Pondération:

Le poids ρ d'un point d'abscisse xj au voisinage de xi est d'autant plus important que sa distance d_{ij} à x_i est petite.

    \rho=(1-(\frac{d_{ij}}{d_{max}})^3)^3

Quand \alpha<1, d_{max} est la distance maximale (en x) séparant deux points du voisinage considéré. La taille des points rouges dans l'animation ci-dessus est proportionnelle à leur poids.

Quand , \alpha \geq 1, d_{max} est égale à α fois la distance maximale observée (en x) entre deux points du jeu de données.

Ajustement:

Une fois le voisinage et les poids déterminés, une courbe correspondant à un polynôme de degré 1 ou 2 est ajustée pour une série de valeurs xi (ici c'est un polynôme de degré 2, dont la courbe apparaît en rouge dans l'animation ci-dessus). La valeur prédite en ce point xi par la régression loess correspond à la valeur prédite par le polynôme ajusté localement.

Bandes de confiance autour de la courbe

Le calcul de la courbe ajustée peut-être assorti d'une mesure de l'erreur standard:

xfit=seq(from=min(x),to=max(x),length.out=100)
prediction=predict(mod,newdata=xfit, se=TRUE)
yfit=prediction$fit
yfit_se=prediction$se
print(yfit_se[1:5])

## [1] 4.754606 4.354402 3.987376 3.654491 3.356802

On peut ainsi utiliser cette erreur standard pour calculer une bande de confiance (ici à 95%):

plot(x,y)
points(xfit,yfit, col="blue",type="l", lwd=2)
points(xfit,yfit+1.96*yfit_se, col="blue",type="l", lwd=1, lty=2)
points(xfit,yfit-1.96*yfit_se, col="blue",type="l", lwd=1, lty=2)