Comment calculer les erreurs standard robustes dans R

L’une des hypothèses de la régression linéaire est que les résidus du modèle sont également dispersés à chaque niveau de la variable prédictive.

Lorsque cette hypothèse n’est pas respectée, on dit que l’hétéroscédasticité est présente dans un modèle de régression.

Lorsque cela se produit, les erreurs types des coefficients de régression du modèle deviennent peu fiables.

Pour tenir compte de cela, nous pouvons calculer des erreurs types robustes , qui sont « robustes » contre l’hétéroscédasticité et peuvent nous donner une meilleure idée des véritables valeurs d’erreur type pour les coefficients de régression.

L’exemple suivant montre comment calculer des erreurs types robustes pour un modèle de régression dans R.

Exemple : calcul des erreurs types robustes dans R

Supposons que nous ayons le bloc de données suivant dans R qui contient des informations sur les heures étudiées et les résultats des examens obtenus par 20 étudiants dans une classe :

#create data frame
df <- data.frame(hours=c(1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4,
                         4, 5, 5, 5, 6, 6, 7, 7, 8),
                 score=c(67, 68, 74, 70, 71, 75, 80, 70, 84, 72,
                         88, 75, 95, 75, 99, 78, 99, 65, 96, 70))

#view head of data frame
head(df)

  hours score
1     1    67
2     1    68
3     1    74
4     1    70
5     2    71
6     2    75

Nous pouvons utiliser la fonction lm() pour ajuster un modèle de régression dans R qui utilise les heures comme variable prédictive et le score comme variable de réponse :

#fit regression model
fit <- lm(score ~ hours, data=df)

#view summary of model
summary(fit)

Call:
lm(formula = score ~ hours, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-19.775  -5.298  -3.521   7.520  18.116 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   71.158      4.708   15.11 1.14e-11 ***
hours          1.945      1.075    1.81    0.087 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.48 on 18 degrees of freedom
Multiple R-squared:  0.154,	Adjusted R-squared:  0.107 
F-statistic: 3.278 on 1 and 18 DF,  p-value: 0.08696

Le moyen le plus simple de vérifier visuellement si l’hétéroscédasticité est un problème dans le modèle de régression est de créer un tracé résiduel :

#create residual vs. fitted plot
plot(fitted(fit), resid(fit))

#add a horizontal line at y=0 
abline(0,0)

L’axe des x montre les valeurs ajustées de la variable de réponse et l’axe des y montre les résidus correspondants.

Sur le graphique, nous pouvons voir que la variance des résidus augmente à mesure que les valeurs ajustées augmentent.

Cela indique que l’hétéroscédasticité est probablement un problème dans le modèle de régression et que les erreurs types du résumé du modèle ne sont pas fiables.

Pour calculer des erreurs standard robustes, nous pouvons utiliser la fonction coeftest() du package lmtest et la fonction vcovHC() du package sandwich comme suit :

library(lmtest)
library(sandwich)

#calculate robust standard errors for model coefficients
coeftest(fit, vcov = vcovHC(fit, type = 'HC0'))

t test of coefficients:

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept)  71.1576     3.3072 21.5160 2.719e-14 ***
hours         1.9454     1.2072  1.6115    0.1245    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Notez que l’erreur standard pour la variable prédictive des heures est passée de 1,075 dans le résumé du modèle précédent à 1,2072 dans ce résumé du modèle.

Étant donné que l’hétéroscédasticité est présente dans le modèle de régression d’origine, cette estimation de l’erreur type est plus fiable et doit être utilisée lors du calcul d’un intervalle de confiance pour la variable prédictive des heures .

Remarque : Le type d’estimation le plus courant à calculer dans la fonction vcovHC() est ‘HC0’, mais vous pouvez vous référer à la documentation pour trouver d’autres types d’estimation.

Ressources additionnelles

Les didacticiels suivants expliquent comment effectuer d’autres tâches courantes dans R :

Comment effectuer le test de White pour l’hétéroscédasticité dans R
Comment interpréter la sortie de régression linéaire dans R
Comment créer un tracé résiduel dans R