• 12 heures
  • Difficile

Ce cours est visible gratuitement en ligne.

course.header.alt.is_video

course.header.alt.is_certifying

J'ai tout compris !

Mis à jour le 23/07/2019

TP : Pratiquez la régression linéaire multiple sur le jeu de données de l'ozone

Connectez-vous ou inscrivez-vous gratuitement pour bénéficier de toutes les fonctionnalités de ce cours !

Appliquons la régression linéaire multiple à l'échantillon ozone.
Modélisons le pic d'ozone journalier en fonction de toutes les autres variables météorologiques.

Jeu de données ozone
Jeu de données ozone

Importez les données

On importe la librairie ggplot2, qui permettra d'afficher les graphiques :

library(ggplot2)

On importe les données, puis on utilise la commande  lm  pour régresser maxO3 en fonction des autres variables de l'échantillon.

ozone <- read.table("ozone.txt", header=TRUE, sep=";", dec=",")

reg_multi <- lm(maxO3~T9+T12+T15+Ne9+Ne12+Ne15+maxO3v,data=ozone)
summary(reg_multi)

Voici le résultat :

Call:
lm(formula = maxO3 ~ T9 + T12 + T15 + Ne9 + Ne12 + Ne15 + maxO3v, 
    data = ozone)

Residuals:
    Min      1Q  Median      3Q     Max 
-57.768  -7.845  -1.359   8.134  38.984 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 12.70548   13.10860   0.969  0.33467    
T9          -0.63596    1.03462  -0.615  0.54011    
T12          2.50600    1.39946   1.791  0.07625 .  
T15          0.71381    1.13674   0.628  0.53142    
Ne9         -2.76057    0.89157  -3.096  0.00252 ** 
Ne12        -0.37193    1.34590  -0.276  0.78283    
Ne15         0.09028    0.99934   0.090  0.92819    
maxO3v       0.37774    0.06121   6.171 1.32e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 14.43 on 104 degrees of freedom
Multiple R-squared:  0.7546,	Adjusted R-squared:  0.738 
F-statistic: 45.68 on 7 and 104 DF,  p-value: < 2.2e-16

On constate ici que certains paramètres ne sont pas significativement différents de 0, car leur p-valeur n'est pas inférieure à 5 %, le niveau de test que nous souhaitons.

Le $\(\operatorname{R}^2\)$ vaut environ 0.75, et le $\(\operatorname{R}^2\)$ ajusté est d'environ 0.74.

Retirez les variables non significatives

On va donc maintenant retirer les variables non significatives. On commence par la moins significative : Ne15, car elle a une p-valeur de 0.93.

reg_multi <- lm(maxO3~T9+T12+T15+Ne9+Ne12+maxO3v,data=ozone)
summary(reg_multi)
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 12.84916   12.95015   0.992  0.32338    
T9          -0.62976    1.02745  -0.613  0.54124    
T12          2.56018    1.25841   2.034  0.04443 *  
T15          0.65787    0.94878   0.693  0.48960    
Ne9         -2.76526    0.88585  -3.122  0.00232 ** 
Ne12        -0.30796    1.13912  -0.270  0.78742    
maxO3v       0.37752    0.06087   6.202 1.12e-08 ***

On voit alors que c'est maintenant  Ne12  , avec une p-valeur de 0.79, qui est la moins significative. On l'enlève donc.

reg_multi <- lm(maxO3~T9+T12+T15+Ne9+maxO3v,data=ozone)
summary(reg_multi)

On constate qu'il faut maintenant retirer la variable  T9  :

reg_multi <- lm(maxO3~T12+T15+Ne9+maxO3v,data=ozone)
summary(reg_multi)

Et l'on retire ensuite  T15  :

reg_multi <- lm(maxO3~T12+Ne9+maxO3v,data=ozone)
summary(reg_multi)
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.76225   11.10038   0.879    0.381    
T12          2.85308    0.48052   5.937 3.57e-08 ***
Ne9         -3.02423    0.64342  -4.700 7.71e-06 ***
maxO3v       0.37571    0.05801   6.477 2.85e-09 ***

On remarque qu'à présent, tous les paramètres sont significatifs. Quant au $\(R^2\)$ , il vaut environ 0.75, tout comme le $\(R^2\)$ ajusté.

Si l'on souhaite prévoir la concentration journalière en ozone, sachant que la température prévue à 12 h sera de 15 °C, que la valeur de  Ne9  sera de 2, et que la concentration maxO3v de la veille vaut 100, alors on saisit les lignes suivantes :

a_prevoir <- data.frame(T12=15,Ne9=2,maxO3v=100)
maxO3_prev <- predict(reg_multi,a_prevoir)
round(maxO3_prev, digits=2)

On obtient une concentration maxO3 de 84.

Pour aller plus loin : analysez vos résultats

Reprenons la régression linéaire multiple que nous avons obtenue :

library(ggplot2)
ozone <- read.table("ozone.txt",header=TRUE,sep=";",dec=",")
reg_mu1ti <- lm(maxO3~T12+Ne9+maxO3v,data=ozone)
summary(reg_mu1ti)
Call:
lm(formula = maxO3 ~ T12 + Ne9 + maxO3v, data = ozone)

Residuals:
    Min      1Q  Median      3Q     Max 
-56.385  -7.872  -1.941   7.899  41.513 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.76225   11.10038   0.879    0.381    
T12          2.85308    0.48052   5.937 3.57e-08 ***
Ne9         -3.02423    0.64342  -4.700 7.71e-06 ***
maxO3v       0.37571    0.05801   6.477 2.85e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 14.23 on 108 degrees of freedom
Multiple R-squared:  0.752,	Adjusted R-squared:  0.7451 
F-statistic: 109.1 on 3 and 108 DF,  p-value: < 2.2e-16

Nous allons ici réaliser les tests à un niveau $\( \alpha=5\%\)$ :

alpha <- 0.05

Récupérons $\(n\)$, le nombre d'individus de l'échantillon, et $\(p\)$ , le nombre de variables.

n <- dim(ozone)[1]
p <- 4

 Nous allons mener des analyses sur les valeurs atypiques et/ou influentes en travaillant sur un dataframe appelé  analyses  .

analyses <- data.frame(obs=1:n)

Calculez les leviers

On peut calculer les leviers comme ceci, en sachant que le seuil des leviers est de $\(2*\frac{p}{n}\)$ .

analyses$levier <- hat(model.matrix(reg_multi))
seuil_levier <- 2*p/n

On peut visualiser les leviers pour chaque point comme ceci :

ggplot(data=analyses,aes(x=obs,y=levier))+
   geom_bar(stat="identity",fill="steelblue")+
   geom_hline(yintercept=seuil_levier,col="red")+
   theme_minimal()+
   xlab("Observation")+
   ylab("Leviers")+
   scale_x_continuous(breaks=seq(0,n,by=5))
Visualisation des leviers
Visualisation des leviers

Pour sélectionner les points pour lesquels le levier est supérieur au seuil, on exécute ces 2 lignes :

idl <- analyses$levier>seuil_levier
idl
analyses$levier[idl]

Calculez les résidus studentisés

Si l'on souhaite maintenant calculer les résidus studentisés, nous écrivons ceci, sachant que le seuil pour les résidus studentisés est une loi de Student à n-p-1 degrés de liberté :

analyses$rstudent <- rstudent(reg_mu1ti)
seuil_rstudent <- qt(1-alpha/2,n-p-1)

Visualisons les résidus studentisés :

ggplot(data=analyses,aes(x=obs,y=rstudent))+
  geom_bar(stat="identity",fill="steelblue")+
  geom_hline(yintercept=-seuil_rstudent,col="red")+
  geom_hline(yintercept=seuil_rstudent,col="red")+
  theme_minimal()+
  xlab("observation")+
  ylab("Résidus studentisés")+
  scale_x_continuous(breaks=seq(0,n,by=5))

Visualisation des résidus studentisés
Visualisation des résidus studentisés

Déterminez la distance de Cook

Pour trouver la distance de Cook, nous exécutons ceci :

influence <- influence.measures(reg_mu1ti)
names(influence)
colnames(influence$infmat)

Le seuil de la distance de Cook est de n-p :

analyses$dcook <- influence$infmat[,"cook.d"]
seuil_dcook <- 4/(n-p)

On peut détecter les observations influentes comme ceci :

ggplot(data=analyses,aes(x=obs,y=dcook))+
  geom_bar(stat="identity",fill="steelblue")+
  geom_hline(yintercept=seuil_dcook,col="red")+
  theme_minimal()+
  xlab("Observation")+
  ylab("Distance de cook")+
  scale_x_continuous(breaks=seq(0,n,by=5))
Visualisation des observations influentes
Visualisation des observations influentes

Vérifier la colinéarité des variables

Une autre chose à vérifier est l'éventuelle colinéarité approchée des variables :

vif(reg_mu1ti)

Ici, tous les coefficients sont inférieurs à 10, il n'y a donc pas de problème de colinéarité.

Testez l’homoscédasticité

On peut également tester l’homoscédasticité (c'est-à-dire la constance de la variance) des résidus :

bptest(reg_mu1ti)

La p-valeur ici n'est pas inférieure à 5 %, on ne rejette pas l'hypothèse $\(H_0\)$ selon laquelle les variances sont constantes (l'hypothèse d’homoscédasticité).

Testez la normalité des résidus

Si l'on veut tester la normalité des résidus, on peut faire un test de Shapiro-Wilk.

shapiro.test(reg_mu1ti$residuals)

Ici, l'hypothèse de normalité est remise en cause.

Néanmoins, l'observation des résidus, le fait qu'ils ne soient pas très différents d'une distribution symétrique, et le fait que l'échantillon soit de taille suffisante (supérieure à 30) permettent de dire que les résultats obtenus par le modèle linéaire gaussien ne sont pas absurdes, même si le résidu n'est pas considéré comme étant gaussien.

Nous aurions pu aussi sélectionner automatiquement un modèle avec l'ensemble des variables à disposition (variables météo et pic d'ozone de la veille) :

reg_null <- lm(maxO3~1,data=ozone)
reg_tot <- lm(maxO3~T9+T12+T15+Ne9+Ne12+Ne15+maxO3v,data=ozone)
reg_backward <- step(reg_tot,direction="backward")

Et voilà, vous avez vu en pratique comment réaliser une régression linéaire multiple pour déterminer et prédire la concentration d'ozone dans l'atmosphère. Dans la prochaine partie, vous aborderez le modèle de régression logistique.

Exemple de certificat de réussite
Exemple de certificat de réussite