0
$\begingroup$

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!

$\endgroup$

migrated from stackoverflow.com Jan 27 at 15:36

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

1
$\begingroup$
  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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.