I run a glm() using robust standard errors. For a subsequent model comparison I calculate the difference of two regression models (coefficients & se). For that calculation I use the summary() function. However, the summary function of the models show different standard errors than the ones I get from coeftest(). Values for coefficients remain identical.
Input:
mod.01 <- glm(dep ~ indep1 + indep2 + indep3,
family = binomial (link = "logit"), data = data)
coeftest(mod.01, vcov. = vcovHC, type= "HC3", df = NULL)
summary(mod.01, robust=T)
Output:
coeftest()
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.72917626 0.16367787 -16.6741 < 2.2e-16 ***
indep1 0.00427870 0.41928906 0.0102 0.991859
indep2 2.00243724 0.19757861 10.1349 < 2.2e-16 ***
indep3 0.36385098 0.32783817 1.1098 0.267159
summary()
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.7291763 0.1758744 -15.518 < 2e-16 ***
indep1 0.0042787 0.3389472 0.013 0.98993
indep2 2.0024372 0.1746829 11.463 < 2e-16 ***
indep3 0.3638510 0.2604196 1.397 0.16236
How do I get the robust se of mod.01 into the summary function? So that eventually my subsequent calculation includes the correct robust se will also be shown shown in the regression table.
Thanks in advance!
The merits of lmtest::coeftest is that it is possible to use a different covariance matrix than computed by lm().
fit <- glm(am ~ hp + cyl, family=binomial(link="logit"), mtcars)
summary(fit)
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) 5.83220 2.06595 2.823 0.00476 **
# hp 0.02775 0.01366 2.031 0.04228 *
# cyl -1.70306 0.60286 -2.825 0.00473 **
If you put the initial covariance matrix from lm() into coeftest, the standard errors are identical.
library(sandwich); library(lmtest)
coeftest(fit, vcov.=vcov(fit))
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) 5.832204 2.065949 2.8230 0.004757 **
# hp 0.027748 0.013664 2.0307 0.042282 *
# cyl -1.703064 0.602862 -2.8250 0.004729 **
Now, using coeftest we may specify a different covariance matrix. (Sidenote: The actual one with "robust" standard errors, also referred to as white standard errors, is specified with type='HC0', the other 'HC*' are refinements, e.g. Stata uses 'HC1').
(ct <- coeftest(fit, vcov.=vcovHC(fit, type='HC0')))
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) 5.832204 2.139307 2.7262 0.006407 **
# hp 0.027748 0.012254 2.2643 0.023553 *
# cyl -1.703064 0.572045 -2.9772 0.002909 **
That's how you get different standard errors from coeftest.
To put those corrected standard errors in a table, here's a way using texreg::screenreg. You may also us texreg::texreg for LaTeX or texreg::htmlreg for HTML, body stays basically the same.
Here with conventional standard errors,
texreg::screenreg(fit)
# =========================
# Model 1
# -------------------------
# (Intercept) 5.83 **
# (2.07)
# hp 0.03 *
# (0.01)
# cyl -1.70 **
# (0.60)
# -------------------------
# AIC 34.63
# BIC 39.02
# Log Likelihood -14.31
# Deviance 28.63
# Num. obs. 32
# =========================
# *** p < 0.001; ** p < 0.01; * p < 0.05
and here with robust standard errors, both, p-values and standard errors need to be overrided.
texreg::screenreg(fit,
override.pvalues=ct[, 4],
override.se=ct[, 3])
# =========================
# Model 1
# -------------------------
# (Intercept) 5.83 **
# (2.73)
# hp 0.03 *
# (2.26)
# cyl -1.70 **
# (-2.98)
# -------------------------
# AIC 34.63
# BIC 39.02
# Log Likelihood -14.31
# Deviance 28.63
# Num. obs. 32
# =========================
# *** p < 0.001; ** p < 0.01; * p < 0.05