title: R - lfe (felm) fit model with fixed effects only
body: <p>I'm using the <code>felm()</code> function from the <code>lfe</code> 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:</p>
  
  <pre><code>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))
  </code></pre>
  
  <p>Which just gives in response:</p>
  
  <pre><code>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)
  </code></pre>
  
  <p>similarly, if I use: <code>names(summary(felm(C ~ 1 | A + B, data = Data)))</code> I get:</p>
  
  <pre><code>[1] "residuals" "p"         "Pp"        "call"   
  </code></pre>
  
  <p>whereas for models where I specify variables other than the FEs, I get many more attributes in the summary, including the R^2.</p>
  
  <p>I also tried just adding a variable to my data that is a set of ones, but this did not work.</p>
  
  <p>I can get this easily with the regular <code>lm()</code> function (<code>summary(lm(C ~ A + B, data = Data))</code>), but that then deprives me the value of the <code>felm()</code> function:</p>