Replicando la opción "robusta" de Stata en R

He estado tratando de replicar los resultados de la opción Stata robusten R. He utilizado el rlmcomando del paquete MASS y también el comando lmrobdel paquete "robustbase". En ambos casos, los resultados son bastante diferentes de la opción "robusta" en Stata. ¿Alguien puede sugerir algo en este contexto?

Aquí están los resultados que obtuve cuando ejecuté la opción robusta en 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

Y esto es lo que obtuve en R con la opción 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 

Charles está casi allí en su respuesta, pero la robustopción del regresscomando (y otros comandos de estimación de regresión) en Stata hace posible el uso de múltiples tipos de estimadores de matriz de varianza-covarianza robustos de autocorrelación y autocorrelación, al igual que la coeftestfunción en el lmtestpaquete, que en turn depende de las respectivas matrices de varianza-covarianza producidas por la vcovHCfunción en el sandwichpaquete.

Sin embargo, las matrices de varianza-covarianza predeterminadas utilizadas por los dos son diferentes:
1. La matriz de varianza-covarianza predeterminada que devuelve vcocHCes la llamada HC3por las razones descritas en la página del manual para vcovHC.
2. La sandwichopción utilizada por Charles coeftestutiliza la HC0matriz robusta de varianza-covarianza.
3. Para reproducir el comportamiento predeterminado de Stata de usar la robustopción en una llamada regress, debe solicitar el vcovHCuso de la HC1matriz robusta de varianza-covarianza.

Lea más sobre esto aquí .

El siguiente ejemplo que demuestra todos los puntos mencionados anteriormente se basa en el ejemplo aquí .

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 última línea de código anterior reproduce los resultados de Stata:

use http://www.ats.ucla.edu/stat/stata/webbooks/reg/elemapi2
regress api00 acs_k3 acs_46 full enroll, robust

Encontré una descripción en el siguiente sitio web que replica la opción "robusta" de Stata en R.

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

Siguiendo las instrucciones, todo lo que necesita hacer es cargar una función en su sesión R y luego establecer el parámetro '' robusto '' en su función de resumen como VERDADERO.

summary(lm.object, robust=TRUE)

A partir de abril de 2018, creo que desea el estimatrpaquete , que proporciona una caída casi total del reemplazo lm. Varios ejemplos extraídos casi de la documentación:

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

El carpaquete de continuación, hace que sea fácil de realizar pruebas de hipótesis ómnibus para estos modelos:

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

Yo editaría la pregunta. Estás confundiendo una regresión robusta con el comando robusto de Stata. Parece que no hay ningún beneficio en introducir esta confusión.

Creo que hay algunos enfoques. No los he mirado a todos y no estoy seguro de cuál es el mejor:

El paquete de sándwich:

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

Pero esto no me da las mismas respuestas que recibo de Stata por alguna razón. Nunca he tratado de averiguar por qué, pero arriba en los comentarios hay una respuesta sugerida: simplemente no uso este paquete.

El paquete rms:

Me resulta un poco difícil trabajar con esto, pero generalmente obtengo buenas respuestas con algo de esfuerzo. Y es lo más útil para mí.

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

Puedes codificarlo desde cero

Consulte esta publicación de blog ( http://thetarzan.wordpress.com/2011/05/28/heteroskedasticity-robust-and-clustered-standard-errors-in-r/ ). Parece la opción más dolorosa, pero notablemente fácil y esta opción a menudo funciona mejor.