Réplication de l'option «robuste» de Stata dans R

39

J'ai essayé de reproduire les résultats de l'option Stata robustdans R. J'ai utilisé la rlmcommande du package MASS ainsi que la commande lmrobdu package "robustbase". Dans les deux cas, les résultats sont assez différents de l’option "robuste" de Stata. Quelqu'un peut-il suggérer quelque chose dans ce contexte?

Voici les résultats que j'ai obtenus lorsque j'ai exécuté l'option robuste dans Stata:

. reg yb7 buildsqb7 no_bed no_bath rain_harv swim_pl pr_terrace, robust

Linear regression                                      Number of obs =    4451
                                                       F(  6,  4444) =  101.12
                                                       Prob > F      =  0.0000
                                                       R-squared     =  0.3682
                                                       Root MSE      =   .5721

------------------------------------------------------------------------------
             |               Robust
         yb7 |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
   buildsqb7 |   .0046285   .0026486     1.75   0.081    -.0005639     .009821
      no_bed |   .3633841   .0684804     5.31   0.000     .2291284    .4976398
     no_bath |   .0832654   .0706737     1.18   0.239    -.0552904    .2218211
   rain_harv |   .3337906   .0395113     8.45   0.000     .2563289    .4112524
     swim_pl |   .1627587   .0601765     2.70   0.007     .0447829    .2807346
  pr_terrace |   .0032754   .0178881     0.18   0.855    -.0317941    .0383449
       _cons |   13.68136   .0827174   165.40   0.000     13.51919    13.84353

Et voici ce que j’ai obtenu dans R avec l’option lmrob:

> modelb7<-lmrob(yb7~Buildsqb7+No_Bed+Rain_Harv+Swim_Pl+Gym+Pr_Terrace, data<-bang7)
> summary(modelb7)

Call:
lmrob(formula = yb7 ~ Buildsqb7 + No_Bed + Rain_Harv + Swim_Pl + Gym + Pr_Terrace, 
    data = data <- bang7)
 \--> method = "MM"
Residuals:
      Min        1Q    Median        3Q       Max 
-51.03802  -0.12240   0.02088   0.18199   8.96699 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 12.648261   0.055078 229.641   <2e-16 ***
Buildsqb7    0.060857   0.002050  29.693   <2e-16 ***
No_Bed       0.005629   0.019797   0.284   0.7762    
Rain_Harv    0.230816   0.018290  12.620   <2e-16 ***
Swim_Pl      0.065199   0.028121   2.319   0.0205 *  
Gym          0.023024   0.014655   1.571   0.1162    
Pr_Terrace   0.015045   0.013951   1.078   0.2809    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

Robust residual standard error: 0.1678 
Multiple R-squared:  0.8062,    Adjusted R-squared:  0.8059 
utilisateur56579
la source
3
Bienvenue sur Cross Validated! J'ai rendu votre titre un peu plus descriptif et ajouté un peu de formatage. En général, les questions de programmation ne sont pas sur le sujet ici, mais je pense que le vôtre, parce que cela implique des problèmes statistiques. En espérant vous voir ...
Matt Krause Le
3
Cela aiderait énormément si vous avez au moins collé le code utilisé pour estimer les modèles dans Stata et R (encore mieux si vous fournissez un exemple entièrement reproductible). Lorsque vous dites "les résultats diffèrent" - si vous estimez le même modèle, seules les erreurs types devraient différer, pas les estimations de coefficients.
Andy W
d'accord ... Voici les résultats que j'ai obtenus avec l'option robuste de STATA:
user56579 le
5
ressemble lmrobn'est pas la même chose que reg y x, robust. Google "Erreurs-types R compatibles avec l'hétéroscédasticité". Vous obtiendrez des pages vous expliquant comment utiliser les bibliothèques lmtestet sandwich.
generic_user
3
Stata utilise un petit facteur de correction de n / (nk). R fait habituellement quelque chose d'autre, alors assurez-vous de l'ajuster.
Dimitriy V. Masterov

Réponses:

47

Charles y est presque dans sa réponse, mais l' robustoption de la regresscommande (et d'autres commandes d'estimation de régression) dans Stata permet d'utiliser plusieurs types d'hétéroscédasticité et d'auto-corrélation de puissants estimateurs matriciels à variance-covariance robustes, de même que la coeftestfonction du lmtestpaquet son tour dépend des matrices de variance-covariance respectives produites par la vcovHCfonction dans lasandwich package.

Cependant, les matrices de variance-covariance par défaut utilisées par les deux sont différentes:
1. La matrice de variance-covariance par défaut renvoyée vcocHCest ainsi nommée HC3pour les raisons décrites dans la page de manuel de vcovHC.
2. L’ sandwichoption utilisée par Charles fait coeftestappel à la HC0robuste matrice de variance-covariance.
3. Pour reproduire le comportement par défaut de Stata consistant à utiliser l' robustoption dans un appel à, regressvous devez demander vcovHCà utiliser la HC1matrice de variance-covariance robuste.

En savoir plus à ce sujet ici .

L'exemple suivant qui illustre tous les points évoqués ci-dessus est basé sur l'exemple ci- dessous .

library(foreign)
library(sandwich)
library(lmtest)

dfAPI = read.dta("http://www.ats.ucla.edu/stat/stata/webbooks/reg/elemapi2.dta")
lmAPI = lm(api00 ~ acs_k3 + acs_46 + full + enroll, data= dfAPI)
summary(lmAPI)                                  # non-robust

# check that "sandwich" returns HC0
coeftest(lmAPI, vcov = sandwich)                # robust; sandwich
coeftest(lmAPI, vcov = vcovHC(lmAPI, "HC0"))    # robust; HC0 

# check that the default robust var-cov matrix is HC3
coeftest(lmAPI, vcov = vcovHC(lmAPI))           # robust; HC3 
coeftest(lmAPI, vcov = vcovHC(lmAPI, "HC3"))    # robust; HC3 (default)

# reproduce the Stata default
coeftest(lmAPI, vcov = vcovHC(lmAPI, "HC1"))    # robust; HC1 (Stata default)

La dernière ligne de code ci-dessus reproduit les résultats de Stata:

use http://www.ats.ucla.edu/stat/stata/webbooks/reg/elemapi2
regress api00 acs_k3 acs_46 full enroll, robust
tchakravarty
la source
Le lien vers les données est mort. Pouvez-vous s'il vous plaît mettre à jour le lien? S'agit-il du même fichier: faculty.smu.edu/tfomby/eco5350/data/Examples/elemapi2.dta ?
vasili111
Comment reproduire aussi les intervalles de confiance?
vasili111
29

J'ai trouvé une description sur le site Web suivant qui reproduit l'option «robuste» de Stata dans R.

https://economictheoryblog.com/2016/08/08/robust-standard-errors-in-r

En suivant les instructions, il vous suffit de charger une fonction dans votre session R, puis de définir le paramètre '' robuste '' de votre fonction de résumé sur TRUE.

summary(lm.object, robust=TRUE)
Alex Rato
la source
10

À partir d'avril 2018, je pense que vous voulez le estimatrpaquet , qui fournit un remplacement presque instantané pour lm. Plusieurs exemples tirés presque de la documentation:

library(estimatr)
library(car)

# HC1 robust standard errors
model <- lm_robust(GPA_year2 ~ gpa0 + ssp, data = alo_star_men,
                   se_type = "stata")
summary(model)
#> 
#> Call:
#> lm_robust(formula = GPA_year2 ~ gpa0 + ssp, data = alo_star_men, 
#>     se_type = "stata")
#> 
#> Standard error type:  HC1 
#> 
#> Coefficients:
#>             Estimate Std. Error  Pr(>|t|) CI Lower CI Upper  DF
#> (Intercept) -3.60625    1.60084 0.0258665 -6.77180  -0.4407 137
#> gpa0         0.06814    0.02024 0.0009868  0.02812   0.1082 137
#> ssp          0.31917    0.18202 0.0817589 -0.04077   0.6791 137
#> 
#> Multiple R-squared:  0.09262 ,   Adjusted R-squared:  0.07937 
#> F-statistic: 6.992 on 2 and 137 DF,  p-value: 0.001284

# HC1 cluster robust standard errors
model2 <- lm_robust(GPA_year2 ~ gpa0 + ssp, cluster = ssp,
                   data = alo_star_men, se_type = "stata")
summary(model2)
#> 
#> Call:
#> lm_robust(formula = GPA_year2 ~ gpa0 + ssp, data = alo_star_men, 
#>     clusters = ssp, se_type = "stata")
#> 
#> Standard error type:  stata 
#> 
#> Coefficients:
#>             Estimate Std. Error Pr(>|t|) CI Lower CI Upper DF
#> (Intercept) -3.60625   1.433195 0.240821 -21.8167  14.6042  1
#> gpa0         0.06814   0.018122 0.165482  -0.1621   0.2984  1
#> ssp          0.31917   0.004768 0.009509   0.2586   0.3798  1
#> 
#> Multiple R-squared:  0.09262 ,   Adjusted R-squared:  0.07937 
#> F-statistic: 6.992 on 2 and 137 DF,  p-value: 0.001284

Le carpackage facilite ensuite la réalisation de tests d’hypothèses omnibus pour ces modèles:

linearHypothesis(model, c("gpa0 = ssp"))
#> Linear hypothesis test
#> 
#> Hypothesis:
#> gpa0 - ssp = 0
#> 
#> Model 1: restricted model
#> Model 2: GPA_year2 ~ gpa0 + ssp
#> 
#>   Res.Df Df  Chisq Pr(>Chisq)
#> 1    138                     
#> 2    137  1 1.8859     0.1697
alexpghayes
la source
4

Je modifierais la question. Vous confondez la régression robuste avec la commande robuste de Stata. Il ne semble y avoir aucun avantage à introduire cette confusion.

Je pense qu'il y a quelques approches. Je ne les ai pas tous regardés et je ne sais pas quel est le meilleur:

Le paquet sandwich:

library(sandwich)    
coeftest(model, vcov=sandwich)

Mais cela ne me donne pas les mêmes réponses que Stata pour une raison quelconque. Je n'ai jamais essayé de comprendre pourquoi - mais ci-dessus dans les commentaires, une réponse est suggérée - je n'utilise tout simplement pas ce paquet.

Le paquetage rms:

Je trouve cela un peu pénible, mais j'obtiens généralement de bonnes réponses avec quelques efforts. Et c'est le plus utile pour moi.

model = ols(a~b, x=TRUE)    
robcov(model)

Vous pouvez le coder à partir de zéro

Voir cet article de blog ( http://thetarzan.wordpress.com/2011/05/28/heteroskedasticity-robust-and-clustered-standard-errard-errors-in-r/ ). Cela semble être l'option la plus douloureuse, mais remarquablement facile, et cette option fonctionne souvent mieux.

Charles
la source
4
Charles a raison sur le point principal, mais pour expliciter ce qui est implicite ailleurs, notez que Stata n'a aucune robustcommande! (Il y a une commande de programmeur _robust, non directement pertinente ici.) Pour obtenir des erreurs standard robustes (Huber-Eicker-White-sandwich), l'approche moderne de Stata consiste à spécifier vce(robust)en option. L'ancienne approche pour spécifier une robustoption fonctionne toujours. Plus généralement, la confusion causée par la différence entre une régression robuste (etc.) et des SE "robustes" est regrettable.
Nick Cox le
Hey. Merci beaucoup. Les codes fonctionnent et fournissent effectivement les résultats obtenus par Stata. Juste une question. Je comprends que la régression robuste est différente des erreurs-types robustes et qu’elle est utilisée lorsque vos données contiennent des valeurs aberrantes. Mais cela résout également le problème de l'hétéroscédasticité. Quelqu'un pourrait-il me dire si l’estimation de type MM fournie par la commande "lmrob" du paquet "robustbase" peut être utilisée simultanément pour résoudre le problème des valeurs aberrantes et de l’hétéroscédasticité?
user56579
@ user56579 Mon hypothèse est que vous souhaitiez poser une question distincte à ce sujet.
tchakravarty