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.
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 vaut environ 0.75, et le 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 , il vaut environ 0.75, tout comme le 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-16Nous allons ici réaliser les tests à un niveau :
alpha <- 0.05Récupérons , le nombre d'individus de l'échantillon, et , 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 .
analyses$levier <- hat(model.matrix(reg_multi))
seuil_levier <- 2*p/nOn 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))
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))
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))
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 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.
