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 :
- Etude des corrélations linéaires entre les variables explicatives deux à deux (pairwise)
- 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)
- Ajustement du modèle complet
- Evaluation des multi-collinéarités par les VIFs
- Evaluation des hypothèses de normalité et d’homoscédasticité des résidus du modèle complet
- Interprétation des résultats du modèle complet
- Sélection du modèle parcimonieux
- Recherche des outliers
- Interprétation du modèle parcimonieux
- 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) 
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 
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) 
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) 
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
24 réponses
Article de très grande qualité ! Meci beaucoup pour ce très joli et très informatif article !
J’aimerais accéder continuellement à ce article très enrichissant.
Merci
super travail. Félicitations
Merci beaucoup Claire, j’ai beaucoup appris avec cet article !!!
Merci beaucoup pour ce travail formidable ca m’a permis de découvrir de nouvelles fonctions très important.
Bravo
Waooh très riche… merci beaucoup pour cet article…
Super comme pédagogie. Merci
Vraiment merci pour cet article.
Très riche en substance !
Article très intéressant. Les explications sont très claires instructives
Merci c’est vraiment bien fait
Super formation. Très intéressant. Merci beaucoup.
Merci !
Merci beaucoup pour ce travail ! C’est très clair et pédagogique, merci de prendre du temps pour le faire ! 🙂
Bonjour Claire,
quelle la règle sur le nombre de variables max à intégrer dans un modèle quand nous avons peu d’observation (n=50)?
# modèle de régression linéaire mixte
# modèle logistique
# modèle de survie
Merci pour votre aide.
Julie
merci pour le partage mais la fonction check_model ne marche pas! je reçois seulement dans ma consolece message : “”Error: The RStudio ‘Plots’ window is too small to show this set of plots.
Please make the window larger.” (j’ai bien sur elargie la fenetre plot mais toujours sans succes
merci pour la réponse
Bonjour Moundir,
Le message vous indique que la fenêtre graphique n’est pas assez grande, il suffit de l’agrandir (surtout en hauteur il me semble).
Bonne continuation.
Bonjour Claire.
Merci pour cet article très complet et très instructif. J’ai cependant un problème avec la fonction check_model(). On me demande à chaque fois si je dois installer le package {see} (il est déjà installé) et j’ai une erreur du type : “Error in `colnames<-`(`*tmp*`, value = `*vtmp*`) :
tentative de modification de 'colnames' sur un objet ayant moins de deux dimensions". Comment y rémédier ? Merci d'avance.
Bonjour Justin,
je n’ai pas de piste, navrée.
Bonne continuation
Bonsoir,
Je sais que je réponds à un post de plus de 6 mois mais sait-on jamais si cela peut servir.
J’ai eu la même à peu près après l’installation de “Performance” :
check_model(mod.rlm3)
# Package `see` required for model diagnostic plots.
# Would you like to install it? [y/n] y
# Error in names(x) <- value: 'names' attribute [4] must be the same length as the vector [3]
En insistant une nouvelle fois :
check_model(mod.rlm3)
# Package `patchwork` required for this function to work.
# Would you like to install it? [y/n] y
Et là la fonction check_model() fonctionne. Peut être un bug dans le package. Peut être qu'aller chercher Patchwork manuellement pourrait aider… Bonne chance.
Merci pour cet article.
Bonsoir madame Claire,
Vous avez développé un cours assez intéressant.
Néanmoins pour le modèle à retenir, il y a un pakage qui permet de procéder par élimination des variables de façon dégressive ou progressive ou mixte et en fin on a un modèle parcimonieux au lieu qu’on procède cette élimination manuellement.
Joël ADANDONON
Biostatisticien Econometre.
Bonjour Joël,
Merci pour votre message.
Pouvez-vous, s’il vous plait, nous indiquer quel package vous employez pour réaliser la sélection de variables.
Dans le domaine médical, nous avons tendance à ne pas réaliser de sélection automatique, mais cela reste bien évidemment intéressant.
Bonjour,
Ma question concerne les tests statistiques (non abordés par votre sujet)
1/ pour etudier l’hypothèse de non colinéarité, le logiciel R fournit il le test des régression auxiliaires?
2/ et au niveau de l’hypothèse de linéarité, il y a le test de Rainbow. Mais malgré mes recherches, je n’ai pas trouvé une presentation détaillé de ce test. Est ce que vous avez plus d’informations à ce propos?
Et merci pour tout ce que vous faites
Bonjour,
Je suis navrée de ne pas pouvoir vous aider davantage, mais je ne connais aucun des 2 points mentionnés.
Je pense que les régressions auxiliaires sont liées aux VIF, mais je ne sais pas davantage.