La régression linéaire simple permet de modéliser la relation linéaire entre deux variables quantitatives dans un objectif explicatif et/ou prévisionnel. Nous disposons d'une variable explicative, notée $x$, et une variable à expliquer, notée $y$, liées par le modèle suivant : $$y=\beta_0+\beta_1x+\varepsilon,$$ où $\varepsilon$ est la variable de bruit ou erreur de mesure. Nous supposons que $\varepsilon$ est une variable aléatoire centrée de variance $\sigma^2$, et éventuellement de loi normale $\mathcal N(0,\sigma^2)$. Nous disposons de $n$ couples $(x_1,y_1),\dots,(x_n,y_n)$ afin d'estimer les paramètres inconnus $\beta_0$ et $\beta_1$. Notons $X=(x_1,\dots,x_n)^T$, $A=[\mathbb 1,X]$ et $Y=(y_1,\dots,y_n)^T$. L'estimateur $(\hat\beta_0,\hat\beta_1)$ par la méthode des moindres carrés est donné par $$(\hat\beta_0,\hat\beta_1)^T=(A^TA)^{-1}A^TY.$$ D'après le cours, on a $$\hat\beta_0=\bar Y-\hat\beta_1\bar X,\qquad \hat\beta_1=\frac{\mathrm{Cov}(X,Y)}{\mathrm{Var}(X)}.$$ Une fois les paramètres estimés, on obtient la droite de régression $$f(x) = \hat\beta_0+\hat\beta_1 x,$$ ce qui permet d'effectuer des prévisions pour une nouvelle variable $x^*$ par $$y^P=f(x^*)=\hat\beta_0+\hat\beta_1 x^*.$$
Les valeurs ajustées sont définies par $$\hat y_i=f(x_i)=\hat\beta_0+\hat\beta_1 x_i,$$ et les résidus estimés par $$\hat\varepsilon_i=y_i-\hat y_i.$$
Voyons sur les données cathedral comment effectuer une régression linéaire simple avec R.
Importons d'abord les données :
cath <- read.table('cathedral.txt',header=TRUE)
names(cath)
attach(cath)
Nous aimerions expliquer la hauteur de cathédrale en fonction de la longueur. Calculons d'abord la corrélation entre variables longueur et hauteur, mesure de la colinéarité entre variables :
cor(haut,long)
La corrélation est positive, mais pas tout près de 1. Traçons le nuage des points pour vérifier davantage la pertinence d'un modèle de régression linéaire pour nos variables. Pour cela deux commandes sont possibles qui produisent exactement le même résultat : soit plot(long,haut) (commande que vous connaissez déjà), soit plot(haut~long) (attention à l'ordre des variables). Essayez les deux commandes :
plot(long,haut)
plot(haut~long)
Pour calculer les estimateurs par moindres carrés dans le modèle linéaire, c'est extrêmement facile car la fonction lm (=linear model) fait tout. Le plus difficile est de savoir lire la sortie :
reg.simple <- lm(haut~long)
summary(reg.simple)
La matrice Coefficients comporte pour chaque paramètre (chaque ligne) quatre colonnes : son estimation (colonne Estimate), son écart-type estimé (Std. Error), la valeur observée de la statistique du test d'hypothèse $H_0:\beta_k=0$ contre $H_1:\beta_k\neq 0$. Enfin, la probabilité critique ou p-valeur (Pr(>|t|)) donne, pour la statistique de test dous $H_0$, la probabilité de dépasser la valeur estimée et permet de conclure.
Les coefficients $\beta_0$ et $\beta_1$ sont estimés par $37.53$ et $0.0874$. Les tests de significativité des coefficients donnent ici des p-valeurs inférieures à $0.001$. Ainsi l'hypothèse nulle de chacun des tests est rejetée au profit de l'hypothèse alternative. La très petite p-valeur pour la constante indique que la constante (l'intercept) doit apparaître dans le modèle. La très petite p-valeur pour la pente indique une liaison linéaire significative entre la hauteur et la longueur de cathédrale.
Le résumé de l'étape d'estimation fait figurer la valeur de l'estimateur $\hat\sigma$ de l'écart-type résiduel $\sigma$, qui vaut ici $11.75$, ainsi que le nombre de degrés de liberté associé $n-2=23$.
Nous pouvons consulter la liste des différents résultats de l'objet reg.simple et summary(reg.simple) avec :
names(reg.simple)
names(summary(reg.simple))
On peut alors récupérer les coefficients avec :
reg.simple$coeff
ou par la fonction
coef(reg.simple)
Pour ajuster un modèle sans la constante, on procède de la manière suivante :
reg.ss.simple <- lm(haut~long-1)
summary(reg.ss.simple)
Pour tracer la droite de régression et la superposer au nuage des points, nous pouvons simplement appliquer la commande abline(reg.simple). De plus, on peut rajouter au graphique les valeurs ajustées (disponible dans reg.simple\$fitted.values) :
plot(haut~long)
abline(reg.simple,col='red')
points(long,reg.simple$fitted.values,pch=15,col='blue')
Les résidus estimés $\hat\varepsilon_i$ sont disponibles dans le vecteur reg.simple\$residuals.
On obtient les résidus studentisés (par validation croisées) $$t_i^*=\frac{\hat\varepsilon_i}{\hat \sigma_{(i)} \sqrt{1-h_{ii}} },\qquad i=1,\dots,n$$ par la fonction rstudent :
residu.simple <- rstudent(reg.simple)
Afin d'identifier d'éventuelles valeurs aberrantes dans les données, on utilise la représentation graphique suivante :
n <- length(haut)
plot(1:n,residu.simple,col='blue',xlab='Index',ylab='Résidus studentisés')
abline(-2,0)
abline(2,0)
Quand l'échantillon ne contient pas de valeurs aberrantes, 95% des résidus studientisés $t_i^*$ se trouvent dans l'intervalle $[-2,2]$ et les autres 5% pas loin de ces limites. C'est le cas dans notre exemple puisque un résidu seulement (sur 25) se trouvent à l'extérieur de l'intervalle $[-2,2]$.
Afin d'identifier des éventuelles observations trop influentes dans le jeu de données, on analyse les poids $h_{ii}$ des observations que l'on obtient par l'instruction suivante :
levier <- hatvalues(reg.simple)
On compare ces points aux seuils $2p/n$ et $3/pn$ par un graphique :
plot(1:n, levier,xlab='Index',ylab='Poids h_ii')
p <- reg.simple$rank
seuil1 <- 2*p/n
seuil2 <- 3*p/n
abline(seuil1,0,lty=2)
abline(seuil2,0,lty=3)
On observe qu'il y a une seule observation qui dépasse les deux seuils, et deux autres au-dessus du premier seuil.
Afin de connaître précisément les observations suspectes on peut rajouter des labels dans le graphique. En effet, les individus suspects ont les indices ID <- (1:n)[levier>seuil1]. On peut les mettre en avant dans le graphique par l'instruction text(ID,levier[ID],ID,col='red',pos=1). Essayez-le.
On peut analyser ces points davantage en calculant leur distance de Cook.
La distance de Cook est une autre mesure pour l'influence d'un individu sur l'estimation de $\beta$. Pour l'observation $i$ elle vaut $$C_i=\frac{h_{ii}}{p(1-h_{ii})}\frac{\hat\varepsilon_i^2}{\hat\sigma^2}.$$ Sous R on l'obtient par la fonction cooks.distance :
cook <- cooks.distance(reg.simple)
Le seuil critique pour la distance de Cook à partir duquel on considère que l'observation est trop influente est le quantile $f_{p,n-p}(0.5)$. Une distance de Cook en-dessous de $f_{p,n-p}(0.1)$ est considéré comme souhaitable. Traçons le graphique pour notre exemple :
plot(1:n, cook, xlab='Index',ylab='Distance de Cook')
s1 <- qf(0.5,p,n-p)
s2 <- qf(0.1,p,n-p)
abline(s2,0,lty=2)
abline(s1,0,lty=3)
On voit qu'aucune des distance de Cook dépasse le seuil qui serait préoccupant.
Comme pour les points leviers, rajouter au graphique les labels des observations suspectes. Comparez avec les valeurs aberrantes et les points leviers.
Ayant une nouvelle observation $x_\text{new}$, il suffit d'utiliser les estimations pour prévoir la valeur de $Y$ correspondante. Cependant, la valeur prédite est de peu d'intérêt sans l'intervalle de confiance associé. Voyons cela sur un exemple. Nous disposons de la valeur de la longueur d'une nouvelle cathédrale qui est de 312 pieds.
Il y a deux façons de calculer la valeur prédite $\hat y^p=\hat\beta_0+\hat\beta_1*x_\text{new}$ :
# à la main :
xnew <- 312
val.pred <- coef(reg.simple)[1]+coef(reg.simple)[2]*xnew
val.pred
# par la fonction predict :
xnew <- data.frame(long=xnew)
predict(reg.simple,xnew,interval="pred")
La fonction predict nécessite en argument un data.frame avec les mêmes noms de variables explicatives (ici long) que les données de départ. En revanche, la fonction predict ne fournit pas seulement la valeur prédite (dans la colonne fit), mais également l'intervalle de prédiction associée de niveau $0.95$. (L'argument level permet de modifier le niveau de confiance de l'intervalle.)
La fonction predict permet aussi de calculer des intervalles de confiance pour $x^T\beta$, en utilisant interval="conf". Par exemple :
predict(reg.simple,xnew,interval='conf',level=0.9)
Parmi les différentes substances chimiques qui polluent l'air, on compte l'ozone qui a un impact préoccupant sur la santé. Afin de mieux comprendre les conditions météorologiques qui favorisent une teneur élevée d'ozone dans l'air, nous allons analyser un jeu de donnes qui contient 13 variables au total dont la concentration en ozone pendant la journée et des informations sur la température, la nébulosité, le vent, la pluie. Nous disposons de 112 données relevées durant l'été 2001 à Rennes.
Dans ce TP nous nous limitons à analyser la concentration maximale en ozone par jour en fonction de la température prévue à midi. Ce sont les variables de nom maxO3 et T12 dans le fichier ozone.txt. Ainsi, nous considérons le modèle $$\text{maxO3} = \beta_0 + \beta_1*\text{T12} + \varepsilon.$$
Les autres variables serons analyser au TP prochain.