Régressions linéaires avec R

Étape 0 - Partons d'un exemple ! Imaginons un échantillon d'une population, dont on a la taille, le poids et l'IMC (Indice de Masse Corporelle) !

Essayons de voir si on peut retrouver la relation entre taille, poids et IMC sachant qu'on connait pour une fois la réponse IMC = poids/(taille/100)²

Simulons ce jeu de données :

nb = 180 ; poids <- round(rnorm(nb,65,11),0)  ; imc <- rnorm(nb,22,5) ; # IMC 

taille <- round(sqrt(poids/imc)*100,0) ; taille_bis <- round(rnorm(nb,170,11),0)

taille_biss <- round(rnorm(nb,170,11),0) ; taille <- round((taille+taille_bis+taille_biss)/3,0)

imc <- round(poids/(taille/100)^2,1)

compil <- data.frame(poids,taille,imc) # Compilation des valeurs

1.2 - Décrire les relations statistiques entre les variables X1, X2... et Y avec les fonctions cor(), cor.test() et le R²

• Établir la significativité de l'influence de chaque paramètre avec cor.test()

# Je teste une relation linéaire entre poids et imc

cor.test(compil$imc,compil$poids)

cor.test(compil$imc,compil$taille)

La p-value renvoyé à chaque fois montre ici un lien significatif entre poids et imc et légèrement significatif entre taille et IMC. 

Remarque : la méthode par défaut ici est la méthode de Pearson. Il en existe d'autres.

• Analyser la forme de la relation avec le R² (coefficient de détermination)

Si R² = 1, le nuage est très allongé, X et Y tracent une droite.

Si R² = 0, le nuage se disperse dans tous les sens et une régression linéaire ne permettra pas de faire des prédictions.

# Je calcule mon R²

cor(compil$imc,compil$poids)^2

cor(compil$imc,compil$taille)^2

A ce stade, la relation poids - IMC est bien établie. On voit en revanche que la relation taille - IMC a un R² proche de 0.

Faut-il renoncer à la taille dans notre régression ? On sait bien que non ! Cela nous invite à rester méfiant et à ne pas exclure une variable si le R² seul ou la p-value nous pousse à renoncer.

Extraire les p-values de chaque coefficient :

summary(reg1)[[4]][,4]

Extraire la p-value finale de la régression

str(summary(reg1))

# Extraction de la p-value

pf(summary(reg1)$fstatistic[1],summary(reg1)$fstatistic[2],summary(reg1)$fstatistic[3],lower.tail=FALSE)

En paraphrasant statistique-et-logiciel-r.com, on sait qu'on peut valider une relation linéaire entre X et Y si les résidus :

1. Sont indépendants

2. Sont distribués selon une loi Normale de moyenne 0

3. Sont distribués de façon homogènes, c’est à dire, avec une variance constante.

C'est quoi un résidu ? C'est la différence entre valeur réelle et sa projection sur la droite de régression : entre pratique et théorie.

Tous les graphiques nécessaires à ce qui va suivre s'affichent en une seule fois avec la commande plot(reg1).

La fonction valreg() pour valider un modèle de régression

On peut s'éviter toutes les étapes qui suivent en utilisant la fonction valreg() qui compile cet ensemble :

• Installer valreg() (package {KefiR}

install.packages("devtools") ; require(devtools) # Risque d'erreur si RTools non installé.

devtools::install_github("Antoine-Masse/KefiR")

• Exécuter le package {KefiR}

library("KefiR")

• Exemple d'utilsation de valreg()

valreg(reg, plot=T, analyse=T)

• plot : permet d'activer le paramètre boot (boostrap) de validation du modèle de régression avec affichage d'un bilan graphique.

• analyse : permet d'afficher les commentaires et les raisons de l'invalidation du modèle de régression.

3.1- Indépendance des résidus

Les résidus ne doivent pas être liés les uns aux autres. En gros : on ne doit pas pouvoir anticiper en résidu par la connaissance d'un autre résidu.

On veut éviter les phénomènes d'auto-corrélation : les résidus augmenteraient ensemble dans une zone donnée et seraient liés les uns aux résidus précédents ou suivants (cas typiques lorsque X est un enregistrement du temps).

Vérification par le test de Durbin-Watson (indépendance : p-value > 0,05)

Cela se vérifie avec le test de Durbin-Watson : il faut une p-value supérieure à 0,05 pour avoir indépendance.

durbinWatsonTest(reg1) # nécessite la librairie car

dwtest(reg1) # nécessite la librairie lmtest

Attention : 0,05 n'est qu'une limite, il ne faut pas renoncer si on a 0,03 !

Autre remarque : ce test n'évalue que l'autocorrélation entre résidus et le résidus qui suit directement (n+1) et non ceux séparés d'un intervalle (résidus+2) qui pourraient être autocorrélés dans un phénomène cyclique.

Vérification par le test de Shapiro-Wilk (normalité : p-value > 0,05)

Un test permet de trancher sur la normalité de la distribution des résidus : c'est le test de normalité de Shapiro-Wilk

shapiro.test(residuals(reg1))

On peut considérer qu'il y a normalité tant que la p-value reste supérieure à 0,05 (mais on le rappelle que c'est une frontière floue : on peut tolérer 0,03).

Vérification par le test de Breush-Pagan (homogénéité : p-value > 0,05)

Le test de  Breush-Pagan permet de vérifier objectivement l'homogénéité des résidus (et donc l'hétéroscédasticité vs homoscédasticité).

install.packages("lmtest") ; library("lmtest") ; 

bptest(reg1) # librairie lmtest

ncvTest(reg1) # librairie car

L'homogénéité est rejeté si la p-value est inférieure à 0,05. Dans cet exemple, la relation poids/imc serait acceptée.

Vérification par le test de Goldfeld-Quandt (homogénéité : p-value > 0,05)

D'autres sources suggèrent de faire le test de Goldfeld-Quandt.

library("lmtest")

gqtest(reg1)

D'autres proposent un 3ème test, celui de White supposé plus capable sur les modèles complexes à détecter l'hétéroscédasticité mais moins puissant que Breusch-Pagan sur les modèles simples.

Voici quelques indices qui permettent de prouver qu'un modèle de régression est bon.

Certains sont pertinents pour les régressions multiples et comparer plusieurs modèles. C'est le cas du RSS, de l'AIC, du BIC, du MSE.

En revanche, biais, variance et R² ajusté peuvent déjà  s'utiliser sur une simple régression à une seule variable.

• Le RSS (Residual Sum of Squares), la somme des résidus au carré, illustre la pertinence d'un modèle de régression.

Plus le RSS est petit, mieux c'est !

Le risque du RSS est qu'il peut être bon si notre régression est absurde est modélisant le bruit. Autre inconvénient, le RSS varie avec le nombre d'individus.

• Le R² ajusté, R squared adjusted. Ce R² doit tendre vers 1. Il se calcule à partir du RSS. Mais attention ! On peut avoir un R² proche de 1 si on a compilé un grand nombre de variable et fait le modèle à partir d'un seul échantillon sans jamais le tester.

Plus le R² ajusté est proche de 1, mieux c'est !  

• L'AIC (Critère d'Information d'Akaike) est plus pertinent que le RSS car elle va simuler un RSS si on avait eu des données inconnues. Elle va aussi pénaliser l'augmentation du nombre de variables.

Plus l'AIC est petit, mieux c'est !

• Le BIC (Critère d'Information Bayésien) pénalise encore plus que l'AIC le nombre de variables. On doit donc trouver un compromis entre le nombre de variables et le BIC le plus bas.

Plus le BIC est petit, mieux c'est !

• Le MSE (Mean Sqare Error). Équivalent au RSS, le MSE est calculé sur un échantillon test que l'on va appliquer au modèle établi à partir d'un échantillon d'apprentissage. La MSE tient compte du nombre d'individus. Critère pertinent !

Plus le MSE est petit, mieux c'est !

• Le biais, distance entre le modèle et la réalité. C'est plus souvent d'ailleurs la différence entre le modèle et l'échantillon.

Plus le biais est petit, mieux c'est !

• La variance, potentialité des coefficients du modèle à bouger d'un échantillon à l'autre.

Moins ils bougent, mieux c'est !

Nous reparlerons plus tard du RSS, de l'AIC, du BIC et du MSE...

On ne parlera dans l'immédiat que du R² ajusté et de la variance.

4.1- Valider la pertinence des coefficients par les p-values et R² ajusté

1. Si le coefficient intercept (la valeur b dans une équation y = ax + b), et les autres coefficients directeurs sont associés à des valeurs Pr petites (et plusieurs étoiles), alors ils ont un impact significatif et sont pertinents.

2. Le R² ajusté vaut ici de 0,5869 ce qui n'est pas énorme. On pourrait améliorer le modèle de régression pour approcher 1.

3. La p-value très basse (2.2e-16) montre aussi que le modèle établit une relation pertinente entre poids et IMC. Ce lien ne peut être nié mais il peut probablement être ajusté.

4.2- La variance du modèle de régression

Lorsque l'on a plusieurs variables, la fonction lm() reste pertinente. Il faut juste en maîtriser la syntaxe.

• lm(y~x) : régression classique y = f(x)  - affine - y = ax+b

• lm(y~x+0) : régression linéaire y = f(x) sans intersection - y = ax

• lm(y~x1+x2) : régression telle que y = B + A1*x1 + A2*x2

• lm(y~x1:x2) : régression telle que y = B + A1*x1*x2 

• lm(y~x1*x2) : régression telle que y = B + A1*x1 + A2*x2 + A3*x1*x2

A cette syntaxe, on peut ajouter d'autres fonctions pour permettre des régression polynomiales avec poly() ou d'autres transformations avec I().

Exemples :

• lm(y~poly(x,3)) : régression polynomiale telle que y = B + A1*x + A2*x² + A3*x^3

• lm(y~I(1/x) : permet d'imposer des transformations à x comme ici y = B + A1*1/x 

Créons donc des modèles de régressions pour essayer de retrouver la relation entre IMC et poids / taille

reg2 <- lm(imc~poids+taille,data=compil)

reg3 <- lm(imc~poids*taille,data=compil)

reg4 <- lm(imc~poids*poly(taille,2),data=compil)

reg5 <- lm(imc~poids:taille,data=compil)

reg6 <- lm(imc~poids*I(taille^2),data=compil)

reg7 <- lm(imc~poids*I(1/(taille/100)^2)+0,data=compil)

reg8 <- lm(imc~poids:I(1/(taille/100)^2)+0,data=compil)# C'est le vrai modèle où IMC = poids/taille_en_cm² - pas de valeur B

La fonction coef() permet d'afficher l'ensemble des coefficients de la régression.

coef(reg4)

6.1- Comparer les R² ajusté. 

On y accède avec summary() sinon on peut regarder les R² non ajustés (très proches, moins fiables).

summary(reg1)  # 0.5892249

summary(reg2)  # 0.9868673

summary(reg3)  # 0.9958284

summary(reg4)  # 0.9998804

summary(reg5)  # 0.354403

summary(reg6)  # 0.992617

summary(reg7)  # 1

summary(reg8)  # 1

On voit ici que le modèle reg1 qui ne dépend que du poids n'est pas pertinent à côté des autres, tout comme reg5.

Il faudrait maintenant comparer les modèles reg2, reg3, reg4, reg6, reg7 et reg8.

Il faut surtout ne pas prendre celui qui a le meilleur R² ajusté ! Car on pourrait avoir un modèle trop complexe derrière. Mieux vaut un modèle un peu moins bon mais simple !

6.2- Comparer les modèles de régression avec la fonction anova().

6.3- Comparer les modèles de régression par l'intermédiaire des RSS.

Calculer les RSS des modèles : ce calcul peut se faire de la façon suivante - sans passer par la fonction anova().

sum(reg7$residuals^2) # sommes des résidus mis au carré pour le modèle reg7

6.4- Comparer les modèles de régression par l'intermédiaire des AIC.

Pour comparer l'AIC, on peut utiliser les fonctions extractAIC() et AIC(). Bien que ces deux fonctions ne renvoient pas la même valeurs d'AIC. Pas d'inquiétude, c'est juste qu'elle n'additionnent pas à la base de leurs calculs la même constante.

L'essentiel est de toujours utiliser la même fonction.

extractAIC() renvoie deux valeurs : la première, le nombre de variables et l'AIC qui, rappelons-le, doit être aussi petit que possible.

extractAIC(reg1) # 2.0000 296.1014 

extractAIC(reg2) # 3.0000 -321.6283

extractAIC(reg3) # 4.0000 -526.0523

extractAIC(reg4) # 6.000 -1161.387

extractAIC(reg5) # 2.0000 377.4847

extractAIC(reg6) # 4.0000 -423.2939

extractAIC(reg7) # 3.000 -1266.462

extractAIC(reg8) # 1.000 -1269.127

On voit que le modèle reg4 et reg7 jugés aussi bon que la réalité (reg8) ont pourtant 4 à 6 variables contre une seule ! (en réalité 2 car la formule poids/taille^2 est interprétée comme une seule variable dans reg8)

6.5- Comparer les modèles de régression par l'intermédiaire des BIC.

On peut comparer les BIC de chaque modèle avec la fonction BIC.

On rappelle que le BIC pénalise encore plus le nombre de variables que l'AIC.

BIC(reg1) # 878.7482

BIC(reg2) # 186.8721

BIC(reg3) # 128.4066

BIC(reg4) # -516.1518

BIC(reg5) # 936.2825

BIC(reg6) # 231.072

BIC(reg7) # -742.7517

BIC(reg8) # -753.0865

Il apparait très clairement que les seuls modèles à même de rivaliser avec la vérité (reg8) sont reg4 et reg7 avec une nette préférence pour reg7 qui est un modèle plus simple que reg4.

6.6- Comparer la colinéarité des variables 

La fonction vif() de la librairie car permet de comparer la colinéarité des variables. Les valeurs renvoyées par vif doivent être basses sinon les variables ont tendance à évoluer dans la même direction.

library(car)

vig(reg2)  # pas de colinéarité

vif(reg4) # colinéarité

vif(reg7) # colinéarité

Le mieux est d'en calculer directement la racine carré avec sqrt() : cela permet de savoir dans quelle mesure l'écart-type d'une variable est influencé par sa colinéarité avec une autre. On peut, par exemple, se fixer un seuil de 10 et éliminer les variables qui dépassent ce seuil.

Ainsi  : Mes analyses précédentes (RSS, AIC, BIC) montrent que reg7 serait un bon modèle.

Les valeurs renvoyées par vif()  suggèrent de se débarrasser de "poids:I(1/(taille/100)^2)" et on retrouverait ainsi la vérité, reg8 !

6.7- Ajouter ou supprimer des variables avec les fonctions add1(), drop1() et update() 

La fonction update permet de modifier un modèle existant en y ajoutant ou en y supprimant des variables.

Un exemple ou je vais supprimer la variable espèce d'une régression : 

data(iris)

lm(Sepal.Length~.,data=iris) -> reg

regressor<-update(reg,~.-Species)

summary(regressor)

La fonction drop1() ou son pendant (add1()) permettent de voir ce qui se passe si on ajoute ou supprime une variable. Il ne change pas réellement le modèle mais nous donnes des infos (AIC, BIC, p-values...).

Son grand avantage est lorsqu'on a fait une régression en fonction d'une variable catégoriel, lm() nous fera une régression en fonction de chaque catégorie : drop1() va regrouper tout cela pour clarifier le problème...

Exemple : j'ai fait une variable sur le jeu de données iris en fonction notamment de la variable Species. Il y a 3 espèces (versicolor, virginica et setosa) : lm()) va faire une régression pour chaque espèce selon la logique "pas versicolor" = 0 et "versicolor" = 1.

De ce fait, Species apparait plusieurs fois dans le modèle. drop1() évite ce problème.

data(iris)

lm(Sepal.Length~.,data=iris) -> reg

summary(reg)

drop1(reg) 

Partons volontairement d'une régression sale ! J'ai fait une régression croisée polynomiale où Y dépend de X1, X1^2, X2, X2^2 et des produits des ces éléments mais aussi des inverses...

reg9 <- lm(imc~poly(poids,2)*poly(taille,2)*poly(I(1/taille),2)*poly(I(1/poids),2))

La fonction step() est capable de nettoyer ces modèles. step() va avancer pas à pas pour supprimer des coefficients et réoptimiser les autres sans que la régression ne perte de sa qualité.

On peut paramétrer strep() de 3 façons :

• forward : on part d'un premier coefficient et on en ajoute peu en voyant jusqu'où s'arrêter

• backward : on part de l'ensemble des coefficients et on en retire peu à peu

• both : un mixe entre forward et backward

summary(reg9)

reg9_prim <-  step(reg9,direction="both")

summary(reg9_prim)

anova(reg9,reg9_prim) # l'anova ne montre pas d'amélioration de reg9_prim par rapport à reg9

extractAIC(reg9)

extractAIC(reg9_prim)

On voit ici que l'AIC a baissé un peu et que j'ai réussi de passer de 25 variables à 9 variables ! J'en ai éliminé 16 !

La fonction update() permet d'ajouter ou de retirer manuellement des variables.

summary(reg9_prim)

reg9_update <- update(reg9_prim,.~.-1) # je retire l'intercept

summary(reg9_update)

reg9_update <- update(reg9_update,.~.-poly(poids, 2)) # je retire les polynomes du poids

reg9_update <- update(reg9_update,.~.-poly(I(1/taille), 2)) # je retire les polynomes de l'inverse de la taille

reg9_update <- update(reg9_update,.~.+I(poids/(taille/100)^2)) # j'introduis la formule réelle du modèle

extractAIC(reg9_prim) ; extractAIC(reg9_update) # j'ai bien baissé mon AIC

reg9_update <- step(reg9_update) # je peux encore utiliser step pour éliminer ce qui peut l'être

summary(reg9_update) # Je suis ainsi arrivé à retrouver mon modèle

Pour valider un modèle, il existe plusieurs approches :

1. Validation avec échantillon d'apprentissage et échantillon test - on s'assure alors que le MSE soit bas. Et que le modèle marche pour l'échantillon test.

2. Validation croisée : un bootstrap permet de refaire plusieurs fois la méthode ci-dessus.

Une fois le modèle validé :

Un bootstrap va permettre de remixer l'échantillon qui a servi pour l'obtention du modèle (en retirant certains points et en dupliquant d'autres).

On va ainsi se faire une idée précise de la variance et du biais.

Comment faire : il suffit de demander à être formé

Voici aussi une petite illustration du potentielle du machine learning : le bootstrap

10.1- Régression en 2D

On peut très facilement afficher une régression linéaire si on a (sans compter l'intercept) 1 seul ou 2 coefficients.

Si j'ai une seule variable et une relation du type y = ax + b, je peux utiliser la fonction abline().

D'autres fonctions ci-dessus sont très utiles telles matplot() ou encore les affichages de régressions obtenues avec Zelig (cf. ci-dessus).

10.3- Et si j'ai plus de 2 variables ?

Impossible de les voir !!! Mais ce problème reste le même sans régression. Il faudra alors penser à l'ACP ou sinon, modéliser pour deux variables et fixer les valeurs des autres variables.

Exemple : si j'avais une variable supplémentaire X3

Grille<-expand.grid(poids=my_seq_poids,taille=my_seq_taille)

Grille <- data.frame(Grille,X3=rep(12,40^2)) # Permet de fixer la variable X3 à 12.

Il ne reste qu'à suivre le code ci-dessus avec la fonction predict()...

Mais attention : si les autres variables ont un rôle à jouer dans les prédictions : cela ne donnera rien.

Étape 11 - Gérer les variables colinéaires

Certains variables sont colinéaires. Parfois aussi, on a plusieurs Y (Y1, Y2, Y3) à mettre en relation avec plusieurs X (X1, X2, X3)...

Comment faire : il suffit de demander à être formé