# R: Confused about robust standard errors using “felm” and “huxreg”

I’m working on panel data and want to use a “least squares dummy variable” (LSDV) approach, aka using factors to control for fixed effects rather than “within” differences.

I use a large set of indicators (many factor levels), so that normal OLS is relatively slow and the regression output is ugly. Therefore I explored the R-package lfe. It provides the function felm which “absorbs” factors (similar to Stats’s areg).

I need to use robust standard errors (HC1 or so) since tests indicate that there might be heteroscedasticity.

This is not so flamboyant after all. However, I obtain odd results for the robust SEs (using felm and huxreg).

Here is an example:

library(sandwich) # for robust SEs
library(lmtest)
library(huxtable) # for tables
library(lfe)      # "absorbe" FE Dummies with "felm"
library(ISLR)     # Auto data
library(jtools)   # for exporting summary


I) Let's compare felm and lm:

fe1 = felm(mpg~+cylinders+displacement | as.factor(year), data=Auto)
fe2 = lm(mpg~+cylinders+displacement + as.factor(year), data=Auto)
huxreg(fe1,fe2)


Output:

───────────────────────────────────────────────────
(1)            (2)
───────────────────────────────
cylinders            -0.297          -0.297
(0.359)         (0.359)
displacement         -0.045 ***      -0.045 ***
(0.006)         (0.006)
(Intercept)                          32.401 ***


So felm and lm yield the same results, just as expected.

As stated in the felm docs (https://www.rdocumentation.org/packages/lfe/versions/2.8-2/topics/felm), we can get robust SEs by calling:

summary(fe1, robust=T)


This is what I get:

Coefficients:
Estimate Robust s.e t value Pr(>|t|)
cylinders    -0.296565   0.350533  -0.846    0.398
displacement -0.045153   0.005994  -7.533  3.7e-13 ***


II) Let's get robust SEs:

Usually when I want to have robust SEs, I do something like this:

fe1_ro1 <- coeftest(fe1, vcov = vcovHC(fe1, type="HC1"))
fe2_ro1 <- coeftest(fe2, vcov = vcovHC(fe2, type="HC1"))
huxreg(fe1_ro1, fe2_ro1, number_format=6)


Output:

──────────────────────────────────────────────────
(1)            (2)
──────────────────────────────
cylinders           -0.296565    -0.296565
(4.593215)   (0.350533)
displacement        -0.045153    -0.045153 ***
(0.078539)   (0.005994)
(Intercept)                      32.400736 ***


Consider the (robust) SEs of displacement:

1. From summary(fe1, robust=T) = 0.005994,

2. From huxreg() using felm = 0.078539,

3. From huxreg() using lm = 0.005994.

I think 0.005994 is correct, but I cannot work out from where the other value (0.07) comes!

III) Try to pass summary() into regression table:

It appears that felm got the robust SEs right when using summary(fe1, robust=T). So I try to pass this result into a proper regression table.

I try export_summs() from jtools (https://www.rdocumentation.org/packages/jtools/versions/1.1.1/topics/export_summs):

export_summs(fe1, fe2, model.names = c("felm","lm"), robust = TRUE, number_format=6)


Output:

───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
felm                                       lm
─────────────────────────────────────────────────────────────────────────────────────
cylinders                                                           -0.296565                                 -0.296565
(0.358947)                                (0.359817)
displacement                                                        -0.045153 ***                             -0.045153 ***
(0.005891)                                (0.006158)
(Intercept)                                                                                                   32.400736 ***


It turnes out that the SEs for displacement using lm in the table above (0.006158) correspond to HC3 errors.

coeftest(fe2, vcov = vcovHC(fe2, type="HC3"))


So far, so good.

However, the SEs for displacement using felm (0.005891) are again different, and I do not see where they come from.

(Trying coeftest(fe1, vcov = vcovHC(fe1, type="HC3")) yields an error Error in 1 - diaghat : non-numeric argument to binary operator)

IV) So here are my questions:

• My stance is that robust SEs should be the same under felm and lm. Is this correct or do I miss some important argument here?
• If correct: How can I display correct robust SEs from felm in a proper (publication ready) regression table (in MS Word format) similar to huxreg?
• Finally and probably most important to me: Is there any alternative to felm (in R), viz. a command which allows to “absorb” a larger set of indicator variables to make LSDV workable?

Thanks!

## migrated from stackoverflow.comJan 27 at 15:36

This question came from our site for professional and enthusiast programmers.

1. The result from huxreg comes directly from the coeftest object:
> fe1_ro1

t test of coefficients:

Estimate Std. Error t value Pr(>|t|)
cylinders    -0.296565   4.593215 -0.0646   0.9486
displacement -0.045153   0.078539 -0.5749   0.5657


So the question is why coeftest is giving one result, and summary.felm another. The felm documentation doesn't say how exactly it calculates robust SEs; you might need to check the code or contact the author. I wouldn't assume that summary.felm is correct and coeftest is wrong (even if summary.felm gives significant results ;-) ). They may both be "right" but using different methods.

1. huxreg uses the tidy method from the broom package to get standard errors. felm has such a method but I don't think that summary.felm does. You can get round this by using the tidy_override function in recent huxtable.

2. You might look at biglm or speedglm for increasing OLS speed, though I don't think they specifically deal with LSDV issues. Another option would be plm, designed for panel data but perhaps useful for you.