Tutoriel · Intermédiaire
Modèles de régression sous R : tests, diagnostics et interprétation
Une démarche complète de régression linéaire appliquée, de l’examen des données à la présentation de résultats statistiquement fondés.
Consommation de carburant dans le jeu de données mtcars
Le modèle étudié est mpg ~ wt + hp. Chaque section explique ce que la procédure teste, comment l’exécuter et quelles conclusions son résultat permet de tirer.
- Observations
- 32
- Variable expliquée
- MPG
- Variables explicatives
- Poids + puissance
Repères graphiques
Voir les idées clés avant les tests
Ces figures schématiques montrent ce que l’on cherche dans les graphiques produits par R.Fondements
Préparer R et examiner les données
Ce tutoriel modélise la consommation de carburant, exprimée en miles par gallon (mpg), à partir du poids (wt, en milliers de livres) et de la puissance (hp) des véhicules du jeu de données mtcars intégré à R. Commencez par vérifier les types des variables, les valeurs manquantes, les statistiques descriptives et les valeurs impossibles.
data(mtcars)
head(mtcars[c("mpg", "wt", "hp")])
str(mtcars[c("mpg", "wt", "hp")])
summary(mtcars[c("mpg", "wt", "hp")])
colSums(is.na(mtcars[c("mpg", "wt", "hp")]))> str(mtcars[c("mpg", "wt", "hp")])
'data.frame': 32 obs. of 3 variables:
$ mpg: num 21 21 22.8 21.4 18.7 ...
$ wt : num 2.62 2.88 2.32 3.21 3.44 ...
$ hp : num 110 110 93 110 175 ...
> colSums(is.na(mtcars[c("mpg", "wt", "hp")]))
mpg wt hp
0 0 0Les trois variables du modèle sont numériques et ne comportent aucune observation manquante. Un résumé satisfaisant ne prouve pas que les données sont exemptes d’erreurs : comparez toujours les plages de valeurs et les unités au contexte de l’étude.
Fondements
Explorer les relations avant l’estimation
Les nuages de points et les corrélations révèlent le sens et la forme approximative des relations, les observations inhabituelles et les liens étroits entre les variables explicatives avant d’imposer une structure linéaire.
pairs(mtcars[c("mpg", "wt", "hp")],
pch = 19, col = adjustcolor("#2563eb", 0.55))
cor(mtcars[c("mpg", "wt", "hp")])
plot(mpg ~ wt, data = mtcars, pch = 19, col = "#2563eb")
abline(lm(mpg ~ wt, data = mtcars), col = "#dc2626", lwd = 2)> round(cor(mtcars[c("mpg", "wt", "hp")]), 3)
mpg wt hp
mpg 1.000 -0.868 -0.776
wt -0.868 1.000 0.659
hp -0.776 0.659 1.000
Sortie graphique : matrice de nuages de points et nuage de mpg en fonction
du poids, avec une droite ajustée décroissante.Dans ces données, mpg est fortement corrélé négativement au poids (−0,868) et à la puissance (−0,776). Le poids et la puissance sont positivement corrélés (0,659) : la multicolinéarité devra donc être examinée.
Fondements
Estimer la régression linéaire multiple
Le modèle estime la moyenne conditionnelle de mpg comme une fonction linéaire du poids et de la puissance : mpgᵢ = β₀ + β₁wtᵢ + β₂hpᵢ + εᵢ.
model <- lm(mpg ~ wt + hp, data = mtcars)
summary(model)
coef(model)
fitted(model)
residuals(model)Call:
lm(formula = mpg ~ wt + hp, data = mtcars)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 37.22727 1.59879 23.285 < 2e-16
wt -3.87783 0.63273 -6.129 1.12e-06
hp -0.03177 0.00903 -3.519 0.00145
Residual standard error: 2.593 on 29 degrees of freedom
Multiple R-squared: 0.8268; Adjusted R-squared: 0.8148
F-statistic: 69.21 on 2 and 29 DF, p-value: 9.109e-12L’équation estimée est mpg = 37,227 − 3,878(wt) − 0,0318(hp). À puissance constante, 1 000 livres supplémentaires sont associées à une baisse d’environ 3,88 miles par gallon. À poids constant, un cheval supplémentaire est associé à une baisse d’environ 0,032 mpg.
Inférence statistique
Tester les coefficients individuels avec les tests t
Le test t de chaque coefficient évalue si un coefficient de régression partiel diffère d’une valeur hypothétique, les autres variables explicatives restant dans le modèle.
Pour la variable j : H₀ : βⱼ = 0 contre H₁ : βⱼ ≠ 0. Au seuil α = 0,05, rejetez H₀ si p < 0,05, tout en examinant la taille de l’effet et les intervalles de confiance.
summary(model)$coefficients
# Tester une valeur non nulle, par exemple H0 : beta_wt = -3
b <- coef(model)["wt"]
se <- summary(model)$coefficients["wt", "Std. Error"]
t_value <- (b - (-3)) / se
p_value <- 2 * pt(abs(t_value), df.residual(model), lower.tail = FALSE)
c(t = t_value, p = p_value)> summary(model)$coefficients
Estimate Std. Error t value Pr(>|t|)
(Intercept) 37.22727012 1.59878754 23.28469 2.565e-20
wt -3.87783074 0.63273349 -6.12870 1.120e-06
hp -0.03177295 0.00902971 -3.51871 1.451e-03
> c(t = t_value, p = p_value)
t p
-1.38736 0.17587Pour le poids, t = −6,129 et p = 1,12×10⁻⁶ ; pour la puissance, t = −3,519 et p = 0,00145. Les deux valeurs p sont inférieures à 0,05, ce qui indique que chacune des associations partielles diffère de zéro dans ce modèle.
Inférence statistique
Calculer et interpréter les intervalles de confiance
Un intervalle de confiance indique une plage de valeurs du coefficient compatibles avec le modèle et l’incertitude d’échantillonnage.
confint(model, level = 0.95)
# Intervalle manuel à 95 % pour le coefficient du poids
b <- coef(model)["wt"]
se <- summary(model)$coefficients["wt", "Std. Error"]
b + c(-1, 1) * qt(0.975, df.residual(model)) * se> confint(model, level = 0.95)
2.5 % 97.5 %
(Intercept) 33.95738245 40.49715778
wt -5.17191604 -2.58374544
hp -0.05024078 -0.01330512L’intervalle à 95 % est [−5,172 ; −2,584] pour le poids et [−0,0502 ; −0,0133] pour la puissance. Aucun ne contient zéro, ce qui concorde avec les tests bilatéraux au seuil de 5 %. Pour le poids, les associations conditionnelles plausibles correspondent à une baisse d’environ 2,58 à 5,17 mpg par tranche de 1 000 livres.
Inférence statistique
Évaluer la significativité globale avec le test F
Le test F global compare le modèle estimé à un modèle avec constante seule et détermine si les variables explicatives améliorent conjointement la part de variation expliquée.
H₀ : βwt = βhp = 0 contre H₁ : au moins une pente est différente de zéro.
summary(model)$fstatistic
null_model <- lm(mpg ~ 1, data = mtcars)
anova(null_model, model)Analysis of Variance Table
Model 1: mpg ~ 1
Model 2: mpg ~ wt + hp
Res.Df RSS Df Sum of Sq F Pr(>F)
1 31 1126.05
2 29 195.05 2 931 69.21 9.109e-12 ***F(2, 29) = 69,21 et p = 9,11×10⁻¹². Rejetez l’hypothèse nulle conjointe : le poids et la puissance améliorent ensemble l’explication de mpg par rapport au modèle avec constante seule.
Inférence statistique
Effectuer des tests F partiels et conjoints
Les tests de modèles emboîtés évaluent si un groupe de variables apporte une information explicative supplémentaire au-delà des variables déjà incluses dans le modèle restreint.
Dans cet exemple, H₀ : les coefficients de hp et am sont conjointement nuls après contrôle de wt.
restricted <- lm(mpg ~ wt, data = mtcars)
unrestricted <- lm(mpg ~ wt + hp + am, data = mtcars)
anova(restricted, unrestricted)
# Hypothèses linéaires générales équivalentes
install.packages("car") # Exécuter une seule fois
library(car)
linearHypothesis(unrestricted, c("hp = 0", "am = 0"))Analysis of Variance Table
Model 1: mpg ~ wt
Model 2: mpg ~ wt + hp + am
Res.Df RSS Df Sum of Sq F Pr(>F)
1 30 278.32
2 28 180.29 2 98.034 7.613 0.00225 **
Linear hypothesis test:
hp = 0
am = 0
F = 7.613, Df = (2, 28), p-value = 0.00225Si la valeur p du test est inférieure au seuil choisi, rejetez les restrictions : le groupe testé apporte une contribution conjointe. Si elle est élevée, l’échantillon ne fournit pas assez d’éléments pour conclure que le modèle élargi améliore l’ajustement.
Inférence statistique
Évaluer l’ajustement sans surinterpréter le R²
Le R² mesure la proportion de variation de la variable expliquée attribuable aux valeurs ajustées dans l’échantillon. Le R² ajusté pénalise l’ajout de variables peu contributives.
summary(model)$r.squared
summary(model)$adj.r.squared
sigma(model)
actual_vs_fitted <- data.frame(
actual = mtcars$mpg,
fitted = fitted(model),
residual = residuals(model)
)> summary(model)$r.squared
[1] 0.8267855
> summary(model)$adj.r.squared
[1] 0.8148396
> sigma(model)
[1] 2.593412R² = 0,8268 et R² ajusté = 0,8148 : le modèle explique environ 82,7 % de la variation de mpg dans l’échantillon. L’écart-type résiduel de 2,593 mpg estime l’ordre de grandeur de l’erreur conditionnelle.
Diagnostics du modèle
Vérifier la linéarité et la forme fonctionnelle
Le graphique des résidus en fonction des valeurs ajustées devrait montrer un nuage sans structure particulière autour de zéro. Une courbure suggère des termes non linéaires omis ou une forme fonctionnelle inadaptée.
RESET de Ramsey : H₀ : la spécification linéaire est adéquate, contre une alternative faisant intervenir des puissances des valeurs ajustées.
plot(model, which = 1)
abline(h = 0, lty = 2, col = "#dc2626")
install.packages("lmtest") # Exécuter une seule fois
library(lmtest)
resettest(model, power = 2:3, type = "fitted")
# Exemple d’extension non linéaire
quadratic_model <- lm(mpg ~ wt + I(wt^2) + hp, data = mtcars)Sortie graphique : les résidus sont représentés autour de la ligne
horizontale zéro en fonction des valeurs ajustées. Une courbe lissée
présentant une courbure suggère une forme linéaire incomplète.
> resettest(model, power = 2:3, type = "fitted")
RESET test
RESET = 0.834, df1 = 2, df2 = 27, p-value = 0.445Sur le graphique, recherchez des courbes systématiques plutôt que des points isolés. Pour RESET, une petite valeur p contredit la forme fonctionnelle retenue ; une valeur élevée signifie que le test n’a pas détecté cette forme particulière de mauvaise spécification.
Diagnostics du modèle
Examiner la normalité des résidus
La normalité intervient surtout dans la validité exacte des tests t et F en petit échantillon. La moyenne conditionnelle peut rester linéaire même si les résidus ne suivent pas une loi normale.
Shapiro–Wilk : H₀ : les résidus suivent une loi normale, contre H₁ : ils ne suivent pas une loi normale.
plot(model, which = 2) # Graphique quantile-quantile normal
hist(residuals(model), breaks = 8, col = "#93c5fd",
main = "Résidus du modèle", xlab = "Résidu")
shapiro.test(residuals(model))> shapiro.test(residuals(model))
Shapiro-Wilk normality test
W = 0.9279, p-value = 0.0337
Sortie graphique : histogramme des résidus et graphique quantile-quantile normal.
La petite valeur p signale un écart détectable à la normalité.Sur le graphique quantile-quantile, un alignement global sur la droite est compatible avec une normalité approximative ; les écarts dans les queues révèlent une asymétrie ou des queues épaisses. Pour Shapiro–Wilk, p < 0,05 contredit la normalité, tandis que p ≥ 0,05 ne la prouve pas.
Diagnostics du modèle
Tester la constance de la variance des erreurs
L’hétéroscédasticité signifie que la variance conditionnelle des résidus varie avec les variables explicatives ou les valeurs ajustées. Sous exogénéité, les coefficients des MCO peuvent rester sans biais, mais les erreurs-types classiques peuvent être incorrectes.
Breusch–Pagan : H₀ : la variance des erreurs est constante, contre H₁ : elle dépend des régresseurs.
plot(model, which = 1)
plot(model, which = 3) # Graphique de dispersion-localisation
library(lmtest)
bptest(model)
# Test auxiliaire plus souple de type White
bptest(model, ~ fitted(model) + I(fitted(model)^2))> bptest(model)
studentized Breusch-Pagan test
BP = 0.8807, df = 2, p-value = 0.6438
Sortie graphique : résidus selon les valeurs ajustées et dispersion-localisation.
Ici, le test ne rejette pas la variance constante au seuil de 5 %.Une forme en éventail ou en entonnoir suggère une variance variable. Pour Breusch–Pagan, p < 0,05 indique une hétéroscédasticité ; p ≥ 0,05 signifie que le test n’a pas trouvé suffisamment d’éléments contre la constance de la variance.
Diagnostics du modèle
Tester l’autocorrélation lorsque les observations sont ordonnées
La corrélation sérielle concerne principalement les observations ordonnées dans le temps ou l’espace, ou regroupées en grappes. Elle viole l’indépendance et peut fausser les erreurs-types classiques.
Durbin–Watson : H₀ : absence d’autocorrélation des résidus d’ordre un. L’alternative bilatérale usuelle est une autocorrélation d’ordre un non nulle.
library(lmtest)
dwtest(model)
# Vérification visuelle lorsque l’ordre a un sens concret
plot(residuals(model), type = "o", pch = 19,
ylab = "Résidu", xlab = "Ordre des observations")
abline(h = 0, lty = 2)> dwtest(model)
Durbin-Watson test
DW = 1.3624, p-value = 0.0714
alternative hypothesis: true autocorrelation is greater than 0
Sortie graphique : résidus reliés selon l’ordre actuel des lignes.
Pour mtcars, cet ordre est illustratif et ne correspond pas au temps.Une statistique de Durbin–Watson proche de 2 est compatible avec une faible autocorrélation d’ordre un. Des valeurs nettement inférieures à 2 suggèrent une corrélation positive ; des valeurs supérieures à 2 suggèrent une corrélation négative. Utilisez la valeur p correspondant à l’alternative annoncée.
Diagnostics du modèle
Diagnostiquer la multicolinéarité avec le VIF
La multicolinéarité apparaît lorsque les variables explicatives apportent des informations redondantes. Elle accroît l’incertitude des coefficients et peut rendre les estimations instables sans nécessairement dégrader la prévision globale.
install.packages("car") # Exécuter une seule fois
library(car)
vif(model)
# VIF manuel pour wt
auxiliary <- lm(wt ~ hp, data = mtcars)
1 / (1 - summary(auxiliary)$r.squared)> vif(model)
wt hp
1.766625 1.766625
> 1 / (1 - summary(auxiliary)$r.squared)
[1] 1.766625
Les deux VIF sont bien inférieurs aux seuils de repérage usuels de 5 ou 10.Le VIF vaut 1 lorsqu’une variable explicative n’est pas liée aux autres. Des valeurs supérieures à 5 méritent un examen et des valeurs supérieures à 10 sont souvent considérées comme préoccupantes ; il s’agit toutefois de conventions de diagnostic, pas de seuils universels de test.
Diagnostics du modèle
Repérer les effets de levier, les valeurs atypiques et les observations influentes
Le levier mesure les combinaisons inhabituelles de variables explicatives, les résidus studentisés signalent les valeurs expliquées inhabituelles conditionnellement aux régresseurs, et la distance de Cook résume l’influence sur le modèle estimé.
influence_table <- data.frame(
car = rownames(mtcars),
leverage = hatvalues(model),
studentized = rstudent(model),
cooks_d = cooks.distance(model)
)
influence_table[order(-influence_table$cooks_d), ][1:5, ]
plot(model, which = 4) # Distance de Cook
plot(model, which = 5) # Résidus en fonction du levier> influence_table[order(-influence_table$cooks_d), ][1:5, ]
car leverage studentized cooks_d
17 Chrysler Imperial 0.0764 2.572 0.4236
20 Toyota Corolla 0.1067 2.129 0.2087
31 Maserati Bora 0.3942 0.929 0.1574
28 Lotus Europa 0.1942 1.413 0.1295
16 Lincoln Continental 0.1097 -1.514 0.0830
Sortie graphique : distance de Cook et résidus en fonction du levier.Dans ce modèle, Chrysler Imperial présente la plus grande distance de Cook (environ 0,424). Les règles D de Cook > 4/n ou levier > 2p/n aident au repérage : elles désignent des observations à examiner, pas à supprimer automatiquement.
Application
Utiliser des erreurs-types robustes à l’hétéroscédasticité
Les estimateurs robustes de covariance conservent les coefficients des MCO, mais ajustent leur incertitude estimée lorsque la constance de la variance est douteuse.
install.packages(c("sandwich", "lmtest")) # Exécuter une seule fois
library(sandwich)
library(lmtest)
coeftest(model, vcov. = vcovHC(model, type = "HC3"))
# Intervalles de confiance robustes
robust_vcov <- vcovHC(model, type = "HC3")
robust_se <- sqrt(diag(robust_vcov))
critical <- qt(0.975, df.residual(model))
cbind(estimate = coef(model),
lower = coef(model) - critical * robust_se,
upper = coef(model) + critical * robust_se)> coeftest(model, vcov. = vcovHC(model, type = "HC3"))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 37.227270 2.229805 16.6953 < 2.2e-16
wt -3.877831 0.768519 -5.0458 2.23e-05
hp -0.031773 0.009385 -3.3855 0.00206
Les coefficients restent inchangés ; seule leur incertitude estimée change.Comparez les erreurs-types robustes et classiques. Des différences importantes montrent que la formule homoscédastique influençait les résultats. Interprétez les coefficients dans les mêmes unités, mais fondez les tests et les intervalles sur la matrice de covariance robuste retenue.
Application
Comparer les modèles de façon raisonnée
Les tests F de modèles emboîtés évaluent des restrictions, tandis que le R² ajusté, l’AIC, le BIC et la validation examinent différents aspects de l’adéquation et de la parcimonie du modèle.
model_small <- lm(mpg ~ wt, data = mtcars)
model_full <- lm(mpg ~ wt + hp, data = mtcars)
anova(model_small, model_full) # Test F de modèles emboîtés
AIC(model_small, model_full)
BIC(model_small, model_full)
cbind(
adjusted_R2 = c(summary(model_small)$adj.r.squared,
summary(model_full)$adj.r.squared),
RMSE = c(sigma(model_small), sigma(model_full))
)> anova(model_small, model_full)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 30 278.32
2 29 195.05 1 83.27 12.381 0.001451 **
> AIC(model_small, model_full)
df AIC
model_small 3 166.029
model_full 4 156.652
> BIC(model_small, model_full)
df BIC
model_small 3 170.426
model_full 4 162.515Pour les modèles emboîtés, une petite valeur p de l’ANOVA soutient l’ajout des termes. Un AIC ou BIC plus faible indique un meilleur compromis entre ajustement et complexité parmi les modèles comparés. Ces critères ne prouvent ni la validité causale, ni la performance externe.
Application
Produire des prévisions et distinguer les intervalles
Un intervalle de confiance estime la réponse moyenne des voitures présentant les caractéristiques spécifiées. Un intervalle de prévision concerne une nouvelle voiture et est plus large, car il inclut la variation résiduelle individuelle.
new_car <- data.frame(wt = 3, hp = 120)
predict(model, newdata = new_car,
interval = "confidence", level = 0.95)
predict(model, newdata = new_car,
interval = "prediction", level = 0.95)> predict(model, newdata = new_car, interval = "confidence")
fit lwr upr
1 21.78102 20.77178 22.79027
> predict(model, newdata = new_car, interval = "prediction")
fit lwr upr
1 21.78102 16.38174 27.18031Pour wt = 3 et hp = 120, la prévision est de 21,78 mpg. L’intervalle de confiance à 95 % de la moyenne conditionnelle est [20,77 ; 22,79], tandis que l’intervalle de prévision à 95 % pour une nouvelle voiture est nettement plus large : [16,38 ; 27,18].
Application
Présenter les résultats avec leur portée statistique et concrète
Un compte rendu complet précise le modèle, l’échantillon, les unités, les estimations, leur incertitude, les diagnostics, les limites et l’interprétation exacte autorisée par le protocole de recherche.
# Tableau compact pour publication
install.packages("modelsummary") # Exécuter une seule fois
library(modelsummary)
modelsummary(model,
statistic = "conf.int",
stars = TRUE,
gof_map = c("nobs", "r.squared", "adj.r.squared", "rmse"))Model summary
────────────────────────────────────────
Estimate 95% CI
(Intercept) 37.227 [33.957, 40.497]
wt -3.878 [ -5.172, -2.584]
hp -0.032 [ -0.050, -0.013]
────────────────────────────────────────
Num.Obs. 32
R2 0.827
R2 Adj. 0.815
RMSE 2.593
────────────────────────────────────────Exemple : « À puissance constante, une hausse du poids de 1 000 livres était associée à une baisse de 3,88 mpg (IC à 95 % [−5,17 ; −2,58], p < 0,001). Le modèle expliquait 82,7 % de la variation de l’échantillon (R² ajusté = 0,815). » Précisez les limites des diagnostics et du protocole avant toute généralisation.