Après avoir présenté rapidement la régression linéaire multiple avec R, et parlé un peu des problèmes de multicolinéarité, on va se pencher sur différentes techniques qu’il est possible d’employer pour sélectionner un modèle. Bien sûr, il en existe beaucoup d’autres. Le but est ici de donner un rapide aperçu.

Dans un premier temps, on va s’attacher à sélectionner les variables qu’il faut intégrer à notre modèle, et celles qu’il n’est pas nécessaire de conserver.

Pour sélectionner le « meilleur modèle » (si tant est qu’il existe), on va s’appuyer sur une mesure qui permet de comparer les modèles entre eux (par exemple, le coefficient de détermination ajusté \(R^2_a\), le critère d’information bayésien (SIC ou BIC), ou encore le critère d’information d’Akaike (AIC)). Il existe deux méthodes souvent employées pour effectuer de la sélection :

  • méthode exhaustive ;
  • méthode pas-à-pas.

Dans la méthode exhaustive, il s’agit de régresser tous les modèles qu’il est possible de construire à partir des variables, puis de les comparer à l’aide d’une mesure par la suite. Lorsqu’on n’a pas beaucoup de régresseurs, cette technique peut être envisageable. Prenons un exemple : si on a 3 régresseurs (a, b et c), les différents modèles que l’on peut estimer sont les suivants :

lm(y ~ 1 + a)
lm(y ~ 1 + b)
lm(y ~ 1 + c)
lm(y ~ 1 + a + b)
lm(y ~ 1 + a + c)
lm(y ~ 1 + b + c)
lm(y ~ 1 + a + b + c)

Cela fait donc 7 estimations à réaliser.
Avec \(n\) régresseurs, il faut estimer \(\sum_{i=1}^{n} C_n^i\) modèles. Ainsi, rien qu’avec 5 variables explicatives, on passe à 31 régressions, et avec 6 variables explicatives, à 63.
Dans la méthode pas-à-pas, trois manières de procéder sont employées :

  • forward : on part du modèle avec uniquement une constante, et on ajoute les variables unes à unes jusqu’à ce que l’ajout d’une variable supplémentaire se solde par un modèle jugé moins bon en fonction du critère de comparaison sélectionné ;
  • backward : on part du modèle avec toutes les variables auxquelles on a pensé, en retirant à chaque étape une seule variable jusqu’à ce que la comparaison des modèles indique qu’il est préférable de ne plus retirer de variable ;
  • stepwise : un mélange des méthodes forward et backward, basée sur le F-partiel. On vérifie que l’ajout d’une variable ne provoque pas la suppression d’une variable déjà introduite

On reprend les mêmes données que précédemment (Régression linéaire avec R : sélection de modèle).

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)

Nous allons procéder en "oubliant" qu'il y a des problèmes de multicolinéarité (et nous nous apercevrons que les résultats obtenus sont un peu étrange, et suggèrent que la personne qui créé le modèle doit rester vigilante avec les techniques automatisées).

1. Sélection en backward en se fiant au critère d'Akaike (AIC)

On va donc partir du modèle où la consommation d'électricité est modélisée à l'aide de toutes les variables de notre data frame, et se fier au critère AIC pour décider de retirer ou non une variable, à chaque itération.

Le critère AIC est défini par :
\[AIC = n \times log\left(\frac{SCR}{n}\right) + 2(l+1),\]
avec \(n\) le nombre d'observations, \(SCR\) la somme des carrés des résidus du modèle estimé et \(l\) le nombre de variables explicatives.

Intuitivement, si le modèle s'ajuste mieux aux données, la somme des carrés des résidus va diminuer et donc le terme \(n \times log\left(\frac{SCR}{n}\right)\) va également diminuer. Le critère AIC devrait donc potentiellement baisser, en fonction de la valeur du terme restant \(2(l+1)\), qui vient pénaliser les modèles avec plus de variables. En effet, la somme des carrés des résidus augmente ou reste stable avec l'ajout de variables dans le modèle. Aussi, il est nécessaire de pénaliser l'ajout de variables supplémentaires.

On peut aisément faire une fonction pour obtenir l'AIC avec R :

aic <- function(reg){
  n <- nrow(reg$model)
  SCR <- sum(reg$residuals^2)
  p <- ncol(reg$model)-1
  resul <- n * log(SCR/n) + 2*(p+1)
  return(resul)
}

L'opération consiste donc à partir du modèle complet, de retirer une variable et de voir quel retrait entraîne la plus forte baisse de l'AIC. Si le retrait d'une variable ne se solde pas par une diminution de l'AIC, on s'arrête. Sinon, on recommence le processus de retrait.

Pour faciliter l'opération, voici une fonction qui retourne la valeur de l'AIC du modèle avec le plus de variable, ainsi que les valeurs du critère AIC pour les modèles auxquels on retire une variable :

AIC.retirer <- function(y, regresseurs){
  laFormule <- paste(y, "~", paste(regresseurs, collapse = "+"), sep = "")
  reg.all <- lm(laFormule)

  obtenirAICReg <- function(y, regresseurs, uneVariable){
    # On retire uneVariable de la liste des regresseurs
    regresseurs <- regresseurs[-which(regresseurs%in%uneVariable)]
    laFormule <- paste(y, "~", paste(regresseurs, collapse = "+"), sep = "")
    reg.tmp <- lm(laFormule, data = df)
    return(c(variable = uneVariable, AIC = aic(reg.tmp)))
  }
  
  res1 <- c(AIC = aic(reg.all))
  res2 <- do.call("rbind", lapply(regresseurs, function(x) obtenirAICReg(y = "kw", regresseurs = regresseurs, uneVariable = x)))
  
  if(any(res2[,"AIC"] < res1)){
    # On peut retirer une variable
    varARetirer <- as.character(res2[which(res2[,"AIC"] == min(res2[,"AIC"])), "variable"])
    cat("La variable faisant le plus baisser l'AIC est :", varARetirer, sep = "\n")
  }else{
    # Tous les AIC obtenus en retirant une variable sont superieurs ou egaux
    # a l'AIC obtenu avec le modele ou figurent toutes les variables
    cat("Le retrait d'une variable fait augmenter l'AIC, on conserve le modele", laFormule ,sep = "\n")
  }
  return(list(all = res1, moinsUn = res2))
}

Voici les différentes étapes :

> # Premiere etape : modele complet
> AIC.retirer(y = "kw", regresseurs = c("surface", "pers", "pavillon", "age", "vol", "sbains"))
La variable faisant le plus baisser l'AIC est :
vol
$all
    AIC 
220.662 

$moinsUn
     variable   AIC               
[1,] "surface"  "219.100682981217"
[2,] "pers"     "224.198835542612"
[3,] "pavillon" "225.012502253576"
[4,] "age"      "218.819173225281"
[5,] "vol"      "218.744689237871"
[6,] "sbains"   "218.869400969036"

> 
> # Seconde etape : on retire "vol"
> AIC.retirer(y = "kw", regresseurs = c("surface", "pers", "pavillon", "age", "sbains"))
La variable faisant le plus baisser l'AIC est :
sbains
$all
     AIC 
218.7447 

$moinsUn
     variable   AIC               
[1,] "surface"  "241.929604226733"
[2,] "pers"     "225.287648259288"
[3,] "pavillon" "223.468684246265"
[4,] "age"      "226.861973666485"
[5,] "sbains"   "216.930365210544"

> 
> # Troisieme etape : on retire "sbains"
> AIC.retirer(y = "kw", regresseurs = c("surface", "pers", "pavillon", "age"))
Le retrait d'une variable fait augmenter l'AIC, on conserve le modele
kw~surface+pers+pavillon+age
$all
     AIC 
216.9304 

$moinsUn
     variable   AIC               
[1,] "surface"  "240.671160290049"
[2,] "pers"     "223.423119145426"
[3,] "pavillon" "221.526143791873"
[4,] "age"      "224.977206891741"

> 
> # Le modele retenu est :
> kw~surface+pers+pavillon+age
kw ~ surface + pers + pavillon + age

La sélection en backward en se fiant au critère AIC désigne le modèle dans lequel la consommation d'électricité est fonction linéaire de la surface, du nombre de personnes, de la dummy pavillon et de l'âge du logement.

> summary(lm(kw~surface+pers+pavillon+age, data = df))

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

Residuals:
    Min      1Q  Median      3Q     Max 
-686.63 -115.60   20.73  265.88  442.32 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -522.588    399.740  -1.307  0.21376    
surface       39.741      6.182   6.428 2.24e-05 ***
pers        -420.354    150.147  -2.800  0.01504 *  
pavillon     584.532    243.692   2.399  0.03217 *  
age           14.403      4.621   3.117  0.00817 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 369 on 13 degrees of freedom
Multiple R-squared:  0.9383,	Adjusted R-squared:  0.9193 
F-statistic: 49.44 on 4 and 13 DF,  p-value: 9.715e-08

Ce qui est intéressant ici, c'est que le modèle proposé souffre vraisemblablement de multicolinéarité, à en croire le billet précédent. Les données semblent étranges aussi : il est difficile d'admettre que l'augmentation d'une personne dans le logement provoque (toutes choses égales par ailleurs) une diminution de la consommation d'électricité de 420 KWH en janvier... Prudence donc, lors de l'utilisation de ce genre d'algo !

2. Sélection de variables en forward à l'aide du F-partiel de Fisher

Cette fois, il s'agit d'ajouter pas-à-pas les variables. On régresse dans un premier temps la variable à expliquer sur chacune des explicatives. On récupère les statistiques de Student pour la variable ajoutée. Si le coefficient est statistiquement non nul (au seuil de 5% par exemple), la variable est candidate pour entrer dans le modèle. Celle qui entre effectivement est celle pour laquelle le F-partiel (le carré de la statistique de Student) est le plus élevé, parmi les variables candidates. S'il n'y a pas de variable candidate, on s'arrête et on conserve le modèle auquel l'ajout d'une variable ne produit pas un meilleur modèle au yeux du F-partiel.

Étape 1

On effectue les régression de la consommation d'électricité sur chacune des variables explicatives.

> reg.surface <- lm(kw ~ surface, data= df)
> reg.pers <- lm(kw ~ pers, data= df)
> reg.pavillon <- lm(kw ~ pavillon, data= df)
> reg.age <- lm(kw ~ age, data= df)
> reg.vol <- lm(kw ~ vol, data= df)
> reg.sbains <- lm(kw ~ sbains, data= df)
> 
> summary(reg.surface)$coef
               Estimate Std. Error   t value     Pr(>|t|)
(Intercept) -1299.83986 531.326506 -2.446405 2.636170e-02
surface        36.98337   4.407128  8.391718 2.966392e-07
> summary(reg.pers)$coef
             Estimate Std. Error    t value    Pr(>|t|)
(Intercept)  79.54321   952.7412 0.08348879 0.934498528
pers        703.49630   220.8456 3.18546646 0.005750654
> summary(reg.pavillon)$coef
            Estimate Std. Error  t value     Pr(>|t|)
(Intercept) 2052.625   347.6887 5.903628 2.223649e-05
pavillon    1724.675   466.4734 3.697263 1.953624e-03
> summary(reg.age)$coef
              Estimate Std. Error  t value     Pr(>|t|)
(Intercept) 2249.30472  379.69992 5.923901 2.139556e-05
age           33.43052   12.14699 2.752166 1.417082e-02
> summary(reg.vol)$coef
              Estimate Std. Error   t value     Pr(>|t|)
(Intercept) -522.26569 368.200252 -1.418428 1.752500e-01
vol           11.19035   1.106711 10.111363 2.354069e-08
> summary(reg.sbains)$coef
            Estimate Std. Error  t value    Pr(>|t|)
(Intercept) 1252.912   650.1045 1.927246 0.071896438
sbains      1265.664   430.7522 2.938264 0.009643102

Au seuil de \(5\%\), le coefficient associé à la variable ajoutée est non nul, pour chacune des régressions. Aussi, toutes les variables sont candidates à l'ajout dans le modèle. On va effectivement ajouter celle pour qui le F-partiel est le plus élevé.

> summary(reg.surface)$coef["surface","t value"]^2
[1] 70.42094
> summary(reg.pers)$coef["pers","t value"]^2
[1] 10.1472
> summary(reg.pavillon)$coef["pavillon","t value"]^2
[1] 13.66976
> summary(reg.age)$coef["age","t value"]^2
[1] 7.574419
> summary(reg.vol)$coef["vol","t value"]^2
[1] 102.2397
> summary(reg.sbains)$coef["sbains","t value"]^2
[1] 8.633397

Parmi les variables candidates, "vol" est celle pour qui le F-partiel est le plus élevé. On l'intègre donc à notre modèle.

Étape 2

Le modèle de base est à présent :
\[\textrm{conso}_i = \alpha_0 + \alpha_1 \textrm{vol}_i + u_i, \quad i=1,\ldots,n\]
avec \(u_i\) un terme d'erreur d'espérance nulle et de variance \(\sigma^2_u\).

On ajoute à ce modèle une variable explicative, et on regarde si le coefficient associé est significativement différent de zéro. S'il l'est, la variable est candidate à l'ajout définitif dans le modèle.

> reg.vol.surface <- lm(kw ~ vol + surface, data= df)
> reg.vol.pers <- lm(kw ~ vol + pers, data= df)
> reg.vol.pavillon <- lm(kw ~ vol + pavillon, data= df)
> reg.vol.age <- lm(kw ~ vol + age, data= df)
> reg.vol.sbains <- lm(kw ~ vol + sbains, data= df)
> 
> summary(reg.vol.surface)$coef
                Estimate Std. Error     t value   Pr(>|t|)
(Intercept) -541.5646525 569.316457 -0.95125417 0.35655258
vol           10.9839880   4.672733  2.35065596 0.03283628
surface        0.7245757  15.908175  0.04554738 0.96427186
> summary(reg.vol.pers)$coef
              Estimate Std. Error   t value     Pr(>|t|)
(Intercept)   23.68536 391.131539  0.060556 9.525123e-01
vol           14.31799   1.601213  8.941965 2.130525e-07
pers        -368.01988 150.257027 -2.449269 2.708023e-02
> summary(reg.vol.pavillon)$coef
             Estimate Std. Error    t value     Pr(>|t|)
(Intercept) -373.3628  379.25546 -0.9844625 3.405021e-01
vol           10.0507    1.40293  7.1640771 3.260525e-06
pavillon     379.6400  296.21307  1.2816452 2.194248e-01
> summary(reg.vol.age)$coef
               Estimate Std. Error    t value     Pr(>|t|)
(Intercept) -479.346785 389.644565 -1.2302155 2.375616e-01
vol           10.834571   1.384441  7.8259542 1.127812e-06
age            3.047261   6.784967  0.4491195 6.597663e-01
> summary(reg.vol.sbains)$coef
               Estimate Std. Error     t value     Pr(>|t|)
(Intercept) -523.112047 386.336870 -1.35403087 1.957759e-01
vol           11.178670   1.480647  7.54985735 1.744042e-06
sbains         3.265424 263.042125  0.01241407 9.902589e-01

Seul le coefficient de la variable "pers" est statistiquement non-nul au seuil de \(5\%\). Il n'y a donc qu'une variable candidate, qui est donc d'office ajoutée au modèle.

Étape 3

Le modèle de base devient :
\[\textrm{conso}_i = \alpha_0 + \alpha_1 \textrm{vol}_i + \alpha_2 \textrm{pers}_i + u_i, \quad i=1,\ldots,n\]
avec \(u_i\) un terme d'erreur d'espérance nulle et de variance \(\sigma^2_u\).

On estime les modèles augmentés d'une variable pour sélectionner d'éventuelles variables à l'ajout au modèle.

> reg.vol.pers.surface <- lm(kw ~ vol + pers + surface, data= df)
> reg.vol.pers.pavillon <- lm(kw ~ vol + pers + pavillon, data= df)
> reg.vol.pers.age <- lm(kw ~ vol + pers + age, data= df)
> reg.vol.pers.sbains <- lm(kw ~ vol + pers + sbains, data= df)
> 
> summary(reg.vol.pers.surface)$coef
              Estimate Std. Error    t value   Pr(>|t|)
(Intercept) -336.40946 477.476743 -0.7045567 0.49264098
vol            9.80491   3.895359  2.5170749 0.02463973
pers        -457.98491 163.570697 -2.7999203 0.01418255
surface       18.53048  14.637546  1.2659550 0.22619041
> summary(reg.vol.pers.pavillon)$coef
              Estimate Std. Error   t value     Pr(>|t|)
(Intercept)  399.08759 364.085621  1.096137 2.915170e-01
vol           13.30992   1.420693  9.368607 2.082373e-07
pers        -461.96845 133.247809 -3.466987 3.774968e-03
pavillon     601.78056 233.852386  2.573335 2.209367e-02
> summary(reg.vol.pers.age)$coef
              Estimate Std. Error    t value     Pr(>|t|)
(Intercept)  110.87874 389.812946  0.2844409 7.802364e-01
vol           16.60489   2.398378  6.9233821 7.056138e-06
pers        -512.74382 186.758874 -2.7454857 1.578426e-02
age           -9.05270   7.175330 -1.2616423 2.276937e-01
> summary(reg.vol.pers.sbains)$coef
              Estimate Std. Error     t value     Pr(>|t|)
(Intercept)   14.14652  403.78803  0.03503452 9.725468e-01
vol           14.08628    1.76728  7.97059847 1.429391e-06
pers        -376.36848  156.47199 -2.40534094 3.055469e-02
sbains        84.58675  231.51360  0.36536405 7.203006e-01

Seul le coefficient associé à la variable "pavillon" est statistiquement non-nul au seuil de \(5\%\). La variable candidate est la seule candidate et est de fait ajoutée au modèle.

Étape 4

Le modèle de base devient :
\[\textrm{conso}_i = \alpha_0 + \alpha_1 \textrm{vol}_i + \alpha_2 \textrm{pers}_i + \alpha_3 \textrm{pavillon}_i + u_i, \quad i=1,\ldots,n\]
avec \(u_i\) un terme d'erreur d'espérance nulle et de variance \(\sigma^2_u\).

> reg.vol.pers.pavillon.surface <- lm(kw ~ vol + pers + pavillon + surface, data= df)
> reg.vol.pers.pavillon.age <- lm(kw ~ vol + pers + pavillon + age, data= df)
> reg.vol.pers.pavillon.sbains <- lm(kw ~ vol + pers + pavillon + sbains, data= df)
> 
> summary(reg.vol.pers.pavillon.surface)$coef
              Estimate Std. Error    t value    Pr(>|t|)
(Intercept)  147.97440 470.587325  0.3144462 0.758171741
vol           10.63635   3.440921  3.0911337 0.008592159
pers        -509.14123 145.425061 -3.5010556 0.003906131
pavillon     551.88934 243.240152  2.2689072 0.040951187
surface       11.32069  13.243258  0.8548264 0.408125227
> summary(reg.vol.pers.pavillon.age)$coef
               Estimate Std. Error    t value     Pr(>|t|)
(Intercept)  412.693132 370.945583  1.1125436 2.860549e-01
vol           14.608089   2.306340  6.3338844 2.601073e-05
pers        -530.441631 165.431482 -3.2064128 6.882100e-03
pavillon     549.590585 248.677316  2.2100552 4.564627e-02
age           -4.792756   6.634644 -0.7223834 4.828523e-01
> summary(reg.vol.pers.pavillon.sbains)$coef
              Estimate Std. Error    t value     Pr(>|t|)
(Intercept)  419.29324  382.70212  1.0956126 2.931311e-01
vol           13.44188    1.53555  8.7537917 8.233328e-07
pers        -459.25569  138.11878 -3.3250778 5.477021e-03
pavillon     623.10476  252.37075  2.4690054 2.818741e-02
sbains       -61.21502  206.83903 -0.2959549 7.719392e-01

Le coefficient de chaque variable ajoutée n'est pas statistiquement différent de zéro au seuil de \(5\%\). Il n'y a plus de variable candidate.

Modèle retenu

En procédant avec la méthode ascendante (ou forward), et en se fiant au critère du F-partiel, le "meilleur" modèle est le suivant :

\[\textrm{conso}_i = \alpha_0 + \alpha_1 \textrm{vol}_i + \alpha_2 \textrm{pers}_i + \alpha_3 \textrm{pavillon}_i + u_i, \quad i=1,\ldots,n\]
avec \(u_i\) un terme d'erreur d'espérance nulle et de variance \(\sigma^2_u\).

Une fonction pour aller un peu plus vite, juste pour le fun

> fPartiel.Ajouter <- function(base, regresseurs){
+   estimationUneVarEnPlus <- function(base, unRegresseur){
+     laFormule <- formula(paste(base,"+",unRegresseur))
+     reg.tmp <- lm(laFormule, data = df)
+     
+     # le t de Student observe
+     tObs <- summary(reg.tmp)$coefficients[unRegresseur,"t value"]
+     
+     # La p-value du test de Student de significativite du coef
+     pVal <- summary(reg.tmp)$coefficients[unRegresseur,"Pr(>|t|)"]
+     
+     return(c(tObs = tObs, pVal = pVal))
+   }
+   testStud <- do.call("rbind", lapply(regresseurs, function(x) estimationUneVarEnPlus(base, x)))
+   testStud <- data.frame(testStud)
+   testStud$fPartiel <- testStud$tObs^2
+   
+   # On va regarder quels sont les coefficients statistiquement differents de 0 au seuil de 5%
+   IND <- which(testStud$pVal < 0.05)
+   
+   if(length(IND)>0){# On a des candidats
+     # Le fPartiel le plus eleve
+     fPartielMax <- max(testStud[IND,"fPartiel"])
+     # Ce max est pour la variable :
+     varRetenue <- regresseurs[which(testStud$fPartiel == fPartielMax)]
+     cat("La variable pour laquelle l'effet est significatif a 5% et qui a le f-partiel le plus eleve est :", varRetenue, sep = "\n")
+   }else{# Aucune des variable ajoutees n'a de coefficient correspondant statistiquement non nul a 5%
+     cat("L'ajout d'une variable ne produit pas un meilleur modele.\n On conserve le modele suivant :",
+         base, sep = "\n")
+   }
+ }
> fPartiel.Ajouter(base  <- "kw ~ 1", regresseurs = c("surface", "pers", "pavillon", "age", "vol", "sbains"))
La variable pour laquelle l'effet est significatif a 5% et qui a le f-partiel le plus eleve est :
vol
> fPartiel.Ajouter(base  <- "kw ~ vol", regresseurs = c("surface", "pers", "pavillon", "age", "sbains"))
La variable pour laquelle l'effet est significatif a 5% et qui a le f-partiel le plus eleve est :
pers
> fPartiel.Ajouter(base  <- "kw ~ vol + pers", regresseurs = c("surface", "pavillon", "age", "sbains"))
La variable pour laquelle l'effet est significatif a 5% et qui a le f-partiel le plus eleve est :
pavillon
> fPartiel.Ajouter(base  <- "kw ~ vol + pers + pavillon", regresseurs = c("surface", "age", "sbains"))
L'ajout d'une variable ne produit pas un meilleur modele.
 On conserve le modele suivant :
kw ~ vol + pers + pavillon

L'estimation du modèle retenu

> summary(reg.fPart <- lm(kw ~ vol + pers + pavillon, data = df))

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

Residuals:
    Min      1Q  Median      3Q     Max 
-649.33 -182.79   71.39  180.75  543.86 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  399.088    364.086   1.096  0.29152    
vol           13.310      1.421   9.369 2.08e-07 ***
pers        -461.968    133.248  -3.467  0.00377 ** 
pavillon     601.781    233.852   2.573  0.02209 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 366.7 on 14 degrees of freedom
Multiple R-squared:  0.9344,	Adjusted R-squared:  0.9203 
F-statistic: 66.45 on 3 and 14 DF,  p-value: 1.598e-08

La lecture de la sortie de régression nous donne les informations suivantes :

  • Le test de Fisher de significativité globale (sous \(H_0\), tous les coefficients sauf celui de la constante égaux à zéro, et au moins un différent de zéro sous \(H_1\)) indique que la régression ne semble pas mauvaise (au seuil de \(5\%\)), comparativement à celle dans laquelle le consommation est régressée en fonction d'une constante seule ;
  • la valeur du \(R^2\) est de \(0.9344\), ce qui indique que \(93.44\%\) des variations de la consommation d'électricité sont expliquées par le modèle ;
  • le coefficient associé à la variable "vol" est statistiquement différent de zéro, au seuil de \(5\%\) (La probabilité, si on émet l'hypothèse que le coefficient associé à la variable "vol" est nul, d'observer une valeur, lors d'un tirage à partir d'une variable aléatoire de distribution Student à \(18-3-1=14\) degrés de liberté au moins aussi grande à la valeur absolue du \(t_obs = 9.369\) est de \(2.08e-07\) ; de fait on rejette l'hypothèse de nullité du coefficient associé à la variable "vol" au seuil de \(5\%\).) Une interprétation de ce coefficient peut être la suivante : lorsque le volume du logement augmente (diminue) d'un mètre cube, toutes choses égales par ailleurs, la consommation augmente (diminue) de \(13.31\) KWH ;
  • le coefficient associé à la variable "pers" est différent de zéro, au seuil de \(5\%\) (on a en effet \(p-value < 0.05\)). L'interprétation laisse cependant penser qu'il y a un problème avec le modèle, ou avec les données : une personne supplémentaire (en moins) dans le logement entraîne, toutes choses égales par ailleurs, une diminution (augmentation) de la consommation d'électricité en janvier de \(461.968\) KWH. Ce résultat ne semble pas du tout intuitif, il faudrait creuser un peu pour comprendre ce qui se passe ici... ;
  • le coefficient associé à la variable "pavillon" (qui est une variable indicatrice, ou encore appelée dummy variable) est significativement différent de zéro, au seuil de \(5\%\). L'interprétation diffère légèrement du cas où la variable est quantitative : le fait d'être dans un pavillon ("\(\textrm{pavillon}=1\)") fait augmenter la consommation d'électricité, toutes choses égales par ailleurs, de \(601.781\) KWH, comparativement au fait d'être dans un appartement ("\(\textrm{pavillon}=0\)").

3. Sélection de variables en stagewise

Il s'agit d'une méthode de type forward, dans laquelle on choisit d'introduire la variable qui va expliquer le mieux la part de la réponse \(y\) non expliquée par les variables déjà introduites (donc qui explique le mieux les résidus !). L'algorithme que l'on va suivre est constitué de 4 étapes :

  1. on initie \(e=y\) ;
  2. on choisit la variable \(x_j\) la plus corrélée avec \(e\). Si le coefficient associé à cette variable est significativement non-nul au sens du \(t_2\) (défini ci-dessous), à un seuil donné (par ex. \(5\%\)), on l'introduit. Sinon, on s'arrête ;
  3. on calcule la part de \(y\) non-expliquée par les variables déjà sélectionnées / intro-
    duites : \(e = y− (\hat{\alpha_0} + \hat{\alpha}_j x_j + \ldots)\) ;
  4. on retourne en étape 2.

Le test à effectuer s'appuie sur la statistique suivante :
\[t_2 = \frac{\rho_{je}}{\sqrt{\frac{1-\rho^2_{je}}{n-(l+1)}}} \sim \mathcal{S}t(n-(l+1),\]
avec \(\rho_{je}\) le coefficient de corrélation de la variable \(x_j\) avec \(e\), \(n\) le nombre d'observations, \(l\) le nombre de variables explicatives dans le modèle.

Étape 1 (initialisation)

On définit \(e = y\).

Étape 2

Calculons les corrélations entre \(e\) et chaque variable explicative :

> 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

La variable la plus corrélée avec \(e\) est "vol". La statistique \(t2\) observée est :

> n <- nrow(df) ; l <- 1
> (rho.je <- cor(kw,vol))
[1] 0.9298826
> (t2 <- rho.je / (sqrt((1-rho.je^2) / (n - (l+1)))))
[1] 10.11136

que l'on compare au quantile d'une Student à 16 degrés de liberté, pour une probabilité de \(1-0.05/2 = 0.975\) :

> # Valeur tabulee
> qt(p = 1-0.05/2, df = n-(l+1))
[1] 2.119905

On rejette donc \(H_0\) au seuil de \(5\%\) (le coefficient de la variable "vol" est significativement différent de zéro au sens du t^2).

Étape 3

On redéfinit \(e = \textrm{conso} - (\hat{\alpha}_0 + \hat{\beta}_{\textrm{vol}} \textrm{vol})\).

reg <- lm(kw ~ vol, data = df)
e <- df$kw - (reg$coef["(Intercept)"] + reg$coef["vol"] * df$vol)
# Ou de maniere equivalente :
e <- residuals(reg) # les residus de la regression

Retour en étape 2

On choisit la variable \(x_j\) la plus corrélée avec notre nouvelle variable \(e\) :

> cor(data.frame(e, surface, pers, pavillon, age, sbains))
                    e     surface       pers  pavillon        age      sbains
e         1.000000000 0.002876307 -0.3224664 0.2429989 0.09446898 0.002474351
surface   0.002876307 1.000000000  0.8373894 0.6764150 0.36355769 0.659700224
pers     -0.322466436 0.837389371  1.0000000 0.6333333 0.15242287 0.574960988
pavillon  0.242998903 0.676414970  0.6333333 1.0000000 0.11605802 0.588984427
age       0.094468976 0.363557687  0.1524229 0.1160580 1.00000000 0.200186507
sbains    0.002474351 0.659700224  0.5749610 0.5889844 0.20018651 1.000000000

La variable la plus corrélée avec \(e\) est "pers". Regardons la valeur du \(t2\) :

> n <- length(e) ; l <- 2
> (rho.je <- cor(e,pers))
[1] -0.3224664
> (t2 <- rho.je / (sqrt((1-rho.je^2) / (n - (l+1)))))
[1] -1.319388

que l'on compare au quantile \(t_{1-0.05/2}^{15}\) :

> qt(p = 1-0.05/2, df = n-(l+1))
[1] 2.13145

On a donc \(\mid t2 \mid < t_{1-0.05/2}^{15}\). De fait, la variable "pers" n'est pas significative au sens du \(t2\). On s'arrête donc ici, et le modèle retenu est alors : \[\textrm{conso}_i = \alpha_0 + \alpha_1 \textrm{vol}_i + u_i, \quad i=1,\ldots,n\] avec \(u_i\) un terme d'erreur d'espérance nulle et de variance \(\sigma^2_u\). Comme on peut le voir, selon le critère retenu, on peut aboutir à différents modèle. Il appartient donc à la discrétion de la personne qui modélise de faire certains choix arbitraires.

2 thoughts on “[L3 Eco-Gestion] Régression linéaire avec R : sélection de modèle

  1. Bonjour,

    Le fait que le coefficient associé à la variable « pers » soit négatif (égal à -461.968 kW/pers dans le modèle lm(kw ~ vol + pers + pavillon, data = df), cf. deuxième méthode (F-partiel de Fisher)), peut s’expliquer assez simplement d’un point de vue physique (thermodynamique) :
    Dans un foyer moyen, la plus grande partie du coût énergétique (en kWh) est due au chauffage.
    Si toutes choses égales par ailleurs, même volume, même surface, même âge (, et même isolation), on ajoute une personne dans la maison, cela fait un corps de mammifère à 37°C qui est nettement plus chaud que le sol et les murs.
    On a un « radiateur » biologique de plus, partant, on a moins besoin de chauffer.

    Mat2boc

  2. Bonjour,

    Est-ce que la méthode de Sélection de variables en forward à l’aide du F-partiel de Fisher permet de donner une solution (au moins partielle) sur la Multi-colinéarité, dire autrement, est-ce que en choisissant les variables de cette façon permet de rassembler les variables les moins colinéaires?

    Merci par avance.

    Z. Dhima

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.