# Chargez les packages
library(readxl)
library(gtools)
library(ggplot2)
library(corrplot)
library(car)
library(gtsummary)
library(GGally)
library(ggstatsplot)
library(sjPlot)Modèles linéaires
Avant-propos
Si ce n’est pas déjà fait, chargez et préparez la table :
# Importez la table (disponible ici : https://gitlab.huma-num.fr/tfeuillet/cours/data)
df <- read_excel("data/table_climato.xlsx")
# Discrétisez la variable "ampl" en terciles et renommez les niveaux
df$amplDisc <- quantcut(df$ampl, 3)
levels(df$amplDisc) <- c("faible","moyenne","forte")
# Discrétise la variable "jrs_neige" en quartiles et renommez les niveaux
df$neigeDisc <- quantcut(df$jrs_neige, 4)
levels(df$neigeDisc) <- c("faible","moyenne","forte","très forte") 1 Régression linéaire simple
1.1 Principes
Il est maintenant temps de passer à la régression linéaire simple, également appelée MCO (moindres carrés ordinaires). Cette méthode consiste à estimer les termes de l’équation suivante :
\(ampl_i = a \times distLitt_i + b + e_i\)
où ampl est l’amplitude thermique annuelle (la variable que l’on cherche à expliquer), a le coefficient associé à la distance au littoral (distLitt), b la constante (intercept) et e les résidus. Cette équation est celle de la droite qui ajuste au mieux le nuage de points des individus (qui sont ici des villes, notées i dans l’équation).
On peut aussi représenter cette droite de régression sur le nuage de points :
plot1 <- ggplot(df, aes(x = dist_litt, y = ampl)) + geom_point() +
labs(title = "Relation amplitude thermique vs distance au littoral")+
xlab("Distance au littoral (km)") + ylab("Amplitude thermique (degrés)")+
theme(plot.title = element_text(face = "bold"))
plot1 + geom_point(col="blue") + geom_smooth(method = "lm", color="red", se=TRUE) + theme_grey()`geom_smooth()` using formula = 'y ~ x'
Ou avec le package car
car::scatterplot(ampl ~ dist_litt, data = df, id = list(n=4)) 
[1] 6 7 9 16# Avec l'identification des 4 principaux outliers
L’intérêt est de prédire les valeurs de Y (amplitude thermique) en fonction des valeurs de la variable explicative (X). La régression linéaire est dite simple quand il n’y a qu’une seule variable explicative.
1.2 Estimation du modèle
# Estimation du modèle avec la fonction lm (pour linear model)
lm1 <- lm(data = df, formula = ampl ~ dist_litt)
summary(lm1)
Call:
lm(formula = ampl ~ dist_litt, data = df)
Residuals:
Min 1Q Median 3Q Max
-2.77968 -0.30598 0.09252 0.45291 1.21509
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.357e+01 2.425e-01 55.96 < 2e-16 ***
dist_litt 8.819e-03 8.186e-04 10.77 5.26e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7934 on 31 degrees of freedom
Multiple R-squared: 0.7892, Adjusted R-squared: 0.7824
F-statistic: 116.1 on 1 and 31 DF, p-value: 5.258e-12
La fonction summary() présente les estimations de façon brute, mais d’autres fonctions permettent des rendus plus esthétiques, comme avec le package gtsummary :
gtsummary::tbl_regression(lm1)| Characteristic | Beta | 95% CI | p-value |
|---|---|---|---|
| dist_litt | 0.01 | 0.01, 0.01 | <0.001 |
| Abbreviation: CI = Confidence Interval | |||
1.3 Interprétation des résultats
On obtient un \(R^{2}\) (coefficient de détermination) de 0.78 : 78% de la variation de l’amplitude thermique est expliquée par la distance au littoral, ce qui est déjà très satisfaisant. Le coefficient de régression associé à la distance au littoral est égal à 0.008819. Cela signifie que quand on s’éloigne de 1 km du littoral, l’amplitude thermique augmente en moyenne de 0.008819 °C (soit presque 1 °C tous les 100 km). La p-value associée à ce coefficient (notée Pr(>|t|) dans le tableau de résultats) est très proche de 0. Cela signifie que l’effet de la distance au littoral sur l’amplitude thermique est statistiquement significatif (on n’a quasiment aucune chance de se tromper en affirmant cela).
1.4 Analyse des résidus (diagnostic)
Il convient ensuite d’examiner les résidus de la régression.
1.4.1 Normalité des résidus
Ces fonctions ci-dessous permettent de représenter graphiquement les résidus studentisés (ie standardisés) pour caractériser leur normalité et leur grandeur. Les valeurs supérieures à 2 sont considérées comme fortes :
par(mfrow=c(1,3))
plot(rstudent(lm1))
hist(rstudent(lm1), breaks = 20, col="darkblue", border="white", main="Analyse visuelle des résidus")
qqPlot(lm1)
[1] 7 16
1.4.2 Outliers
On constate que 2 stations s’éloignent assez nettement du modèle : la station 7 surtout, et la station 16. Identifions ces deux stations :
df[c(7,16),]# A tibble: 2 × 11
nom precip alt temp ampl jrs_neige dist_litt X Y amplDisc
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
1 Cherbourg 896 15 11.3 10.8 5.1 1 366208. 6.96e6 faible
2 Marseille 588 50 14.2 16.8 1.3 532 891476. 6.25e6 forte
# ℹ 1 more variable: neigeDisc <fct>Il s’agit donc de Cherbourg et de Marseille. Cherbourg présente une amplitude particulièrement faible par rapport aux autres stations littorales. Marseille a également une amplitude relativement faible par rapport à se distance (importante) à l’océan.
Voyons les conséquences si on relance la régression sans l’individu 7 (Cherbourg) :
lm2 <- update(lm1, subset=-7)
summary(lm2)
Call:
lm(formula = ampl ~ dist_litt, data = df, subset = -7)
Residuals:
Min 1Q Median 3Q Max
-1.32725 -0.32533 0.05123 0.33192 1.04749
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.386e+01 1.943e-01 71.32 < 2e-16 ***
dist_litt 8.028e-03 6.459e-04 12.43 2.3e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6054 on 30 degrees of freedom
Multiple R-squared: 0.8374, Adjusted R-squared: 0.832
F-statistic: 154.5 on 1 and 30 DF, p-value: 2.305e-13par(mfrow=c(1,2))
plot(rstudent(lm2))
qqPlot(lm2)
6 16
6 15 Le \(R^{2}\) passe de 0.78 à 0.83. On peut aussi comparer les coefficients de deux modèles :
car::compareCoefs(lm1, lm2, pvals = TRUE)Calls:
1: lm(formula = ampl ~ dist_litt, data = df)
2: lm(formula = ampl ~ dist_litt, data = df, subset = -7)
Model 1 Model 2
(Intercept) 13.571 13.856
SE 0.242 0.194
Pr(>|z|) <2e-16 <2e-16
dist_litt 0.008819 0.008028
SE 0.000819 0.000646
Pr(>|z|) <2e-16 <2e-16
La différence est mince, l’effet de la distance au littoral est à peu près identique. Notez qu’il peut être délicat de supprimer des individus, notamment en termes de justification théorique. Par contre, il faut le faire si l’analyse des résidus indique un véritable problème avec un ou des individus (inversion de valeurs, valeurs aberrantes ou fausses, etc.).
2 Régression linéaire multiple
Très souvent, une seule variable indépendante ne suffit pas à epliquer correctement la variable dépendante étudiée. Dans l’exemple précédent, nous n’expliquions que 79% de la variance de l’amplitude thermique. L’intégration d’autres variables permettraient donc d’augmenter ce pourcentage. Lorsqu’un modèle de régression contient plus d’une variable explicative, on parlera de régression multiple.
Nous allons dans l’exemple qui suit tenter de modéliser la température moyenne annuelle dans les 33 stations en fonction de l’altitude, de la latitude et de la distance au littoral, 3 variables supposées être des facteurs pertinents.
Le modèle s’écrit donc de la façon suivante :
2.1 Exploration bivariée des variables en jeu
Commençons par obersver la distribution de la variable dépendante :
hist(df$temp, freq=FALSE, breaks = 15, col = "darkblue", border = "white",
main = "Température moyenne annuelle")
La distribution est dissymétrique à gauche, mais cela peut être dû à un échantillon trop petit et peu représentatif. Nous allons poursuivre l’analyse en observant une matrice de graphiques cartésiens (entre la température et chacune des variables exogènes candidates), ce que permet une fonction très pratique du package car, la fonction scatterplotMatrix() :
scatterplotMatrix(~temp+dist_litt+alt+Y,
data=df,smoother=loessLine,
main="Matrice de graphiques cartésiens")
On constate que la température semble liée à l’altitude et à la latitude, mais pas à la distance au littoral (la droite MCO est presque horizontale). On complète cette exploration bivariée par le calcul de la matrice de corrélations. Le package Hmisc (fonction rcorr(), qui fonctionne sur des matrices) le permet :
library(Hmisc)
rcorr(as.matrix(df[,c("temp","dist_litt","alt","Y")])) temp dist_litt alt Y
temp 1.00 -0.05 -0.42 -0.68
dist_litt -0.05 1.00 0.55 -0.35
alt -0.42 0.55 1.00 -0.32
Y -0.68 -0.35 -0.32 1.00
n= 33
P
temp dist_litt alt Y
temp 0.8026 0.0159 0.0000
dist_litt 0.8026 0.0010 0.0443
alt 0.0159 0.0010 0.0698
Y 0.0000 0.0443 0.0698
Cette matrice confirme que la distance au littoral n’est pas associée à la température moyenne, contrairement aux deux autres variables. On peut aussi représenter visuellement ces corrélations bivariées, notamment grâce au package corrplot (voir http://www.sthda.com/french/wiki/matrice-de-correlation-avec-r-analyse-et-visualisation).
mcor <- cor(df[,c("temp","dist_litt","alt","Y")])
corrplot::corrplot(mcor, type="upper", order="hclust", tl.col="black")
2.2 Multicolinéarité
L’un des points les plus importants en régression multiple est de respecter l’indépendance entre les variables exogènes (absence de colinéarité), le risque étant sinon de biaiser les estimations des erreurs types des coefficients de régression. La VIF (variance inflation factor) est une bonne méthode pour détecter ces problèmes de colinéarité. La fonction vif() du package car permet de la calculer. Il faut d’abord estimer le modèle complet, puis appliquer la fonction à l’objet issu de lm():
vif(lm(data = df, formula = temp ~ dist_litt + alt + Y))dist_litt alt Y
1.500893 1.464131 1.172594 Les VIF sont toutes < 2, ce qui est tout à fait raisonnable. Les variables exogènes sont très peu corrélées entre elles.
2.3 Estimation du modèle
On peut donc finalement passer au modèle complet :
lmMult <- lm(data = df, formula = temp ~ dist_litt + alt + Y)
summary(lmMult)
Call:
lm(formula = temp ~ dist_litt + alt + Y, data = df)
Residuals:
Min 1Q Median 3Q Max
-0.7950 -0.1923 -0.1002 0.1401 1.1057
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.521e+01 2.202e+00 20.532 < 2e-16 ***
dist_litt 2.450e-04 5.168e-04 0.474 0.639
alt -5.975e-03 5.442e-04 -10.979 7.61e-12 ***
Y -4.984e-06 3.253e-07 -15.320 1.94e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4088 on 29 degrees of freedom
Multiple R-squared: 0.9143, Adjusted R-squared: 0.9055
F-statistic: 103.2 on 3 and 29 DF, p-value: 1.419e-15
Avec le package gtsummary :
gtsummary::tbl_regression(lmMult) %>% add_vif()| Characteristic | Beta | 95% CI | p-value | VIF |
|---|---|---|---|---|
| dist_litt | 0.00 | 0.00, 0.00 | 0.6 | 1.5 |
| alt | -0.01 | -0.01, 0.00 | <0.001 | 1.5 |
| Y | 0.00 | 0.00, 0.00 | <0.001 | 1.2 |
| Abbreviations: CI = Confidence Interval, VIF = Variance Inflation Factor | ||||
Ou présenté de façon graphique avec le package GGally :
GGally::ggcoef_model(lmMult)
Ou encore avec le package ggstatsplot, qui précise les coefficients estimés, directement sur le graphique :
ggstatsplot::ggcoefstats(lmMult, exclude.intercept = TRUE)
Nous arrivons avec ce modèle à expliquer plus de 90% de la variation de la température moyenne, ce qui est très satisfaisant. On constate que la latitude et l’altitude sont très fortement associés à la température (p-values très proches de 0). Quand l’altitude augmente de 100 m, la température diminue de 0,59 °C en moyenne. De même, quand on monte de 100 km en latitude, la température diminue de 0,49 °C en moyenne. La distance au littoral en rechanche n’apparait pas du tout significative (p-value = 0,64).
2.3.1 Effets partiels
Un autre élément intéressant consiste à observer les effets partiels, c’est-à-dire l’effet d’un des prédicteurs sur Y en tenant compte des effets des autres prédicteurs. On appelle cela le graphique des résidus partiels (à ne pas confondre avec le graphique de régression partielle). Le graphique des résidus partiels consiste à placer en X le prédicteur d’intérêt et en Y les résidus du modèle auxquels on ajoute la valeur du prédicteur multiplié par son coefficient de régression.
Le package car a une fonction très utile pour cela :
car::crPlots(model = lmMult, terms = ~.)
2.3.2 Régression pas à pas
La fonction step permet de réaliser une régression pas à pas descendante (ou ascendante) ainsi :
- Régression pas à pas descendante :
step(lmMult, direction = "backward")Start: AIC=-55.3
temp ~ dist_litt + alt + Y
Df Sum of Sq RSS AIC
- dist_litt 1 0.038 4.885 -57.042
<none> 4.847 -55.297
- alt 1 20.147 24.995 -3.169
- Y 1 39.232 44.080 15.553
Step: AIC=-57.04
temp ~ alt + Y
Df Sum of Sq RSS AIC
<none> 4.885 -57.042
- alt 1 25.377 30.262 1.142
- Y 1 41.878 46.763 15.503
Call:
lm(formula = temp ~ alt + Y, data = df)
Coefficients:
(Intercept) alt Y
4.548e+01 -5.848e-03 -5.019e-06
Sans surprise, on constate que ne pas intégrer la variable distance au littoral n’a que peu d’impact sur la performance du modèle (cela confirme l’observation des effets partiels). Selon le principe de parcimonie, il convient de ne garder dans le modèle final que les deux variables significatives.
2.3.3 Export des résultats
Il est souvent utile de pouvoir exporter la table des estimations dans un format .doc pour pouvoir être incorporée directement dans un rapport, un mémoire ou un article. Vous pouvez bien entendu le faire manuellement, mais il existe aussi des fonctions très utiles pour le faire automatiquement, comme la fonction tab_model() du package sjPlot :
tabReg <- sjPlot::tab_model(lmMult,show.stat = TRUE, show.se = TRUE,
show.ci = FALSE, file = "tabReg.doc")