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

Question

He tratado de replicar los resultados de la opción de Stata robust en R. He utilizado el rlm del paquete MASS y también el comando lmrob del paquete "robustbase". En ambos casos, los resultados son bastante diferentes de la opción "robusta" de Stata. ¿Puede alguien sugerir algo en este contexto?

Estos son los resultados que obtuve al ejecutar 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 he obtenido 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

Answer 1

Charles está casi ahí en su respuesta, pero robust opción del regress (y otros comandos de estimación de la regresión) en Stata permite utilizar múltiples tipos de estimadores de la matriz de varianza-covarianza robustos a la heteroscedasticidad y a la autocorrelación, al igual que el comando coeftest en la función lmtest que a su vez depende de las respectivas matrices de varianza-covarianza producidas por el paquete vcovHC en la función sandwich paquete.

Sin embargo, las matrices de varianza-covarianza utilizadas por defecto por ambos son diferentes:
1. La matriz de varianza-covarianza por defecto devuelta por vcocHC es el llamado HC3 por las razones descritas en la página de manual de vcovHC .
2. El sandwich opción utilizada por Charles hace coeftest utilizar el HC0 matriz de varianza-covarianza robusta.
3. Para reproducir el comportamiento por defecto de Stata de utilizar el robust en una llamada a regress debe solicitar vcovHC para utilizar el HC1 matriz de varianza-covarianza robusta.

Más información ici .

El siguiente ejemplo que demuestra todos los puntos expuestos anteriormente se basa en el ejemplo ici .

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

Answer 2

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 tiene que hacer es cargar una función en su sesión de R y luego establecer el parámetro ''robusto'' en su función de resumen a TRUE.

summary(lm.object, robust=TRUE)

Answer 3

A partir de abril de 2018 creo que quiere el estimatr paquete que proporciona un reemplazo casi directo para lm . Varios ejemplos sacados 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 car facilita la realización de pruebas de hipótesis globales 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

Answer 4

Yo editaría la pregunta. Estás confundiendo la regresión robusta con el comando robusto de Stata. No parece haber ningún beneficio en introducir esta confusión.

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

El paquete de sándwiches:

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

Pero esto no me da las mismas respuestas que obtengo 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:

A mí me resulta un poco pesado trabajar con esto, pero suelo obtener buenas respuestas con algo de esfuerzo. Y es el más útil para mí.

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

Puedes codificarla desde cero

Véase esta entrada del blog ( http://thetarzan.wordpress.com/2011/05/28/heteroskedasticity-robust-and-clustered-standard-errors-in-r/ ). Parece la opción más dolorosa, pero es notablemente fácil y esta opción suele ser la que mejor funciona.