[MOD1] : Régression linéaire

GEO UNIV’R Tunisie 2024

Auteur·rice

Claude Grasland, Malika Madelin

Date de publication

2024-05-13

Objectif

On se propose dans ce TD de modéliser la relation entre latitude (X) et température moyenne (Y) en Tunisie

Contrairement à la corrélation linéaire qui fait jouer un rôle symétrique au variables X et Y (\(r_{XY} = r_{YX}\)), la régression linéaire va introduire une dissymétrie en donnant à chacune des variables X et Y un rôle différent et en introduisant une hypothèse de causalité ou de dépendance :

  • la variable Y est la variable dépendante, c’est-à-dire celle que l’on veut expliquer ou prédire.
  • la variable X est la variable indépendante, c’est-à-dire la variable explicative ou du moins celle qui permet de prédire lesvaleurs de Y.

Dans notre exemple, il semble logique de considérer que les températures moyennes (Y) sont une conséquences de la position en latitude (X). Nous cherchons donc un modèle de la forme \(Y = f(X)\) dans lequel la fonction \(f\) peut prendre différentes formes.

Nous commencerons par le cas le plus simple d’une relation linéaire prenant la forme \(Y = a.X+b\) On commencera donc par utiliser un modèle de régression linéaire simple.

1. PREPARATION DES DONNEES

On utilise le tableau habituel de 21 stations de Tunisie

Télécharger le jeux de données

1.1 Importation des données

Tableau de données
code nom reg_code reg_nom lat lon alt tmin tmax tmoy prec vent rosee
BJ Beja NO Nord-ouest 36.73333 9.183333 159.00 11.45894 24.62033 18.55692 590.000 6.436291 10.881502
BZ Bizerte NE Nord-est 37.24545 9.791453 6.09 13.36735 23.96475 18.73467 639.000 13.080826 13.507263
GB Gabes SE Sud-est 33.87692 10.103333 7.92 15.62121 25.44585 20.64563 192.000 11.564158 12.396331
GF Gafsa SO Sud-ouest 34.42202 8.822503 323.08 13.73596 26.88865 20.42736 163.000 12.759514 7.894512
JE Jendouba NO Nord-ouest 36.48333 8.800000 144.00 12.03502 25.75739 19.17381 458.000 8.187385 11.549468
KR Kairouan CO Centre-ouest 35.66667 10.100000 68.00 15.05704 27.37991 20.88537 303.000 5.826256 10.894726
KS Kasserine CO Centre-ouest 35.16667 8.833333 707.00 11.47720 24.13852 18.46683 340.000 12.853035 8.366253
KB Kebili SO Sud-ouest 33.70000 8.966667 46.00 15.34825 28.66706 22.71979 89.000 14.061362 10.373775
KF Le Kef NO Nord-ouest 36.13333 8.700000 518.00 10.30506 23.19958 17.24170 528.000 8.246948 8.059934
MH Mahdia CE Centre-est 35.50000 11.066667 12.00 16.50674 23.29379 20.16491 290.000 10.854368 14.624500
ME Medenine SE Sud-est 33.35000 10.483333 117.00 16.17117 27.56439 22.57362 159.000 8.743863 11.404804
MS Monastir CE Centre-est 35.66667 10.750000 2.00 15.64783 24.26410 19.76152 322.062 13.888661 13.587094
NB Nabeul NE Nord-est 36.46667 10.700000 78.00 15.70167 23.43309 19.81197 450.000 12.190401 14.033335
SF Sfax CE Centre-est 34.71795 10.690972 25.90 14.88812 25.24082 20.05744 221.000 11.408961 12.427977
SZ Sidi Bou Zid CO Centre-ouest 35.00000 9.483333 355.00 13.03144 26.20685 20.25616 280.000 8.283933 9.461638
SL Siliana NO Nord-ouest 36.06667 9.366667 445.00 11.36109 24.50220 19.02255 389.000 9.242761 9.455365
SS Sousse CE Centre-est 35.70000 10.600000 5.00 15.70000 24.40000 19.90000 310.000 12.500000 13.600000
TA Tataouine SE Sud-est 32.91667 10.450000 215.00 15.90320 27.20100 21.41218 110.000 10.298155 9.182343
TO Tozeur SO Sud-ouest 33.93972 8.110556 87.47 16.83189 28.64269 22.63410 97.000 14.732460 8.697420
TU Tunis NE Nord-est 36.85103 10.227217 6.70 14.26899 24.61744 19.14238 453.000 13.405618 12.626422
ZA Zaghouan NE Nord-est 36.43333 10.083333 156.00 12.43361 25.06573 18.57657 501.000 10.331711 11.907365

On charge ensuite le fichier des métadonnées:

# Importe les métadonnées
meta<-read_xlsx(path = "data/TUN-CLIMAT/tun_climat.xlsx",
              sheet = "meta")
kable(meta, caption = "Tableau de métadonnées")
Tableau de métadonnées
variable def
code code de la station
nom nom de la station
reg_code code de la région
reg_nom nom de la region
gouv_nom nom du gouvernorat
lat Latitude
lon Longitude
alt Altitude (en mètres)
tmin Moyenne de Tmin
tmax Moyenne de Tmax
tmoy Moyenne de Tmoy
prec Moyenne de précipitations totales
vent Moyenne de vent km/h
rose Moyenne de point de rosée

1.2 Sélection des variables

On décide de garder les deux variables et de les renommer X et Y conformément à nos hypothèses.

  • X : Latitude en degrés
  • Y : Température moyenne en degrés Celsius

On procède donc à l’extraction de ces variables en y ajoutant le nom et le code iso des stations.

# Création des variables X et X
don$X<-don$lat
don$Y<-don$tmoy

# Sélection des colonnes
tab<-don[,c("code","nom","X","Y")]

1.3 Astuce : stockage des textes d’habillage

On prépare un ensemble de textes que l’on pourra utiliser pour l’habillage de nos graphiques. Cela évitera de devoir ensuite les retaper à chaque fois.

On décide ici que les textes seront en français :

nomX <- "Latitude "
nomY <- "Température moyenne en degrés" 
titre <- "Le climat de Tunisie"
note <- "Source : Salem Dahech, 2024"

2. ANALYSE DE LA VARIABLE Y

2.1 Calculer les paramètres principaux

summary(tab$Y)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  17.24   19.02   19.90   20.01   20.65   22.72 
sd(tab$Y)
[1] 1.459039
  • Commentaire : Les températures moyennes vont de 17.24 à 22.72 degrés avec une moyenne de 20.01 degrés et un écart type de 1.46.

2.2 Faire un histogramme

  • Histogramme rapide
hist(tab$Y)

  • Histogramme amélioré
hist(tab$Y, 
     xlab=nomY,
     breaks=quantile(tab$Y, c(0,0.25,0.5,0.75,1)),
     xlim=c(17,23),
     main = titre,
     sub = note,
     col = "lightyellow")
lines(density(tab$Y),col="red")
rug(tab$Y)

  • Commentaire : La distribution est globalement symétrique et unimodale malgré un petit mode secondaire

2.3 Tester la normalité

Pour savoir si une distribution est gaussienne (normale) on peut utiliser un test statistique (test de Shapiro-Wilks) à l’aide de la fonction shapiro.test() et tracer un graphique d’écart à la loi gaussienne à l’aide des fonctions qqnorm() et qqline() :

En statistiques, un Q-Q (quantile-quantile) plot est une méthode graphique pour comparer deux distributions de probabilité en affichant leur quantiles contre quantiles. Un point (x,y) du graphique représente un quantile de la seconde distribution (axe y) contre le même quantile de la première distribution (axe x). Ainsi la droite est une courbe paramétrique dont le paramètre est le nombre d’intervalle des quantiles.

Interprétation d’un qqplot

  • Normal qqplot: La distribution normale est symétrique, donc aucun biais (skew) et la moyenne est égale à la médiane.

  • Right skewed qqplot : Right-skew aussi appelé positive skew signifie que la distribution comporte des valeurs exceptionnelles à droite et que la moyenne est supérieure à la médiane.

  • Left skewed qqplot: Left-skew aussi appelé negative skew signifie que la distribution comporte des valeurs exceptionnelles à gauche et que la moyenne est inférieure à la médiane

  • Light tailed qqplot: Cela veut dire que comparé à la distribution normale il y a un peu plus de données dans les extrémités que dans le centre de la distribution.

  • Heavy tailed qqplot: Cela veut dire que comparé à la distribution normale il y a un beaucoup plus de données dans les extrémités que dans le centre de la distribution.

  • Biomodel qqplot: illustre une distribution bimodale comportant deux zones de concentration avec donc deux pics sur l’histogramme.

Source : Zach Bogart & Joyce Robbins, 2019, EDAV-Info

# Graphique 
qqnorm(tab$Y)
qqline(tab$Y, col = "red")

# test
shapiro.test(tab$Y)

    Shapiro-Wilk normality test

data:  tab$Y
W = 0.95026, p-value = 0.3446
  • Commentaire : Le graphique montre que la distribution suit approximativement une loi gaussienne, ce qui est confirmé par le test de Shapiro-Wilks (p > 0.344)

2.4 Examiner la présence de valeurs exceptionnelles

La solution la plus courante est d’utiliser une boxplot :

boxplot(tab$X, 
        horizontal = T,
        xlab = nomX,
        main = titre,
        sub = note)

  • Commentaire : La boxplot ne montre la présence d’aucune valeur exceptionnelle.

3. CORRELATION

3.1 Visualiser la relation entre X et Y

  • Graphique rapide
plot(tab$X,tab$Y)

  • Graphique amélioré
plot(don$X,don$Y,
     main = titre,   # titre
     cex.main = 1,      # police du titre
     cex.sub = 0.6,     # police du sous-titre
     xlab = nomX,    # nom de l'axe X
     xlim = c(32.7,37.5),   # intervalle de l'axe X
     ylab = nomY,    # nom de l'axe Y
    ylim = c(17,24),    # intervalle de l'axe Y
     cex.axis = 0.8,    # police des gradations d'axes
     cex.lab = 0.8,     # police des noms d'axes
     cex = 0.6,         # taille des symboles
     col = "blue")       # couleur des symboles


# Ajout d'une ligne horizontale  correspondant à la moyenne de Y
abline(h=mean(don$Y),col="red",lwd = 1, lty = 2)
# Ajout d'une ligne verticlae  correspondant à la moyenne de X
abline(v=mean(don$X),col="red",lwd = 1, lty = 2)

text(x = don$X,
     y = don$Y,
     label = don$code,
     cex = 0.7,
     pos=3,
     col = "blue")

  • Commentaire : : La relation semble négative et linéaire

3.2 Tester la significativité de la relation entre X et Y

Coefficient de Pearson

cor.test(tab$X,tab$Y)

    Pearson's product-moment correlation

data:  tab$X and tab$Y
t = -5.7451, df = 19, p-value = 1.549e-05
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.9140080 -0.5562694
sample estimates:
       cor 
-0.7966527 
cor(tab$X,tab$Y)**2
[1] 0.6346554
  • Commentaire : Selon le test du coefficient de Pearson, la relation est très significative (p < 0.001) et le pouvoir explicatif de X par rapport à Y mesuré par la coefficient de détermination (\(r_{XY}^2\)) sera élevé (63%).

Coefficien de Spearman

cor.test(tab$X,tab$Y, method = "spearman")

    Spearman's rank correlation rho

data:  tab$X and tab$Y
S = 2738, p-value = 4.83e-05
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.7779221
  • Commentaire : Le coefficient de corrélation de Spearman (-0.77) est sensiblement égal à celui celui de Pearson (+0.80). Ceci est en général bon signe et confirme que la distribution est sans doute linéaire.

4. REGRESSION LINEAIRE

4.1 Calculer l’équation de la droite Y = aX+B

modreglin <- lm(tab$Y~tab$X)
summary(modreglin)

Call:
lm(formula = tab$Y ~ tab$X)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.02273 -0.43088  0.06247  0.54638  1.32667 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  52.9170     5.7316   9.232 1.87e-08 ***
tab$X        -0.9313     0.1621  -5.745 1.55e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9048 on 19 degrees of freedom
Multiple R-squared:  0.6347,    Adjusted R-squared:  0.6154 
F-statistic: 33.01 on 1 and 19 DF,  p-value: 1.549e-05
  • Commentaire : L’équation de la droite est donc \(Y =-0.931\times X + 52.92\). Le coefficient de pente de la droite indique que les températures diminuent de 0.93 degrés chaque fois que la latitude augmente de 1. Comme 1 degré vaut environ 100 km, cela signifie que la température baisse du sud vers le nord d’environ 0.01 degré par km.

4.2 Visualiser la droite

plot(don$X,don$Y,
     main = titre,   # titre
     cex.main = 1,      # police du titre
     cex.sub = 0.6,     # police du sous-titre
     xlab = nomX,    # nom de l'axe X
     xlim = c(32.7,37.5),   # intervalle de l'axe X
     ylab = nomY,    # nom de l'axe Y
    ylim = c(17,24),    # intervalle de l'axe Y
     cex.axis = 0.8,    # police des gradations d'axes
     cex.lab = 0.8,     # police des noms d'axes
     cex = 0.6,         # taille des symboles
     col = "blue")       # couleur des symboles


# Ajout d'une ligne horizontale  correspondant à la moyenne de Y
abline(h=mean(don$Y),col="red",lwd = 1, lty = 2)
# Ajout d'une ligne verticlae  correspondant à la moyenne de X
abline(v=mean(don$X),col="red",lwd = 1, lty = 2)

text(x = don$X,
     y = don$Y,
     label = don$code,
     cex = 0.7,
     pos=3,
     col = "blue")

abline(modreglin, col ="black", lwd =2)

  • Commentaire: La droite s’ajuste assez bien au nuage de points mais certains points en sont assez éloignés.

4.3 Calculer les valeurs estimées et les résidus

# Extraction des valeurs estimées et résiduelles
tab$Yest <- modreglin$fitted.values
tab$Yres <- modreglin$residuals

# Affichage du tableau trié
kable(tab[order(tab$Yres),])
code nom X Y Yest Yres
KF Le Kef 36.13333 17.24170 19.26442 -2.0227251
KS Kasserine 35.16667 18.46683 20.16472 -1.6978858
TA Tataouine 32.91667 21.41218 22.26024 -0.8480634
GB Gabes 33.87692 20.64563 21.36592 -0.7202877
SF Sfax 34.71795 20.05744 20.58263 -0.5251847
GF Gafsa 34.42202 20.42736 20.85824 -0.4308751
ZA Zaghouan 36.43333 18.57657 18.98502 -0.4084519
SL Siliana 36.06667 19.02255 19.32651 -0.3039556
BJ Beja 36.73333 18.55692 18.70561 -0.1486901
SZ Sidi Bou Zid 35.00000 20.25616 20.31994 -0.0637850
MS Monastir 35.66667 19.76152 19.69905 0.0624704
SS Sousse 35.70000 19.90000 19.66800 0.2319971
JE Jendouba 36.48333 19.17381 18.93845 0.2353609
MH Mahdia 35.50000 20.16491 19.85427 0.3106395
BZ Bizerte 37.24545 18.73467 18.22866 0.5060068
TU Tunis 36.85103 19.14238 18.59600 0.5463832
ME Medenine 33.35000 22.57362 21.85666 0.7169607
NB Nabeul 36.46667 19.81197 18.95397 0.8579988
KR Kairouan 35.66667 20.88537 19.69905 1.1863196
KB Kebili 33.70000 22.71979 21.53069 1.1890956
TO Tozeur 33.93972 22.63410 21.30743 1.3266717

Commentaire : Le tableau permet de repérer les stations qui s’éloignent le plus de la droite en raison d’une surestimation ou d’une sous-estimation de leurs température par la latitude.

  • Les résidus négatifs correspondent à des stations dont la température est moins chaude que ce que laisserait prévoir leur latitude. C’est par exemple le cas de la station d’El Kef dont la latitude (36.13) laissait prévoir une températude de 19.2 degrés mais qui en pratique a une température de 17.2° soit un résidu de 2 degrés de moins que prévu.

  • Les résidus positifs correspondent à des stations dont la température est plus chaude que ce que laisserait prévoir leur latitude. C’est par exemple le cas de la station de Tozeur dont la latitude (33.94) laissait prévoir une températude de 21.3 degrés mais qui en pratique a une température de 22.6 degrés soit un résidu de 1.3 degrés de plus que prévu.

4.4 Sauvegarder les résultats du modèle

On peut si on le souhaite sauvegarder les résultats au format .csv

write.table(x = tab,
            file = "result.csv",
            row.names = FALSE)

5. EVALUATION DU MODELE

Avant de tirer des conclusions hâtives sur les résidus, il est préférable de vérifier si les hypothèses fondamentales du modèle de régression ont bien été respectées. On va utiliser pour cela quatre graphiques de bases fournis par R et des tests présents dans le package car (acronyme de “Companion for Applied Regression”).

5.1 Autocorrélation des résidus

library(car)
plot(modreglin,
     which = 1,
     labels.id = tab$nom,
     col="red")

durbinWatsonTest(modreglin)
 lag Autocorrelation D-W Statistic p-value
   1      -0.5133805      3.014614   0.014
 Alternative hypothesis: rho != 0
  • Commentaire : le graphique permet de voir que les résidus sont dans l’ensemble indépendants des valeurs estimées de Y, ce qui signifie que les points sont bien répartis autour de la droite On peut s’en assurer à l’aide du test de Durbin Watson qui pose l’hypothèse H0 : Il existe une autocorrélation des résidus. Cette hypothèse peut être rejetée (p < 0.05) donc il n’existe pas d’autocorrélation susceptible de fausser les résultats

5.2 Normalité des résidus

plot(modreglin,
     which = 2,
     labels.id = tab$nom,
     col="red")

shapiro.test(tab$Yres)

    Shapiro-Wilk normality test

data:  tab$Yres
W = 0.95959, p-value = 0.5078
  • Commentaire : La normalité de la distribution des résidus est également une condition importante de validité du modèle de régression linéaire puisqu’elle permet de définir un intervalle de confiance des estimations en se servant de l’écart-type de ces résidus (e.g. + ou - 2 écarts-type pour un intervalle de confiance à 95%). Au vu du diagramme QQ plot on voit que la condition de normalité des résidus semble bien vérifiée, ce que confirme le test de Shapiro (p > 0.05)

5.3 Homogénéité des résidus

plot(modreglin,
     which = 3,
     labels.id = tab$nom,
     col="red")

ncvTest(modreglin)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 0.3034993, Df = 1, p = 0.5817
  • Commentaire : En liaison avec ce qui précède, l’analyse de l’homogénéité des résidus permet de vérifier si la variance des résidus est constante et donc si l’intervalle de confiance sera le même pour l’ensemble des valeurs estimées. Ici, c’est à peu près le cas même si les résidus varient un peu en fonction de Y. On peut vérifier l’homogénéité en appliquant le test de Breush-Pagan qui examine l’hypothèse “H0 : la distribution des résidus est homogène”. Dans notre exemple H0 ne peut pas être rejetée (p < 0.001) ce qui signifie que l’hypothèse d’homogénéité des résdus est vérifiée

5.4 Absence de valeurs exceptionnellement influentes

plot(modreglin,
     which = 4,
     labels.id = tab$nom,
     col="red")

outlierTest(modreglin, labels = tab$nom)
No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:
        rstudent unadjusted p-value Bonferroni p
Le Kef -2.660486           0.015932      0.33457
  • Commentaire : Le dernier test consiste à vérifier si la relation observée est bien le résultat d’un ensemble d’observations indépendantes et non pas l’effet de la présence d’une ou deux valeurs exceptionnelles. Plusieurs tests sont ici possibles qui visent au même objectif : déterminer à quel point le retrait d’unevaleur unique modifie le résultat de l’analyse, c’est à dire le coefficient de détermination \(r_{XY}^2\) et les paramètres \(a\) et \(b\) de l’équation \(Y=aX+b\). Le graphique proposé par R utilise la distance de Cook pour mettre en valeur l’influence potentielle des valeurs exceptionnelles et on y retrouve les trois stations de Kebeli, El Kef et Tataouine. On peut arriver à un résultat similaire en utilisant le test de Bonferroni qui signale le caractère influent de la station d’El Kef (p < 0.05)

5.5 Tous les tests d’un coup

Une fois que l’on a bien compris les tests précédents, on peut afficher les quatre graphiques correspondant en une seule commande :

par(mfrow=c(2,2))
plot(modreglin,
     which = c(1,2,3,4),
     labels.id = tab$nom,
     col="red")

CONCLUSION

Le modèle apparaît finalement valide à tous points de vue, la seule réserve étant le rôle très influent de la station d’El Kef. Notre modèle n’explique cependant que 64% de la variance des températures ce qui laisse supposer que d’autres facteurs sont à l’oeuvre pour expliquer les 36% restants. On peut penser à l’altitude, la distance à la mer, la situation d’abri face au vent, …

Pour aboutir à un modèle plus complet, on utilisera alors la régression multiple.