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

Question

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!

Answer 1

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.

2. 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.

3. 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.

Answer 2

This is an old thread, but I wanted to clarify some issues in case anyone stumbles upon it.

First an aside: You'd be much better off posting a question about R and R packages on StackOverflow than StackExchange.

1. lfe::felm() has its own idiomatic way of calculating robust (and clustered) standard errors; more on this below. You should use that instead of trying to pass it through lmtest::coeftest(), which does not provide a method for handling felm objects. It might seem like it's working because the code executes (at least, sometimes), but the coeftest() documentation makes it pretty clear that it is intended for specific objects.

2. Okay, so how to report robust SEs from felm models in a regression table? As the OP correctly states, huxtable uses broom::tidy() underneath the hood. The latest version of broom (v0.7.1 at the time of writing) actually includes some bug fixes and and enhancements specific to felm objects. You can read the gory details here and the ?tidy.felm help documentation provides a bunch of examples. The TL;DR version is to install the latest version of broom and then use the se.type argument.

Here's an example using the modelsummary package, which has similar functionality to huxreg, but I personally find a easier to control output and customization.

library(lmtest)
library(sandwich)
library(lfe)
library(ISLR)         ## Auto data
library(modelsummary) ## NB: make sure you've updated to broom>=v.0.7.1 first

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

mods = list(fe1, coeftest(fe2, vcov = vcovHC(fe2, type="HC1")))

msummary(
    mods, 
    output = 'table.docx', ## Optional if saving to MS Word as OP requested
    se.type = 'robust'     ## Argument being passed to tidy.felm()
    )

^{Created on 2020-10-05 by the reprex package (v0.3.0)}