Feature request: Implementing `fitstat` for `fixest_multi`

[igelstorm]

It seems to me that it would be useful to be able to call fitstat on fixest_multi objects, and get a data frame with the statistics for each model. This would be analogous to/consistent with how coef, confint, coeftable, etc. currently behave. I imagine it looking something like this:

library(fixest)

m <- feols(Sepal.Length ~ Sepal.Width, iris, split = ~Species)
fitstat(m, "f", simplify = TRUE)
#>   id sample.var     sample     stat            p df1 df2
#> 1  1    Species     setosa 58.99373 6.709843e-10   1  48
#> 2  2    Species versicolor 18.35169 8.771860e-05   1  48
#> 3  3    Species  virginica 12.68707 8.434625e-04   1  48

We could also produce this output manually like this:

model_stats <- m |>
  lapply(fitstat, "f", simplify = TRUE) |>
  lapply(as.data.frame) |>
  do.call(rbind, args = _)
rownames(model_stats) <- NULL
cbind(models(m), model_stats)

But this relies on some assumptions about the internals of fixest_multi (specifically that lapply(m) lists the models in the same order as models(m)), which made me feel a bit uneasy. And it would of course be nice to have this built in 😁