R - lfe (felm) fit model with fixed effects only

Question

I'm using the felm() function from the lfe package to fit linear models with large numbers of fixed effects. I would like to be able to fit a model using only fixed effects. For instance, I'd like to be able to know the R^2 of such a model, and potentially compare it to that of a model with a larger set of predictors. Consider the example below:

library(lfe)

N = 1000

A = sample(1:3, N, replace = TRUE)
B = sample(1:5, N, replace = TRUE)

C = A + B + rnorm(N)

Data = data.frame(A, B, C)
Data$A = as.factor(A)
Data$B = as.factor(B)

summary(felm(C ~ 1 | A + B, data = Data))

Which just gives in response:

Call:
   felm(formula = C ~ 1 | A + B, data = Data) 

Residuals:
    Min      1Q  Median      3Q     Max 
-3.8101 -0.6750  0.0014  0.6765  4.4254 

Coefficients:
(No coefficients)

similarly, if I use: names(summary(felm(C ~ 1 | A + B, data = Data))) I get:

[1] "residuals" "p"         "Pp"        "call"   

whereas for models where I specify variables other than the FEs, I get many more attributes in the summary, including the R^2.

I also tried just adding a variable to my data that is a set of ones, but this did not work.

I can get this easily with the regular lm() function (summary(lm(C ~ A + B, data = Data))), but that then deprives me the value of the felm() function:

Answer 1

the fixest package (which handles high-dimensional fixed effects at least as smoothly as lfe) returns a set of statistics, e.g. the adjusted R-squared:

> summary(fixest::feols(C ~ 1 | A + B, data = Data))
# OLS estimation, Dep. Var.: C
# Observations: 1,000 
# Fixed-effects: A: 3,  B: 5
# RMSE: 0.994166   Adj. R2: 0.733952