Tutoriel : la régression linéaire multiple avec R

Pour faire suite à mon article d’introduction à la régression linéaire multiple, nous allons voir ici, dans ce tutoriel,  comment réaliser ce type d’analyse statistique, avec le logiciel R.

J’emploie généralement une approche en 10 étapes :

  1. Etude des corrélations linéaires entre les variables explicatives deux à deux (pairwise)
  2. Evaluation de la linéarité entre chacune des variables explicatives (ou indépendantes) de type numérique et la variable réponse (ou dépendante)
  3. Ajustement du modèle complet
  4. Evaluation des multi-collinéarités par les VIFs
  5. Evaluation des hypothèses de normalité et d’homoscédasticité des résidus du modèle complet
  6. Interprétation des résultats du modèle complet
  7. Sélection du modèle parcimonieux
  8. Recherche des outliers
  9. Interprétation du modèle parcimonieux
  10. Synthèse des résultats

Pour illustrer cette approche, nous allons employer le jeu de données mtcars (présent dans R).

Il s’agit de données extraites du magazine américain Motor Trend de 1974, elles comprennent la consommation de carburant et 10 aspects de la conception et des performances automobiles pour 32 automobiles (modèles 1973-74) :

  • mpg : Miles/(US) gallon
  • cyl : Number of cylinders
  • disp : Displacement (cu.in.)
  • hp : Gross horsepower
  • drat : Rear axle ratio
  • wt : Weight (1000 lbs)
  • qsec : 1/4 mile time (le temps pour parcourir 1/4 de mile)
  • vs : Engine (0 = V-shaped, 1 = straight)
  • am : Transmission (0 = automatic, 1 = manual)
  • gear : Number of forward gears
  • carb : Number of carburetors.

Nous allons employer une régression linéaire multiple  afin d’évaluer quelles variables ont un effet sur la distance parcourue avec un gallon d’essence (c’est la variable réponse), indépendamment des autres variables explicatives.

Etude des corrélations linéaires entre les variables explicatives deux à deux

Lorsque deux variables explicatives très fortement corrélées sont incluses ensemble dans le modèle de régression linéaire multiple, cela peut le rendre instable : des relations non significatives peuvent le devenir subitement lorsqu’une variable est retirée.

Lorsque des fortes corrélations linéaires  sont mises en évidence (j’emploie le seuil de 0.85, en termes de coefficient de corrélation de Pearson), une des deux variables (celle qui a le moins de sens “métier”, ou celle qui est le moins liée à la réponse), ne devra pas être incluse dans le modèle de régression linéaire multiple.

Pour étudier les corrélations linéaires entre les variables explicatives deux à deux, j’emploie généralement la fonction ggpairs du package GGally, parce sa syntaxe est simplissime :

library(GGally)
ggpairs(mtcars) 
corrélation deux à deux

Ici, nous pouvons voir qu’il existe plusieurs fortes corrélations linéaires (r > 0.85):

  • entre disp et cyl : r= 0.902
  • entre wt et cyl : r= 0.88

Il existe encore une relation avec un coefficient de 0.83 entre cyl et hp.

Je fais alors le choix de retirer disp car cette variable a, pour moi qui ne connais pas grand-chose aux voitures, moins de sens que le poids (wt).

library(tidyverse)

mtcars2 <- mtcars %>% 
    select(-disp) 

Je ne fais pas de choix ente cyl et hp car le coefficient est <0.85, et si besoin, je retirerai d’autres variables lors de l’étape 4, avec la vérification des VIFs

Evaluation de la linéarité entre la réponse et les variables explicatives numériques

Pour cela, j’utilise généralement la fonction scatterplotMatrix du package car. Mais avant cela, je vais limiter les variables explicatives aux seules variables numériques :

# retire les variables binaires
mtcars3 <- mtcars %>% 
    select(-vs, -am)

library(car)
scatterplotMatrix(mtcars3 
Etude des linéarités

On s’intéresse aux graphiques de la première ligne, qui représentent la relation entre la variable réponse mpg en y, et chacune des variables explicatives, en x.

Remarque : pour que les graphiques à analyser soient sur la première ligne il faut que la variable réponse (ici mpg) soit sur la première colonne du data frame.

Il s’agit de regarder (en utilisant le zoom) si les relations sont globalement linéaires. Je dirais que c’est plutôt le cas ici, avec quelques défauts (qui me semblent négligeables) au niveau des variables disp et gear.

Si les relations sont curvilinéaires (courbées) vous pouvez essayer de les linéariser en utilisant une transformation log de la variable prédictive. Une autre solution peut consister à catégoriser cette variable prédictive. Ou encore d’ajouter le terme au carré pour modéliser cette forme curvilinéaire (mais cela est ensuite plus difficile à intrepréter).

Ajustement du modèle complet.

Avant d’ajuster le modèle de régression linéaire multiple complet (contenant toutes les variables explicatives), je passe les variables vs et am en facteur avec des niveaux V-shaped et straight pour la variable vs, et automatic et manual pour la variable am. C’est seulement pour que les sorties que nous obtiendrons soient plus faciles à interpréter.

library(tidyverse)
# remplace les valeurs 0 et 1
mtcars2 <- mtcars2 %>% 
    mutate(vs =ifelse(vs==0, "V-shaped","straight"),
           am = ifelse(am==0, "automatic","manual"))

# passage en facteur
mtcars2 <- mtcars2 %>% 
    mutate_if(is.character, as.factor)
str(mtcars2)
## 'data.frame':    32 obs. of  10 variables:
##  $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
##  $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
##  $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
##  $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
##  $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
##  $ qsec: num  16.5 17 18.6 19.4 17 ...
##  $ vs  : Factor w/ 2 levels "straight","V-shaped": 2 2 1 1 2 1 2 1 1 1 ...
##  $ am  : Factor w/ 2 levels "automatic","manual": 2 2 2 1 1 1 1 1 1 1 ...
##  $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
##  $ carb: num  4 4 1 1 2 1 4 2 2 4 ... 

Une fois ces modifications réalisées, nous pouvons ajuster le modèle complet :

#ajustement du modèle complet
mod.rlm1 <- lm(mpg~., data=mtcars2) 

Le terme ~. permet d’inclure toutes les variables explicatives présentent dans les données mtcars2.

Pour vérifier que tout s’est bien déroulé, nous pouvons afficher la table de régression , en utilisant la fonction summary()

summary(mod.rlm1)
## 
## Call:
## lm(formula = mpg ~ ., data = mtcars2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7863 -1.4055 -0.2635  1.2029  4.4753 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 12.71088   18.83577   0.675   0.5068  
## cyl          0.09627    0.99715   0.097   0.9240  
## hp          -0.01295    0.01834  -0.706   0.4876  
## drat         0.92864    1.60794   0.578   0.5694  
## wt          -2.62694    1.19800  -2.193   0.0392 *
## qsec         0.66523    0.69335   0.959   0.3478  
## vsV-shaped  -0.16035    2.07277  -0.077   0.9390  
## ammanual     2.47882    2.03513   1.218   0.2361  
## gear         0.74300    1.47360   0.504   0.6191  
## carb        -0.61686    0.60566  -1.018   0.3195  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.623 on 22 degrees of freedom
## Multiple R-squared:  0.8655, Adjusted R-squared:  0.8105 
## F-statistic: 15.73 on 9 and 22 DF,  p-value: 1.183e-07 

Il est nécessaire de vérifier :

  • que tous les paramètres ont été estimés : il ne doit pas avoir de NA dans cette table de régression
  • que les erreurs standard des pentes partielles ne sont pas très grands par rapport aux pentes partielles. Si c’est le cas, c’est qu’e vous avez ’il existe un phénomène de multi-collinéarité. Celui-ci pourra être objectivé, à la prochaine étape, celle de la vérification des VIFs.

Ici tout est OK.

Evaluation des multicollinéarités par les VIF

La page d’aide de la fonction check_collinearity() du package performance fournit une description détaillée de la multicolinéarité. Voici une traduction, un peu simplifiée.

La multicolinéarité ne doit pas être confondue avec une forte corrélation brute entre les prédicteurs. Ce qui compte, c’est l’association entre une ou plusieurs variables prédictives, conditionnellement aux autres variables du modèle

Alors, pour être honnête; jusqu’à ce que je lise cette phrase, pour moi forte corrélation brute en les prédicteurs et multicollinéarité était la même chose…

La multicolinéarité signifie qu’une fois que l’on connaît l’effet d’une variable explicative sur la réponse, la valeur à connaître, au niveau de l’autre variable prédictive, est faible.Autrement dit, en cas de multicolinéarité, l’une des variables prédictive n’aide pas beaucoup en termes de meilleure compréhension du modèle ou de prédiction du résultat.En cas de multicolinéarité, le modèle semble suggérer que les variables explicatives en question ne sont pas associées à la réponse (pentes partielles faibles, erreurs types élevées), alors que ces variables sont en réalité fortement associés à la réponse, c’est-à-dire qu’elles pourraient effectivement avoir un effet important (McElreath 2020, chapitre 6.1).Les VIFs (Variance Inflation Factor) permettent de mesurer l’ampleur de la multicolinéarité des variables explicatives incluses dans le modèle. La racine carrée du VIF indique de combien de fois l’erreur standard d’une pente partielle est augmentée à cause de la multicollinéarité.Un VIF inférieur à 5 indique une faible corrélation de cette variable explicative avec d’autres variables explicatives. Une valeur comprise entre 5 et 10 indique une corrélation modérée, tandis que des valeurs VIF supérieures à 10 sont le signe d’une corrélation élevée et non tolérable (James et al. 2013).

Lorsque qu’une variable à un VIF > 10, il est nécessaire de la retirer du modèle,  puis de recalculer les VIFs, et de retirer une seconde variable si nécessaire, etc… jusqu’à n’obtenir que des VIFs <5.

Le seuil de 10 ne fait pas forcément consensus, j’utilise plutôt 5.

La fonction check_collinearity(), du package performance permet d’obtenir ces VIFs très facilement :

library(performance)
check_collinearity(mod.rlm1)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##  Term  VIF Increased SE Tolerance
##  drat 3.33         1.82      0.30
##    vs 4.92         2.22      0.20
##    am 4.65         2.16      0.22
##  carb 4.31         2.08      0.23
## 
## Moderate Correlation
## 
##  Term  VIF Increased SE Tolerance
##    hp 7.12         2.67      0.14
##    wt 6.19         2.49      0.16
##  qsec 6.91         2.63      0.14
##  gear 5.32         2.31      0.19
## 
## High Correlation
## 
##  Term   VIF Increased SE Tolerance
##   cyl 14.28         3.78      0.07 

On retire donc la variable cyl du modèle (nous avions vu à l’étape 1 qu’elle était fortement liée aux variables wt et hp , ce n’est donc pas une surprise.).

Pour retirer la variable du modèle nous pouvons employer la fonction update(), comme ceci :

mod.rlm2 <- update(mod.rlm1, .~.-cyl)
check_collinearity(mod.rlm2)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##  Term  VIF Increased SE Tolerance
##  drat 3.11         1.76      0.32
##    vs 4.27         2.07      0.23
##    am 4.29         2.07      0.23
##  gear 4.69         2.17      0.21
##  carb 4.29         2.07      0.23
## 
## Moderate Correlation
## 
##  Term  VIF Increased SE Tolerance
##    hp 6.02         2.45      0.17
##    wt 6.05         2.46      0.17
##  qsec 5.92         2.43      0.17 

Puisque la variable wt a du sens pour moi, je vais d’abord retirer la variable hp :

mod.rlm3 <- update(mod.rlm2, .~.-hp)
check_collinearity(mod.rlm3)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##  Term  VIF Increased SE Tolerance
##  drat 3.04         1.74      0.33
##  qsec 4.14         2.03      0.24
##    vs 4.19         2.05      0.24
##    am 4.26         2.06      0.23
##  gear 4.69         2.17      0.21
##  carb 3.83         1.96      0.26
## 
## Moderate Correlation
## 
##  Term  VIF Increased SE Tolerance
##    wt 5.10         2.26      0.20 

Le VIF de la variable wt est tout juste > 5, je décide de la garder dans le modèle.

Je vérifie que les variables ont bien été retirées en affichant la table de régression :

summary(mod.rlm3)
## 
## Call:
## lm(formula = mpg ~ drat + wt + qsec + vs + am + gear + carb, 
##     data = mtcars2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9187 -1.1587 -0.1858  1.3021  4.3141 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   8.3579    11.7109   0.714  0.48230   
## drat          1.0565     1.4897   0.709  0.48504   
## wt           -2.9502     1.0543  -2.798  0.00997 **
## qsec          0.8955     0.5198   1.723  0.09782 . 
## vsV-shaped    0.1033     1.8548   0.056  0.95605   
## ammanual      2.5377     1.8883   1.344  0.19155   
## gear          0.6730     1.3400   0.502  0.62006   
## carb         -0.7573     0.5530  -1.370  0.18350   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.542 on 24 degrees of freedom
## Multiple R-squared:  0.8623, Adjusted R-squared:  0.8221 
## F-statistic: 21.46 on 7 and 24 DF,  p-value: 6.989e-09 

Evaluation des hypothèses de normalité et d’homoscédasticité des résidus

Là encore le package performance propose des fonctions très utiles :

check_model() : qui réalise un diagnostic de régression à l’aide de 6 graphiques, et qui permet d’évaluer les hypothèses de linéarité, d’homoscédasticité et de normalité des résidus, ainsi que les multi-collinéarité et les valeurs influentes.

check_model(mod.rlm3) 
diagnostic de régression

Comme nous l’avions vu lors d’une précédente étape, la linéarité souffre d’un léger défaut, mais il me semble acceptable.

La normalité et l’homscédasticité des résidus ne sont pas parfaites, mais me semblent également raisonnable.

Vous pouvez réaliser un test de Shapiro-Wilk pour tester la normalité des résidus en employant la fonction check_normality() :

check_normality(mod.rlm3) 

La normalité n’est pas rejetée (p=042).

Vous pouvez aussi réaliser un test de Breusch-pagan pour tester l’hypothèse d’homoscédasticité des résidus en employant la fonction check_heteroscedasticity():

check_heteroscedasticity(mod.rlm3)
 

L’homoscédasticité n’est pas rejetée (p=0.073).

Interprétation des résultats du modèle complet

La table de régression

La table de régression peut être obtenue avec la fonction summary() :

summary(mod.rlm3)
## 
## Call:
## lm(formula = mpg ~ drat + wt + qsec + vs + am + gear + carb, 
##     data = mtcars2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9187 -1.1587 -0.1858  1.3021  4.3141 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   8.3579    11.7109   0.714  0.48230   
## drat          1.0565     1.4897   0.709  0.48504   
## wt           -2.9502     1.0543  -2.798  0.00997 **
## qsec          0.8955     0.5198   1.723  0.09782 . 
## vsV-shaped    0.1033     1.8548   0.056  0.95605   
## ammanual      2.5377     1.8883   1.344  0.19155   
## gear          0.6730     1.3400   0.502  0.62006   
## carb         -0.7573     0.5530  -1.370  0.18350   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.542 on 24 degrees of freedom
## Multiple R-squared:  0.8623, Adjusted R-squared:  0.8221 
## F-statistic: 21.46 on 7 and 24 DF,  p-value: 6.989e-09 

Une seule variable,  le poids (wt),  apparait significativement liée à la distance parcourue avec un gallon d’essence (mpg) : plus la voiture est lourde, moins la distance parcourue est importante (la pente est négative, beta=-2.9). Ainsi, lorsque toutes les autres variables sont fixées, la distance parcourue diminue de pratiquement 3 miles lorsque le poids de la voiture augmente d’une unité (ici 1000 livres).

A noter que la p-value associée à la variable qsec (je n’ai pas trop compris de quoi il s’agit…) est < 0.1. Cette variable se retrouvera probablement dans le modèle de régression parcimonieux.

La table d'analyse de variance

Lorsque les variables sont de type catégoriel avec k modalités, la table de régression fournit k-1 p-values. Pour obtenir la p-value de l’effet global de la variable, nous devons obtenir une table d’analyse de la variance, comme ceci :

car::Anova(mod.rlm3)
## Anova Table (Type II tests)
## 
## Response: mpg
##            Sum Sq Df F value   Pr(>F)   
## drat        3.251  1  0.5030 0.485035   
## wt         50.601  1  7.8294 0.009973 **
## qsec       19.178  1  2.9674 0.097821 . 
## vs          0.020  1  0.0031 0.956051   
## am         11.672  1  1.8060 0.191553   
## gear        1.630  1  0.2523 0.620059   
## carb       12.123  1  1.8758 0.183496   
## Residuals 155.111 24                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Remarque : attention à ne pas employer la fonction anova() car elle fournit des résultats qui dépendent de l’ordre d’apparition des variables dans le modèle.

Sélection et interprétation du modèle parcimonieux

Pour cela, nous retirons de manière itérative, les variables les moins significativement liées à la réponse. Dans le modèle global, c’est la variable vs car sa pvalue est la plus élevée (p=0.95)

mod.rlm4 <- update(mod.rlm3, .~.-vs)
Anova(mod.rlm4)
## Anova Table (Type II tests)
## 
## Response: mpg
##            Sum Sq Df F value   Pr(>F)   
## drat        3.233  1  0.5210 0.477115   
## wt         55.254  1  8.9044 0.006274 **
## qsec       34.151  1  5.5036 0.027211 * 
## am         12.012  1  1.9357 0.176392   
## gear        1.621  1  0.2612 0.613804   
## carb       12.108  1  1.9513 0.174727   
## Residuals 155.131 25                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Puis, nous retirons la variable gear:

mod.rlm5 <- update(mod.rlm4, .~.-gear)
Anova(mod.rlm5)
## Anova Table (Type II tests)
## 
## Response: mpg
##            Sum Sq Df F value   Pr(>F)   
## drat        4.498  1  0.7461 0.395620   
## wt         71.232  1 11.8151 0.001988 **
## qsec       38.839  1  6.4421 0.017477 * 
## am         20.097  1  3.3335 0.079391 . 
## carb       11.134  1  1.8468 0.185826   
## Residuals 156.752 26                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Puis, nous retirons la variable drat:

mod.rlm6 <- update(mod.rlm5, .~.-drat)
Anova(mod.rlm6)
## Anova Table (Type II tests)
## 
## Response: mpg
##            Sum Sq Df F value    Pr(>F)    
## wt        104.754  1 17.5401 0.0002686 ***
## qsec       54.371  1  9.1040 0.0055070 ** 
## am         33.281  1  5.5726 0.0257208 *  
## carb        8.036  1  1.3456 0.2562120    
## Residuals 161.250 27                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Puis, nous retirons la variable carb:

mod.rlm7 <- update(mod.rlm6, .~.-carb)
Anova(mod.rlm7)
## Anova Table (Type II tests)
## 
## Response: mpg
##            Sum Sq Df F value    Pr(>F)    
## wt        183.347  1 30.3258 6.953e-06 ***
## qsec      109.034  1 18.0343 0.0002162 ***
## am         26.178  1  4.3298 0.0467155 *  
## Residuals 169.286 28                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Toutes les variables restantes sont significativement liées à la distance parcourue avec un gallon d’essence. Nous pouvons afficher les estimations des pentes partielles : 

summary(mod.rlm7)
## 
## Call:
## lm(formula = mpg ~ wt + qsec + am, data = mtcars2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4811 -1.5555 -0.7257  1.4110  4.6610 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.6178     6.9596   1.382 0.177915    
## wt           -3.9165     0.7112  -5.507 6.95e-06 ***
## qsec          1.2259     0.2887   4.247 0.000216 ***
## ammanual      2.9358     1.4109   2.081 0.046716 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.459 on 28 degrees of freedom
## Multiple R-squared:  0.8497, Adjusted R-squared:  0.8336 
## F-statistic: 52.75 on 3 and 28 DF,  p-value: 1.21e-11 

Nous pouvons alors vérifier les conditions d’application sur ce modèle parcimonieux :

check_model(mod.rlm7) 
évaluation des conditions d'application
check_heteroscedasticity(mod.rlm7)
## OK: Error variance appears to be homoscedastic (p = 0.212). 
check_normality(mod.rlm7)
## OK: residuals appear as normally distributed (p = 0.087). 

Tout est satisfaisant.

Recherche des outliers

Il peut être important de vérifier qu’il n’y a pas d’outliers dans ce modèle de régression, afin d’évaluer si les résultats n’ont pas été influencés par certains points.

Pour cela, nous pouvons employer la fonction check_outliers() du packages performance. Plusieurs approches sont employées ; pour plus d’informations, consultez l’aide de cette fonction.

check_outliers(mod.rlm7)
## OK: No outliers detected 

Si des outliers sont mis en évidence, il est nécessaire de refaire l’analyse en les retirant afin d’évaluer la robustesse des résultats (est ce que les résultats restent identiques ?).

Pour cela vous pouvez utiliser la fonction update() , en indiquant les indices des outliers dans l’argument subset, comme cela par exemple : update(mod.rlm, subset=c(3,12,15)).

 

Interprétation du modèle parcimonieux

summary(mod.rlm7)
## 
## Call:
## lm(formula = mpg ~ wt + qsec + am, data = mtcars2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4811 -1.5555 -0.7257  1.4110  4.6610 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.6178     6.9596   1.382 0.177915    
## wt           -3.9165     0.7112  -5.507 6.95e-06 ***
## qsec          1.2259     0.2887   4.247 0.000216 ***
## ammanual      2.9358     1.4109   2.081 0.046716 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.459 on 28 degrees of freedom
## Multiple R-squared:  0.8497, Adjusted R-squared:  0.8336 
## F-statistic: 52.75 on 3 and 28 DF,  p-value: 1.21e-11 

Trois variables sont significativement liées à la distance parcourue :

  • le poids : lorsque celui-ci augmente d’une unité (1000 livres), la distance parcourue avec un gallon d’essence diminue d’environ 4 miles.
  • qsec : la distance parcourue augmente en moyenne de 1.23 miles quand le temps pour parcourir 1/4 de mile) augmente d’une unité
  • am  (la transmission) : la distance parcourue augmente d’environ 3 miles lorsque la transmission est manuelle par rapport à une transmission automatique

Présentation des résultats

Selon les besoins, les résultats peuvent être présentés de différentes façons. Dans certaines situations, on ne présentera que la table de régression du modèle parcimonieux. Dans d’autres (de plus en plus fréquentes), les pentes (partielles), leurs intervalles de confiance et les pvalue correspondantes sont présentés pour ces différents modèles de régressions :

  • les modèles de régression univariés (chaque variable explicative séparément)
  • le modèle de régression linéaire multiple complet
  • le modèle de régression linéaire parcimonieux.

Pour obtenir ces estimations, vous pouvez employer les fonctions lmuni(), et lmmulti() du package finalfit, comme ceci :

library(finalfit)

explanatory = c("drat", "qsec", "vs", "am", "gear","carb","wt")
dependent = 'mpg'

# estimation des pentes, IC et pvalues des modèles de régression univarié (une seule variable explicative)
mtcars2 %>%
  lmuni(dependent, explanatory)  %>% 
    fit2df() 
##   explanatory                      Coefficient
## 1        drat    7.68 (4.60 to 10.76, p<0.001)
## 2        qsec     1.41 (0.27 to 2.55, p=0.017)
## 3  vsV-shaped -7.94 (-11.27 to -4.61, p<0.001)
## 4    ammanual    7.24 (3.64 to 10.85, p<0.001)
## 5        gear     3.92 (1.25 to 6.59, p=0.005)
## 6        carb  -2.06 (-3.22 to -0.89, p=0.001)
## 7          wt  -5.34 (-6.49 to -4.20, p<0.001)

# estimation des pentes partielle, IC et pvalues dans le modèle de régression multiple complet (toutes les variables explicatives)
mtcars2 %>%
  lmmulti(dependent, explanatory)  %>% 
    fit2df()
##   explanatory                     Coefficient
## 1        drat   1.06 (-2.02 to 4.13, p=0.485)
## 2        qsec   0.90 (-0.18 to 1.97, p=0.098)
## 3  vsV-shaped   0.10 (-3.72 to 3.93, p=0.956)
## 4    ammanual   2.54 (-1.36 to 6.43, p=0.192)
## 5        gear   0.67 (-2.09 to 3.44, p=0.620)
## 6        carb  -0.76 (-1.90 to 0.38, p=0.183)
## 7          wt -2.95 (-5.13 to -0.77, p=0.010)

# estimation des pentes partielles, IC et pvalues dans le modèle de régression parcimonieux  (uniquement les variables explicatives sélectionnées)
explanatory_final = c("qsec", "am", "wt")
dependent = 'mpg'
mtcars2 %>%
  lmmulti(dependent, explanatory_final)  %>% 
    fit2df()
##   explanatory                     Coefficient
## 1        qsec    1.23 (0.63 to 1.82, p<0.001)
## 2    ammanual    2.94 (0.05 to 5.83, p=0.047)
## 3          wt -3.92 (-5.37 to -2.46, p<0.001) 

Conclusion

A ma connaissance, il n’ya pas de consensus absolu dans la manière de mener une régression linéaire multiple. Dans ce tutoriel je vous ai montré les 10 étapes que je mène. Mais si vous avez une autre routine, n’hésitez pas à la partager en commentaire.

 

Et comme d’habitude, si cet article vous a plus, ou vous a été utile, et si vous le souhaitez, vous pouvez soutenir ce blog en faisant un don sur sa page Tipeee 

12 réponses

Laisser un commentaire

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

Bonjour !

vous venez souvent ?

Identifiez-vous pour avoir accès à toutes les fontionnalités !

Aide mémoire off'R ;)

Enregistrez vous pour recevoir gratuitement mes fiches “aide mémoire” (ou cheat sheets) qui vous permettront de réaliser facilement les principales analyses biostatistiques avec le logiciel R et pour être informés des mises à jour du site.