How to have both lm and fixest model with different vcov standard errors in modelsummary package?

user_14772152 0
Viewed 58

I'd like to create an LM with HC3 corrected standard errors and a fixest model with cluster robust standard errors in the same table.

see my MRE below:

df <- mtcars

models <- list()
models[[1]] <- lm(cyl~disp, data = df)
models[[2]] <- feols(cyl~disp|as.factor(gear), data = df)
library(modelsummary)

# this works
modelsummary::modelsummary(models)

# but these do not
modelsummary::modelsummary(models, vcov = c("HC3", "cluster"))
modelsummary::modelsummary(models, vcov = c(HC3, cluster))
modelsummary::modelsummary(models, vcov = list(HC3, cluster))
modelsummary::modelsummary(models, vcov = list(vcovHC, cluster))
modelsummary::modelsummary(models, vcov = list(vcovHC, vcovHC))
modelsummary::modelsummary(models, vcov = c(vcovHC, vcovHC))

Okay--I figured out a hack, but still leaving question open in case someone finds out a more slick solution.

df <- mtcars
models <- list()
fit <- lm(cyl~disp, data = df)
models[[1]] <- coeftest(fit, vcovHC(fit, type = "HC3"))
models[[2]] <- summary(feols(cyl~disp|as.factor(gear), data = df), "cluster")
library(modelsummary)

# this works
modelsummary::modelsummary(models)
fixest modelsummary r
1 Answers

By default fixest automatically computes cluster-robust standard errors when you include fixed effects. If you set vcov=NULL, modelsummary will return the default standard errors, which will then be cluster-robust.

Alternatively, you can set vcov=~gear to compute the standard errors using the sandwich::vcovCL function under the hood.

Using version 0.10.0 of modelsummary and 0.10.4 of fixest, we can do:

library(modelsummary)
library(fixest)

models <- list(
    lm(cyl~disp, data = mtcars),
    feols(cyl ~ disp | gear, data = mtcars),
    feols(cyl ~ disp | gear, data = mtcars))

modelsummary(models,
             fmt = 6,
             vcov = list("HC3", NULL, ~gear))
Model 1 Model 2 Model 3
(Intercept) 3.188568
(0.289130)
disp 0.012998 0.012061 0.012061
(0.001310) (0.002803) (0.002803)
Num.Obs. 32 32 32
R2 0.814 0.819 0.819
R2 Adj. 0.807 0.800 0.800
R2 Within 0.614 0.614
R2 Pseudo
AIC 79.1 80.2 80.2
BIC 83.5 86.1 86.1
Log.Lik. -36.573 -36.107 -36.107
F 98.503
Std.Errors HC3 by: gear by: gear
FE: gear X X

Notice that the results are the same in the 2nd and 3rd column, and also equal to the plain fixest summary:

summary(models[[2]])
#> OLS estimation, Dep. Var.: cyl
#> Observations: 32 
#> Fixed-effects: gear: 3
#> Standard-errors: Clustered (gear) 
#>      Estimate Std. Error t value Pr(>|t|)    
#> disp 0.012061   0.002803 4.30299 0.049993 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> RMSE: 0.747808     Adj. R2: 0.799623
#>                  Within R2: 0.614334

In some past versions of modelsummary there were issues with labelling of the standard errors at the bottom of the table. I will soon be working on a more robust system to make sure the label matches the standard errors. In most cases, modelsummary calculates the robust standard errors itself, so it has full control over labelling and computation. Packages like fixest and estimatr make things a bit more tricky because they sometimes hold several SEs, and because the default is not always “classical”.

user_3423310
Related