Wild cluster bootstrap after linear model with multiple fixed effects I am running a linear model with multiple fixed effects. I suspect that there is spatial correlation in my data so I cluster the SEs (46 clusters). But I am worried that the standard errors generated by felm::lfe (and reghdfe in Stata) might be too large. So I want to alternatively also compute wild cluster bootstrap (WCB) based confidence intervals for the main parameters. 
This does not seem to be an easy thing to do for panel models with multiple fixed effects. Here's a summary of what I found:


*

*Stata: Current procedures to do WCB on stata (boottest, clustse,
cgmwildboot) do not work with procedures that implement multiple
fixed effects (reghdfe, areg). reghdfe itself does not currently
support WCB.

*R: The felm package that seems to be the best option
for linear panel models with multiple FE and clustered SE does not support WCB. In fact,
it does not even support predict after the regression.


So, there seem to be two possible solutions:


*

*Workaround: WCB procedures on stata work with one level of FE (for example, boottest). So, converting the reghdfe regression to include dummies and absorbing the one FE with largest set would probably work with boottest. But this is likely to be crazily expensive on memory and processing capacity. If I have access to high-performance servers I could run this code on those servers.

*Code own WCB function in R: The second option is to get predictions from felm regression (using this process from a previous answer on Stack Overflow). Then one could write a function to do the wild cluster bootstrap and get the confidence interval for estimated coefficients.


There are probably better solutions out there. What would the community suggest as best practice to get wild cluster bootstrap based CI after a panel regression with multiple fixed effects?
 A: The fwildclusterboot R package implements a wild cluster bootstrap and allows for multiple fixed effects. The package is a port of the boottest package in Stata.
fwildclusterboot accepts objects of type felm and fixest where multiple fixed effects are out-projected in the estimation stage. In the bootstrap inference stage, it converts all fixed effects to sets of dummy variables but allows the user to specify one fixed effect to be projected out in the bootstrap - just as in the boottest.ado Stata routine.
library(fwildclusterboot)
library(fixest)

# load data set voters included in fwildclusterboot
data(voters)

# estimate the regression model via feols
feols_fit <- feols(proposition_vote ~ treatment + ideology1 + log_income | Q1_immigration , data = voters)
# boottest on an object of type fixest, no fixed effect used for inference
boot_feols <- boottest(feols_fit, clustid = "group_id1", param = "treatment", B = 99999, tol = 1e-08)
# now, project out "Q1_immigration" during bootstrap
boot_feols_fe <- boottest(feols_fit, clustid = "group_id1", param = "treatment", B = 99999, fe = "Q1_immigration", tol = 1e-08)

library(modelsummary)
msummary(list(boot_feols, boot_feols_fe), 
         estimate = "{estimate} ({p.value})", 
         statistic = "[{conf.low}, {conf.high}]", 
         digits = 3, 
         output = "markdown")


#|          |    Model 1     |    Model 2     |
#|:---------|:--------------:|:--------------:|
#|treatment | 0.073 (0.001)  | 0.073 (0.001)  |
#|          | [0.033, 0.113] | [0.033, 0.113] |
#|Num.Obs.  |      300       |      300       |
#|R2        |     0.316      |     0.316      |
#|R2 Adj.   |     0.288      |     0.288      |
#|R2 Within |     0.052      |     0.052      |
#|R2 Pseudo |                |                |
#|AIC       |     -84.1      |     -84.1      |
#|BIC       |     -35.9      |     -35.9      |
#|Log.Lik.  |     55.025     |     55.025     |
```

