predict() SEs for fixest models with FEs - when to re-estimate?

[arcruz0]

Hi and thanks for the great package.

I realize that there is no natural way to compute SEs for predictions in fixest models with FEs (for reference, see these fixest and pyfixest issues).

One can get prediction SEs using marginaleffects, which I reckon uses the Delta method. Occasionally, I get quite different results with this approach vs estimating the model with factor variables, especially downstream in e.g. plot_predictions() (see the example below).

I would like to understand this process better, and in particular when would it make sense to re-estimate the model with factor variables instead of using the fixest object. I'd very much appreciate any pointers, thoughts, etc.

library(fixest)
library(marginaleffects)
library(patchwork)

m1 <- lm(Ozone ~ Wind + factor(Month), airquality)
m2 <- feols(Ozone ~ Wind | Month, airquality, vcov = "iid")
#> NOTE: 37 observations removed because of NA values (LHS: 37).

m1 |> predictions() |> head()
#> 
#>  Estimate Std. Error    z Pr(>|z|)    S 2.5 % 97.5 %  Df
#>     42.46       5.66 7.50   <0.001 43.8 31.36   53.6 Inf
#>     39.67       5.46 7.26   <0.001 41.3 28.97   50.4 Inf
#>     18.31       4.96 3.69   <0.001 12.1  8.59   28.0 Inf
#>     23.42       4.89 4.79   <0.001 19.2 13.83   33.0 Inf
#>      7.63       5.46 1.40    0.162  2.6 -3.06   18.3 Inf
#>     36.88       5.29 6.98   <0.001 38.3 26.52   47.2 Inf
#> 
#> Type: response
m2 |> predictions() |> head()
#> 
#>  Estimate Std. Error    z Pr(>|z|)    S   2.5 % 97.5 %  Df
#>     42.46       5.19 8.17  < 0.001 51.6  32.276   52.6 Inf
#>     39.67       5.61 7.07  < 0.001 39.2  28.665   50.7 Inf
#>     18.31       8.84 2.07  0.03839  4.7   0.979   35.6 Inf
#>     23.42       8.07 2.90  0.00371  8.1   7.600   39.2 Inf
#>      7.63      10.46 0.73  0.46548  1.1 -12.864   28.1 Inf
#>     36.88       6.04 6.11  < 0.001 29.9  25.054   48.7 Inf
#> 
#> Type: response

plot_predictions(m1, "Wind") + ggplot2::ggtitle("stats::lm()") + 
  plot_predictions(m2, "Wind") + ggplot2::ggtitle("fixest::feols()")

^{Created on 2025-06-09 with reprex v2.1.1}

(The example is adapted from #531).

[vincentarelbundock]

This is a very question, but I'm not sure I have a great answer. I opened an issue here with some details for further discussion: #1487