Université Cadi Ayyad, 2009 - 2010
Faculté des Sciences et Techniques,
Marrakech.
Travaux Pratiques
Initiation à R
Le but de ce TP est d’illustrer la régression linéaire simple à partir des variables du fichier
[Link]. Ce fichier contient en particulier deux variables : la concentration en ozone
notée O3 et la température notée T 12, et l’on souhaite connaître la relation entre ces deux
quantités.
1 Quelques rappels sur la régression linéaire simple
Dans le cas Gaussien, le modèle de régression linéaire simple pour des observations
{(xi , yi ) i = 1, . . . , n} s’écrit :
yi = β0 + β1 xi + ϵi , i = 1, . . . , n
où les erreurs ϵi , i = 1, . . . , n sont supposées i.i.d N (0, σ 2 ). Pour l’estimation des paramètres
β0 et β1 par moindres carrées, on définit les quantités suivantes :
n n n n
X X 1 X 1 X
x̄ = xi , ȳ = yi , s2x = (xi − x̄)2 , s2y = (xi − ȳ)2 ,
i=1 i=1
n − 1 i=1 n − 1 i=1
n
1 X sxy
sxy = (xi − x̄)(yi − ȳ) et r = .
n − 1 i=1 sx sy
Les estimateurs des moindres carrés sont alors donnés par :
sxy
β̂1 = et β̂0 = ȳ − β̂1 x̄.
s2x
On définit ensuite :
• les valeurs estimées : ŷi = β̂0 + β̂1 xi ,
• les résidus : êi = yi − ŷi ,
1
n
2 1 X 2
• la variance estimée : σ̂ = ê .
n − 2 i=1 i
1/2 2 1/2
x̄2
Posons s0 = σ̂ n1 + (n−1)s 2 et s 1 = σ̂ x̄
(n−1)s2x
, on peut alors tester l’hypothèse de
x
nullité des coefficients β0 et β1 à partir de la loi des statistiques suivantes :
β̂0 − β0 β̂1 − β1
∼ t(n − 2) et ∼ t(n − 2).
s0 s1
où t(n − 2) désigne la loi de Student à (n − 2) degrés de liberté. Afin de déterminer la qualité
du modèle, on définit les quantités suivantes :
• variance totale : SST = (n − 1)s2y ,
s2xy
• variance expliquée : SSR = (n − 1) ,
s2x
• variance résiduelle : SSE = (n − 2)σ̂ 2 .
On peut vérifier que SST = SSR + SSE. Le coefficient de détermination est alors donné
par :
2
s2xy SSR
r = 2 2 = .
sx sy SST
On rappelle que plus r est proche de 1, plus le modèle peut être considéré comme bien adapté
aux données.
2 La concentration en ozone
Nous allons traiter les 50 données journalières de concentration en ozone. La variable à
expliquer est la concentration en ozone notée O3 et la variable explicative est la température
notée T 12.
Pour une régression simple, nous commençons toujours par représenter les données.
> ozone <- [Link]("[Link]",header=T,sep=";")
> plot(ozone$T12,ozone$O3,xlab="T12",ylab="O3")
2
140
●
●
● ●
●
●
120
●
●
●
●
●
● ●
100
● ●
O3
● ● ●
● ● ● ● ●
●
● ● ●
● ●
80
● ● ●
● ●
●
● ●
● ●
●
● ●●
●
60
● ●
●
●
●
●
40
10 15 20 25 30
T12
Figure 1: 50 données journalières de température et O3.
Remarque : Nous obtenons un graphique qui permet de vérifier visuellement si une régres-
sion linéaire est pertinente. Autrement dit il suffit de regarder si le nuage de point s’étire le
long d’une droite. Bien qu’ici il semble que le nuage s’étire sur une première droite jusqu’à
22 ou 23 °C puis selon une autre droite pour les hautes valeurs de températures, nous pouvons
tenter une régression linéaire simple.
Nous effectuons ensuite la régression linéaire, c’est-à-dire la phase d’estimation.
> reg <-lm(O3∼T12,data=ozone)
Afin de consulter les résultats, nous effectuons
> resume <- summary(reg)
> resume
Remarque : Les sorties du logiciel donnent une matrice (sous le mot Coefficients) qui com-
porte pour chaque paramètre (chaque ligne) 5 colonnes.
(1) estimations des paramètres
(2) les écarts-types estimés des paramètres
(3) la valeur observée de la statistique de test d’hypothèse H0 : βi = 0
(4) P (tn−2 > |t|); H0 est rejetée au niveau α si P < α.
La dernière colonne est une version graphique du test :
*** signifie que le test rejette H0 pour α ≥ 0.001
** signifie que le test rejette H0 pour α ≥ 0.01
* signifie que le test rejette H0 pour α ≥ 0.05
. signifie que le test rejette H0 pour α ≥ 0.1.
3
Comparez les sorties de « anova(reg)» avec «summary(reg)».
Représentation des données et de la droite ajustée
> plot(ozone$T12,ozone$O3,xlab="T12",ylab="O3")
> # le paramètre pch permet de changer le motif représenté pour
chaque point
> abline(reg)
140
●
●
● ●
●
●
120
●
●
●
●
● ●
100
● ●
O3
● ● ●
● ● ● ● ●
●
● ● ●
● ●
80
● ● ●
● ●
●
● ●
● ●
●
● ●●
●
60
● ●
●
●
●
●
40
10 15 20 25 30
T12
Figure 2: 50 données journalières de température et O3 et la droite ajusté.
Afin d’examiner la qualité du modèle et des observations, nous traçons la droite ajustée et
les observations. Comme il existe une incertitude dans les estimations, nous traçons aussi un
intervalle de confiance de la droite (à 95 %).
> plot(O3∼T12,data=ozone)
> T12=seq(min(ozone[,"T12"]),max(ozone[,"T12"]),length=100)
> grille <- [Link](T12)
> ICdte <- predict(reg,new=grille,interval="confidence",level=0.95)
> matlines(grille$T12,cbind(ICdte), lty=c(1,2,2),col=1)
4
140
●
●
● ●
●
●
120
●
●
●
●
●
● ●
100
● ●
O3
● ● ●
● ● ● ● ●
●
● ● ●
● ●
80
● ● ●
● ●
●
● ●
● ●
●
● ●●
●
60
● ●
●
●
●
●
40
10 15 20 25 30
T12
Figure 3: 50 données journalières de température et O3 et l’intervalle de confiance pour
E(Y ).
Dans une optique de prévision, il est nécessaire de s’intéresser à la qualité de prévision.
Cette qualité peut être envisagée de manière succincte grâce à l’intervalle de confiance des
prévisions. Afin de bien le distinguer de celui de la droite, nous figurons les deux sur le
même graphique.
> plot(O3∼T12,data=ozone,ylim=c(0,150))
> T12 <- seq(min(ozone[,"T12"]),max(ozone[,"T12"]),length=100)
> grille <- [Link](T12)
> ICdte <- predict(reg,new=grille,interval="conf",level=0.95)
> ICprev <- predict(reg,new=grille,interval="pred",level=0.95)
> matlines(T12,cbind(ICdte,ICprev[,-1]),lty=c(1,2,2,3,3),col=1)
> legend(8,145,lty=2:3,c("E(Y )","Y "))
150
●
E(Y) ●
Y ● ●
●
● ●
●
●
●
●
● ●
100
● ●
● ● ●
● ● ● ● ●
●
● ●
● ● ●
● ● ●
O3
● ●
● ● ● ● ●
● ●
●●
● ● ●
●
●
50
●
0
10 15 20 25 30
T12
Figure 4: Droite de régression et intervalles de confiance pour Y et pour E(Y ).
5
Si nous nous intéressons au rôle des variables, nous pouvons calculer les intervalles de con-
fiance des paramètres.
> seuil <- qt(0.975,df=reg$[Link])
> beta0min <- coef(resume)[1,1]-seuil*coef(resume)[1,2]
> beta0max <- coef(resume)[1,1]+seuil*coef(resume)[1,2]
> c(beta0min,beta0max)
> beta1min <- coef(resume)[2,1]-seuil*coef(resume)[2,2]
> beta1max <- coef(resume)[2,1]+seuil*coef(resume)[2,2]
> c(beta1min,beta1max)
Pour aller plus loin, il est possible de tracer la région de confiance simultanée des deux
paramètres, ce qui est rarement fait en pratique. Nous pouvons la comparer aux intervalles de
confiance au même degré de confiance. Cette comparaison illustre uniquement la différence
entre intervalle simple et région de confiance. En général l’utilisateur de la méthode choisit
l’une ou l’autre forme. Pour cette comparaison, nous utilisons les commandes suivantes :
> library(ellipse)
> plot(ellipse(reg,level=0.95),type="l",xlab="beta0",ylab="beta1")
> points(coef(reg)[1], coef(reg)[2],pch=3)
> # comparaison avec IC
> lines(c(beta0min,beta0min,beta0max,beta0max,beta0min),
c(beta1min,beta1max,beta1max,beta1min,beta1min),lty=2)
1.0 1.5 2.0 2.5 3.0 3.5 4.0
beta1
0 10 20 30 40 50 60
beta0
Figure 5: Région de confiance simultanée des deux paramètres.
Représentation des résidus en fonction de la variable explicative
6
> plot(T12,residuals(reg),pch=16,xlab="T12" ,ylab="résidus de la
régression linéaire")
40
●
● ●
● ●
●
résidus de la régression linéaire
● ●
20
● ● ●
● ● ●
● ●
●
●
● ●
● ●
●
0
●
● ● ● ●
● ● ● ●
● ●
●
●
● ●
●
−20
●
● ● ●
●●
● ●
●
●
−40
10 15 20 25 30
T12
Figure 6: Les résidus en fonction de T 12.
Représentation graphique par défaut (examen des résidus et des données influentes)
> par(mfrow=c(2,2)) # permet d’afficher 2×2 graphes sur la même
fenêtre
> plot(reg)
Residuals vs Fitted Normal Q−Q
Standardized residuals
● 19980504 19980504 ●
−2 −1 0 1 2
● 19970422 ● 19970422● ●
● ● ● ●●●
●
0 20
● ●●
●● ● ● ● ● ●
Residuals
●●●
●●●
● ● ●
●●
● ●● ●
●
● ●● ●
● ● ●
●●
● ● ●● ● ● ● ●
●
● ● ● ●●
●●
●●
●
● ●● ●
●●
●● ●●
●
● ●
●● ●
● ● ● ● ●●
●●
●
●●●
−40
● ●
●
19960627 ●
● 19960627
60 70 80 90 110 −2 −1 0 1 2
Fitted values Theoretical Quantiles
Scale−Location Residuals vs Leverage
0.0 0.5 1.0 1.5
Standardized residuals
Standardized residuals
19960627 ●
● 19980504 ● 19980504 0.5
0 1 2
● 19970422● ● ●● 19970422
● ●
● ● ● ● 19960506 ●
● ● ●
●● ●
●● ● ● ● ● ●
●● ● ●
●
● ● ●● ● ● ●
●● ●
● ●
● ●
●● ●● ●
●
● ●●
●
● ● ●●●● ● ● ●
● ● ●
●● ● ● ●
● ● ● ●● ●
● ● ● ●
●
●
● ●● ●●●
● ●
−2
● ● Cook's distance 0.5
60 70 80 90 110 0.00 0.05 0.10 0.15
Fitted values Leverage