Why does summary() show different standard errors than coeftest()?

Question

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!

jay.sf · Accepted Answer · 2021-12-26T11:33:47.870

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 **

Compute robust standard errors

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.

Publishing table with robust standard errors

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

Why does summary() show different standard errors than coeftest()?

1 Answers1

Compute robust standard errors

Publishing table with robust standard errors