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

Question

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)

Answer 1

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