Réalisez des modélisations de données performantes

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 $\$\backslash (R^2\backslash )\$$ , 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 $\$\backslash (y=x\backslash )\$$ , 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.