Mettons maintenant en œuvre la régression linéaire sur le pic d'ozone, expliqué par la température à midi. Voyons la phase d'estimation du modèle, puis mettons en œuvre une prévision.
Importez les données
On commence par charger la librairie ggplot2, qui permettra d'afficher les graphiques :
library(ggplot2)
On importe ensuite le fichier ozone, qui contient 112 données recueillies dans la ville de Rennes durant l'été 2001.
On trouve dans ce fichier des variables telles que :
MaxO3, qui est la valeur maximale d'ozone observée sur une journée ;
T9, T12 et T15 qui sont les températures prises respectivement à 9 h, 12 h et 15 h ;
Ne9, Ne12, Ne15 qui sont des nébulosités prises à 9 h, 12 h et 15 h ;
Vx9, Vx12 et Vx15 qui sont les composantes est-ouest du vent mesurées à 9 h, 12 h et 15 h ;
MaxO3V, qui donne la teneur maximale en ozone observée la veille ;
vent, l'orientation du vent à 12 h ;
pluie, la présence ou non de pluie.
Visualisez le jeu de données
Visualisons l'ensemble des données avec cette commande.
ozone <- read.table("ozone.txt", header=TRUE, sep=";", dec=",")
On peut représenter graphiquement le nuage de points maxO3 en fonction de T12 :
ggplot(ozone, aes(x=T12, y=maxO3))+
geom_point()+
xlab("T12")+
ylab("maxO3")
Réalisez une régression linéaire simple
Essayons de lancer une régression linéaire simple sur ce nuage de points :
reg_simp <- lm(maxO3~T12,data=ozone)
Voici les résultats en sortie de cette commande :
summary(reg_simp)
Call:
lm(formula = maxO3 ~ T12, data = ozone)
Residuals:
Min 1Q Median 3Q Max
-38.079 -12.735 0.257 11.003 44.671
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -27.4196 9.0335 -3.035 0.003 **
T12 5.4687 0.4125 13.258 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 17.57 on 110 degrees of freedom
Multiple R-squared: 0.6151, Adjusted R-squared: 0.6116
F-statistic: 175.8 on 1 and 110 DF, p-value: < 2.2e-16
Nous obtenons des statistiques sur les résidus, avec le minimum, le maximum et les 3 quartiles, ainsi que des statistiques sur les coefficients obtenus : leur valeur, leur écart-type, la statistique de test de Student, et la p-valeur (le test effectué sur le paramètre est ici le test de significativité : le paramètre vaut 0 versus le paramètre est différent de 0).
Quant au R2 , il est de l'ordre de 0.6. Ce n'est pas très élevé, mais ceci est logique au vu de la dispersion du nuage de points originel.
Visualisez la droite de régression
Alors, voyons à quoi ressemble notre droite.
ggplot(ozone,aes(x=T12,y=maxO3))+
geom_point()+
stat_smooth(method="lm",se=FALSE)+
xlab("T12")+
ylab("MaxO3")
On peut également représenter les valeurs ajustées en fonction des valeurs observées :
ozone$maxO3_ajust_s <- reg_simp$fit
ggplot(ozone, aes(x=maxO3,y=maxO3_ajust_s))+
geom_point()+
geom_abline(intercept=0,slope=1,color="red")+
xlab("MaxO3")+
ylab("MaxO3 ajusté")
La droite qui s'affiche est la première bissectrice. Si le modèle était parfait,
les valeurs réelles et les valeurs ajustées seraient égales, donc sur un tel graphique, les points seraient alignés sur la droite d'équation y=x , soit la première bissectrice.
Représentez les résidus du modèle
On peut obtenir les résidus du modèle à l'aide de cette commande :
ozone$residu_s <- reg_simp$residuals
À partir de ceux-ci, on peut représenter l'histogramme de ces résidus.
ggplot(ozone,aes(x=residu_s))+
geom_histogram(binwidth=10,aes(y=..density..))+
ggtitle("Histogramme")+
xlab("Résidus")+
ylab("")
L'allure de l'histogramme est assez classique : centrée et à peu près symétrique.
Prévoyez la concentration d'ozone
Prévoyons maintenant la concentration en ozone d'une journée. Sachant que la température prévue de cette journée est de 19 °C, on peut utiliser notre modèle de régression à des fins de prévision !
a_prevoir <- data.frame(T12=19)
maxO3_prev <- predict(reg_simp,a_prevoir)
round(maxO3_prev, digits=2)
On obtient une concentration d'ozone d'environ 76.5.
Eh bien, vous avez effectué votre première prédiction grâce à une régression. Félicitations. Le prochain chapitre est un quiz, qui permettra de vous tester sur votre compréhension du modèle.