r studio

Regresión con errores estándar corregidos de heteroskedasticidad

Question

Me gustaría encontrar la implementación de R que más se parezca a la salida de Stata para ajustar una función de Regresión de mínimos cuadrados con errores estándar corregidos heteroskedastic. Específicamente, me gustaría que los errores estándar corregidos estén en el "resumen" y no tengan que hacer cálculos adicionales para mi ronda inicial de pruebas de hipótesis. Estoy buscando una solución que sea tan "limpia" como lo que proporcionan Eviews y Stata.

Hasta ahora, usando el paquete "lmtest" lo mejor que puedo encontrar es:

model <- lm(...)
coeftest(model, vcov = hccm) 

Esto me da la salida que quiero, pero no parece estar usando "coeftest" para su propósito declarado. También tendría que usar el resumen con los errores estándar incorrectos para leer las estadísticas R ^ 2 y F, etc. Siento que debería existir una solución "de una línea" para este problema, dada la dinámica de R.

Gracias

Answer 1

Ahora hay una solución de una línea que utiliza lm_robust del paquete estimatr , que puede instalar desde los paquetes de instalación CRAN install.packages(estimatr) .

> library(estimatr)
> lmro <- lm_robust(mpg ~ hp, data = mtcars, se_type = "stata")
> summary(lmro)

Call:
lm_robust(formula = mpg ~ hp, data = mtcars, se_type = "stata")

Standard error type:  HC1 

Coefficients:
            Estimate Std. Error  Pr(>|t|) CI Lower CI Upper DF
(Intercept) 30.09886    2.07661 4.348e-15 25.85785 34.33987 30
hp          -0.06823    0.01356 2.132e-05 -0.09592 -0.04053 30

Multiple R-squared:  0.6024 ,   Adjusted R-squared:  0.5892 
F-statistic: 45.46 on 1 and 30 DF,  p-value: 1.788e-07

También puede obtener salida ordenada:

> tidy(lmro)
         term    estimate std.error      p.value    ci.lower
1 (Intercept) 30.09886054 2.0766149 4.347723e-15 25.85784704
2          hp -0.06822828 0.0135604 2.131785e-05 -0.09592231
     ci.upper df outcome
1 34.33987404 30     mpg
2 -0.04053425 30     mpg

Los errores estándar "stata" predeterminan a los errores estándar "HC1", que son los errores estándar de rob predeterminados en Stata. También puede obtener "classical", "HC0", "HC1", "HC2", "HC3" y también varios errores estándar agrupados (incluidos los que coinciden con Stata).

Answer 2

Creo que estás en el camino correcto con coeftest en el paquete lmtest. Eche un vistazo al paquete sándwich que incluye esta funcionalidad y está diseñado para trabajar de la mano con el paquete lmtest que ya ha encontrado.

> # generate linear regression relationship
> # with Homoskedastic variances
> x <- sin(1:100)
> y <- 1 + x + rnorm(100)
> ## model fit and HC3 covariance
> fm <- lm(y ~ x)
> vcovHC(fm)
            (Intercept)           x
(Intercept) 0.010809366 0.001209603
x           0.001209603 0.018353076
> coeftest(fm, vcov. = vcovHC)

t test of coefficients:

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept)  1.01973    0.10397  9.8081 3.159e-16 ***
x            0.93992    0.13547  6.9381 4.313e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Para obtener la prueba F, mire la función waldtest() :

> waldtest(fm, vcov = vcovHC)
Wald test

Model 1: y ~ x
Model 2: y ~ 1
  Res.Df Df      F    Pr(>F)    
1     98                        
2     99 -1 48.137 4.313e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Siempre puedes preparar una función simple para combinar estos dos para ti si quisieras el one-liner ...

Hay muchos ejemplos en la viñeta Econometric Computing con HC y HAC Covariance Matrix Estimators que viene con el paquete sandwich de lmtest y sandwich para hacer lo que quiera.

Edición: una sola línea podría ser tan simple como:

mySummary <- function(model, VCOV) {
    print(coeftest(model, vcov. = VCOV))
    print(waldtest(model, vcov = VCOV))
}

Que podemos usar de esta manera (en los ejemplos anteriores):

> mySummary(fm, vcovHC)

t test of coefficients:

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept)  1.01973    0.10397  9.8081 3.159e-16 ***
x            0.93992    0.13547  6.9381 4.313e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Wald test

Model 1: y ~ x
Model 2: y ~ 1
  Res.Df Df      F    Pr(>F)    
1     98                        
2     99 -1 48.137 4.313e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Answer 3

Encontré una función R que hace exactamente lo que estás buscando. Te da errores estándar robustos sin tener que hacer cálculos adicionales. Ejecuta summary() en un objeto lm.objetivo y si configura el parámetro robust=T le devuelve los errores estándar consistentes de heteroscedasticidad de tipo Stata.

summary(lm.object, robust=T)

Puede encontrar la función en https://economictheoryblog.com/2016/08/08/robust-standard-errors-in-r/