Pour faire suite aux billets sur la régression linéaire simple et la régression linéaire multiple par la méthode des moindres carrés, nous allons nous intéresser à présent au problème de multicolinéarité.

Il me semble qu’un des premiers réflexes à avoir face aux données, c’est de passer pas mal de temps à les regarder, s’intéresser aux corrélations, aux types de relations entre les différentes variables, tenter de déceler des points qui nous semblent étranges, etc.

Comme dans les autres billets liés à la régression sur ce blog, l’objectif n’étant pas de donner une recette de cuisine (la méthode idéale n’existe pas), nous allons passer un peu rapidement sur cette étape.

Considérons les données suivantes :

# Description des variables :
# KW : Nombre de KWH consommés pendant le mois de janvier
# SURFACE : Surface du logement en m^2
# PERS : Nombre de  personnes habitant le logement
# PAVILLON :  Pavillon (1); appartement (0)
# AGE : Age du logement
# VOL : Volume interieur du logement en m^3
# SBAINS :  Nombre  de salles de bains

kw <- c(4805,3783,2689,5683,3750,2684,1478,1685,1980,1075,2423,4253,
        1754,1873,3487,2954,4762,3076)
surface <- c(130,123,98,178,134,100,78,100,95,78,110,130,73,87,152,
             128,180,124)
pers <- c(4,4,3,6,4,4,3,4,3,4,5,4,2,4,5,5,7,4)
pavillon <- c(1,1,0,1,1,0,0,0,0,0,1,1,0,1,1,1,1,0)
age <- c(65,5,18,77,5,34,7,10,8,5,12,25,56,2,12,20,27,22)
vol <- c(410,307,254,570,335,280,180,250,237,180,286,351,220,217,400,
         356,520,330)
sbains <- c(1,2,1,3,2,1,1,1,1,1,1,1,1,2,2,1,2,1)

df <- data.frame(kw, surface, pers, pavillon, age, vol, sbains)

Tout d'abord, on peut remarquer que la variable comptant le nombre de salle de bains est discrète. On pourrait vouloir l'utiliser en tant que variable catégorielle, avec des niveaux ordonnés. Cependant, le nombre de données à notre disposition est ici très faible, et il n'y a qu'une observation dans laquelle il y a trois salles de bains.

On peut construire une fonction très simple pour effectuer rapidement les graphiques de la consommation d'électricité en fonction de chaque variable. On y ajoute un lissage pour aider également à la visualisation d'une éventuelle tendance :

graphConso <- function(uneVariable){
  ggplot(data = df, aes_string(x = uneVariable, y = "kw")) + geom_point() + geom_smooth()
}
p1 <- graphConso("surface")
p2 <- graphConso("pers")
p3 <- graphConso("pavillon")
p4 <- graphConso("age")
p5 <- graphConso("vol")
p6 <- graphConso("sbains")

library(gridExtra)
grid.arrange(p1, p2, p3, p4, p5, p6, ncol = 2)

Consommation d'électricité en fonction des covariables

Il se dégage des relations qui ont l'air d'être plus ou moins linéaires entre la consommation et les autres covariables. La variable pavillon est binaire, il peut s'avérer intéressant de regarder comment évolue la consommation d'électricité en fonction des autres variables, et de la variable pavillon.

ggplot(data = df, aes(x =pers , y = kw, group = as.factor(pavillon), col = as.factor(pavillon))) +
  geom_point() + geom_smooth() +
  scale_colour_discrete("Type de logement", labels = c("Appartement", "Pavillon"))

Il peut également s'avérer utile de regarder les boxplots pour chaque variable, pour avoir une idée rapide des quartiles.

graphConso.bp <- function(uneVariable){
  ggplot(data = df, aes_string(x = "factor(0)", y = uneVariable)) +
    geom_boxplot() + xlab("")
}

p1.bp <- graphConso.bp("surface")
p2.bp <- graphConso.bp("pers")
p3.bp <- graphConso.bp("pavillon")
p4.bp <- graphConso.bp("age")
p5.bp <- graphConso.bp("vol")
p6.bp <- graphConso.bp("sbains")

grid.arrange(p1.bp, p2.bp, p3.bp, p4.bp, p5.bp, p6.bp, ncol = 2)

Bien évidemment, il faudrait aller plus loin dans les statistiques descriptives, mais le but de l'exercice est de mettre en avant les problèmes de multicolinéarité. Regardons ce que nous donne la régression par MCO de la consommation d'électricité sur les autres variables.

> summary(reg <- lm(kw ~ surface + pers + pavillon + age + vol + sbains, data = df))

Call:
lm(formula = kw ~ surface + pers + pavillon + age + vol + sbains, 
    data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-653.60 -152.55   26.47  231.58  485.82 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) -226.408   1344.210  -0.168    0.869  
surface       28.290     54.305   0.521    0.613  
pers        -456.262    229.229  -1.990    0.072 .
pavillon     595.984    276.274   2.157    0.054 .
age            8.407     27.067   0.311    0.762  
vol            4.521     20.086   0.225    0.826  
sbains       -78.456    219.732  -0.357    0.728  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 398.2 on 11 degrees of freedom
Multiple R-squared:  0.9392,	Adjusted R-squared:  0.9061 
F-statistic: 28.33 on 6 and 11 DF,  p-value: 4.483e-06

Au seuil de 5%, aucun des coefficients n'est statistiquement différent de zéro. Il faut se fixer un seuil de 10% pour rejeter le test de Student de nullité du coefficient associé à la variable pers et à pavillon. Si on regarde le coefficient de détermination, on s'aperçoit qu'il est élevé : 0.94. Le modèle paraît prédictif...

La variable surface n'a pas d'influence significative sur la consommation d'électricité. Pourtant, si on régresse la consommation sur la surface et une constante, voici ce que l'on observe :

> summary(lm(kw ~ surface, data = df))

Call:
lm(formula = kw ~ surface, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-834.63 -446.36  -75.79  361.87 1297.00 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1299.840    531.327  -2.446   0.0264 *  
surface        36.983      4.407   8.392 2.97e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 576.2 on 16 degrees of freedom
Multiple R-squared:  0.8149,	Adjusted R-squared:  0.8033 
F-statistic: 70.42 on 1 and 16 DF,  p-value: 2.966e-07

Si on effectue une régression linéaire simple de la consommation d'électricité sur chacune des autres variables de la data frame, on peut remarquer que dans tous les cas, le test de nullité du coefficient est rejeté à 5%. Une fois qu'on les combine en un seul modèle, les variables qui prises séparément ont un pouvoir explicatif sur la consommation d'électricité le perdent par la suite. Les estimations des écarts-types de chaque estimateur augmentent dans le cas où on intègre toutes les variables dans un seul modèle, comparativement au cas où on les place seules.

La multicolinéarité peut être à l'origine de ce phénomène. Il existe plusieurs moyens qui ont été abordés pendant le cours pour la détecter. Nous allons dans un premier temps regarder quelques méthodes très simples et assez arbitraires, puis nous continuerons avec une méthode l'étant un poil moins.

Analyse de la corrélation

On regarde l'intensité de la liaison entre les différentes variables, en prenant appui sur le mesure de corrélation :
\[\rho_{XY} = \frac{\mathbb{C}ov(X,Y)}{\sigma_X \sigma_Y},\]
avec \(\sigma_X\) et \(\sigma_Y\) les écarts-types de \(X\) et \(Y\) respectivement.

> cor(df)
                kw   surface      pers  pavillon       age       vol    sbains
kw       1.0000000 0.9026958 0.6229604 0.6787712 0.5668317 0.9298826 0.5920095
surface  0.9026958 1.0000000 0.8373894 0.6764150 0.3635577 0.9696253 0.6597002
pers     0.6229604 0.8373894 1.0000000 0.6333333 0.1524229 0.7975005 0.5749610
pavillon 0.6787712 0.6764150 0.6333333 1.0000000 0.1160580 0.6338246 0.5889844
age      0.5668317 0.3635577 0.1524229 0.1160580 1.0000000 0.5722021 0.2001865
vol      0.9298826 0.9696253 0.7975005 0.6338246 0.5722021 1.0000000 0.6356709
sbains   0.5920095 0.6597002 0.5749610 0.5889844 0.2001865 0.6356709 1.0000000

Première règle

Si la valeur absolue du coefficient de corrélation entre deux variables excède 0.8, on peut soupçonner la colinéarité.
Un moyen visuel plutôt sympathique m'a été proposé par @philandrin : utiliser la fonction corrplot du package corrplot.

library(corrplot)
library(RColorBrewer)
corrplot.mixed(cor(df), col = rev(brewer.pal(10, "Spectral")))

Selon cette règle, on peut soupçonner un problème de colinéarité entre les variables vol et surface, ainsi que pers et surface.

Deuxième règle (règle de Klein)

Si le carré du coefficient de corrélation est supérieur au R^2, on peut soupçonner de la colinéarité.

> cor(df)^2 > summary(reg)$r.squared
            kw surface  pers pavillon   age   vol sbains
kw        TRUE   FALSE FALSE    FALSE FALSE FALSE  FALSE
surface  FALSE    TRUE FALSE    FALSE FALSE  TRUE  FALSE
pers     FALSE   FALSE  TRUE    FALSE FALSE FALSE  FALSE
pavillon FALSE   FALSE FALSE     TRUE FALSE FALSE  FALSE
age      FALSE   FALSE FALSE    FALSE  TRUE FALSE  FALSE
vol      FALSE    TRUE FALSE    FALSE FALSE  TRUE  FALSE
sbains   FALSE   FALSE FALSE    FALSE FALSE FALSE   TRUE

Cette règle met en avant un éventuel problème de colinéarité entre les variables vol et surface.

Troisième règle

Le signe de la correlation brute et celui de la correlation expliquée devraient être les memes.

> cbind(coef = reg$coef, corr = cor(df)[,"kw"])[-1,]
                coef      corr
surface    28.290102 0.9026958
pers     -456.261605 0.6229604
pavillon  595.983723 0.6787712
age         8.407310 0.5668317
vol         4.520655 0.9298826
sbains    -78.455549 0.5920095

Il y a manifestement des contre-sens concernant les relations entre le nombre de personne et la consommation, ainsi que le nombre de salles de bain et la consommation.

Une règle un peu plus technique

Facteur d'inflation (VIF)

La variance du \(j^e\) estimateur est donnée par :
\[\mathbb{V}(\hat{\beta}_j) = \frac{\sigma^2}{n \mathbb{V}(\boldsymbol{x}_j)} \frac{1}{1-R^2_j},\]
où \(n\) est le nombre d'observation, \(\sigma^2\) la variance des erreurs, \(R^2_j\) le coefficient de détermination obtenu en régression la \(j^e\) variable sur tous les autres régresseurs.

Le VIF est défini par :
\[VIF_j = \frac{1}{1-R^2_j}.\]

On voit donc que plus \(R^2_j\) est proche de \(1\) (plus le \(j^e\) régresseur est colinéaire aux autres régresseurs), plus la valeur du VIF est élevée, ce qui se traduit par une augmentation de l'estimation de la variance du coefficient \(\hat{\beta}_j\). Comme \(0\leq R^2_j \leq 1\), la valeur minimale que peut prendre le VIF est \(1\).

Dans le cours, il a été mentionné qu'il y a présence de colinéarité lorsque \(VIF_j >4\), ou encore de forte colinéarité lorsque \(VIF_j > 10\).

Avec R, on peut obtenir les VIFs à la main :

> reg.surface <- lm(surface ~ .-kw, data = df)
> reg.pers <- lm(pers ~ .-kw, data = df)
> reg.pavillon <- lm(pavillon ~ .-kw, data = df)
> reg.age <- lm(age ~ .-kw, data = df)
> reg.vol <- lm(vol ~ .-kw, data = df)
> reg.sbains <- lm(sbains ~ .-kw, data = df)
> 
> summary(reg.surface)$r.squared
[1] 0.9968553
> summary(reg.pers)$r.squared
[1] 0.8659122
> summary(reg.pavillon)$r.squared
[1] 0.5326795
> summary(reg.age)$r.squared
[1] 0.9737686
> summary(reg.vol)$r.squared
[1] 0.9980167
> summary(reg.sbains)$r.squared
[1] 0.4769785
> 
> c(surface = 1/(1-summary(reg.surface)$r.squared),
+   pers = 1/(1-summary(reg.pers)$r.squared),
+   pavillon = 1/(1-summary(reg.pavillon)$r.squared),
+   age = 1/(1-summary(reg.age)$r.squared),
+   vol = 1/(1-summary(reg.vol)$r.squared),
+   sbains = 1/(1-summary(reg.sbains)$r.squared))
   surface       pers   pavillon        age        vol     sbains 
317.999609   7.457799   2.139859  38.122307 504.210822   1.911967

Ou alors, on peut, de manière plus rapide, utiliser la fonction vif du package car :

> library(car)
> vif(reg)
   surface       pers   pavillon        age        vol     sbains 
317.999609   7.457799   2.139859  38.122307 504.210822   1.911967 

On note une forte colinéarité entre surface et les autres régresseurs, entre age et les autres régresseurs et entre volume et les autres régresseurs.

La racine carrée de \(VIF_j\) indique de combien de fois l'écart-type est modifié par rapport à la situation dans laquelle la variable \(\boldsymbol{x}_j\) serait décorrélée aux autres régresseurs.

> sqrt(vif(reg))
  surface      pers  pavillon       age       vol    sbains 
17.832544  2.730897  1.462826  6.174326 22.454639  1.382739 

Par exemple, l'écart-type du coefficient de la variable surface est \(17.83\) fois plus grand qu'il ne le serait si surface n'était pas corrélée aux autres variables.

Conclusion

Il n'y a pas de règle, comme je le disais au début de ce billet. Ce que j'aurais tendance à dire, c'est qu'il est préférable de réfléchir avant de faire les tests, et d'essayer d'imaginer ce que pourraient être les résultats avant de procéder aux tests. Dans notre exemple, il semble assez évident que des variables comme le volume et la surface vont être très liées pour des appartements ou des pavillons. Cependant, j'ai aussi tendance à me méfier des évidences... Selon nos croyances, on peut partir avec des hypothèses fortes nous empêchant de voir la partie submergée de l'iceberg...
Mais alors, que faire ici, dans cet exemple, pour continuer dans la recherche d'un modèle de régression ? J'aurais envie d'écarter les variables pour lesquelles le VIF est très élevé, et de pousser un peu plus loin pour la sélection de modèle.
Ce sera le sujet du prochain billet.

One thought on “[L3 Eco-Gestion] Régression linéaire avec R : problèmes de multicolinéarité

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *

Time limit is exhausted. Please reload CAPTCHA.