¿Cómo cambiar errores estandards normales a robustos en un modelo lineal?

Usemos la data de gapminder:

library(rio)
data<-import("https://raw.githubusercontent.com/resbaz/r-novice-gapminder-files/master/data/gapminder-FiveYearData.csv")

Hagamos una regresión simple con expectativa de vida como dependiente y gdpPercap y población (en logaritmo) como independiente

m1=lm(lifeExp~gdpPercap+log(pop), data=data)

Ahora pedimos los resultados con summary

summary(m1)

## 
## Call:
## lm(formula = lifeExp ~ gdpPercap + log(pop), data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -77.315  -7.558   1.479   7.803  20.199 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.121e+01  2.454e+00  12.719   <2e-16 ***
## gdpPercap   7.612e-04  2.516e-05  30.252   <2e-16 ***
## log(pop)    1.445e+00  1.546e-01   9.345   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.23 on 1701 degrees of freedom
## Multiple R-squared:  0.3729,	Adjusted R-squared:  0.3722 
## F-statistic: 505.8 on 2 and 1701 DF,  p-value: < 2.2e-16

El resumen muestra los errores normales. Para estimar los errores robustos (robust standard errors) podemos usar el paquete sandwich y lmtest

Con la función vcovHC podemos obtener la matriz de varianzas y covarianzas.

library(sandwich)
vcovHC(m1, type = "HC0") 

##               (Intercept)     gdpPercap      log(pop)
## (Intercept)  6.5237234376  1.090675e-04 -4.477830e-01
## gdpPercap    0.0001090675  8.181630e-09 -1.024504e-05
## log(pop)    -0.4477830082 -1.024504e-05  3.213692e-02

Y si complementamos eso con coeftest podemos probar si los coeficientes son significativos usando los errores standard robustos.

library(lmtest)
coef2<-coeftest(m1, vcov = vcovHC(m1, "HC0")) 

Remplacemos el resumen del modelo usando estos errores estandards.

sum1<-summary(m1)
sum1$coefficients[,2]<-coef2[,2]

Comparemos:

summary(m1)

## 
## Call:
## lm(formula = lifeExp ~ gdpPercap + log(pop), data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -77.315  -7.558   1.479   7.803  20.199 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.121e+01  2.454e+00  12.719   <2e-16 ***
## gdpPercap   7.612e-04  2.516e-05  30.252   <2e-16 ***
## log(pop)    1.445e+00  1.546e-01   9.345   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.23 on 1701 degrees of freedom
## Multiple R-squared:  0.3729,	Adjusted R-squared:  0.3722 
## F-statistic: 505.8 on 2 and 1701 DF,  p-value: < 2.2e-16

sum1

## 
## Call:
## lm(formula = lifeExp ~ gdpPercap + log(pop), data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -77.315  -7.558   1.479   7.803  20.199 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.121e+01  2.554e+00  12.719   <2e-16 ***
## gdpPercap   7.612e-04  9.045e-05  30.252   <2e-16 ***
## log(pop)    1.445e+00  1.793e-01   9.345   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.23 on 1701 degrees of freedom
## Multiple R-squared:  0.3729,	Adjusted R-squared:  0.3722 
## F-statistic: 505.8 on 2 and 1701 DF,  p-value: < 2.2e-16

Tengo un tutorial de R en español. Lo puede revisar aquí: http://www.joseincio.com/R_tutorial/