title: R + FE Regression with many fixed effects where I extract the fe coefficients and standard errors
body: <p>I have a dataset with around 2 million individiuals. A simplified version of the dataset contains three types of variables: (1) dependent variable (y), (2) covariates (x1, x2), (3) regional fixed effects (reg, around 2000 thousand).</p>
  <p>(For a second regression at the regional level) I need to extract the regional fixed effects and their standard errors. The main problem is that my computer gets stuck, apparently because of a lack of memory (I have 16GB, 4-core macbook pro) and at the moment I don't have access to a more powerful machine.</p>
  <p>Strategies used:</p>
  <ol>
  <li><p>The simple <code>lm(y ~ x1 + x2 + reg,data=DT)</code> gets stuck</p>
  <ul>
  <li>I've tried reducing the number of regional dummies (e.g. dropping individuals in different regions) and this works when I have around 500 dummies. Nonetheless this implies dropping half of my dataset.</li>
  </ul>
  </li>
  <li><p>Using lm_robust:</p>
  </li>
  </ol>
  <ul>
  <li><code>lm_robust(y ~ x1 + x2, fixed_effects= ~ reg, se_type='stata',data = DT)</code> works, is relatively fast, gives the coefficients of the fixed effects (fixed_effects), but not the standard errors</li>
  <li><code>lm_robust(y ~ x1 + x2 + reg, se_type='stata',data = DT)</code> gives me the fixed effects coefficients and standard errors, BUT gets stuck when I use more than 500 regions</li>
  </ul>
  <ol start="3">
  <li><p>Using either felm or feols</p>
  <ul>
  <li><code>felm(y ~ x1 + x2 | reg, data=DT)</code></li>
  <li><code>feols(y ~ x1 + x2 | reg, data=DT)</code></li>
  </ul>
  <p>Both methods work, are very fast, BUT DON'T HAVE IN THE OUTPUT THE FE COEFFICIENTS AND STANDARD ERRORS</p>
  </li>
  </ol>
  <p>One strategy I've considered is to use felm/feols and use the outputs &quot;residuals&quot; (residuals of the full system, with dummies) and &quot;r.residuals&quot; (residuals resulting from predicting without the dummies) to back out the regional residuals (i.e. for region i,  hat(reg)_i = (y-hat(b1)*x1 - hat(b2)*x2 - hat(reg_i)) - (y-hat(b1)*x1 - hat(b2)*x2 ). The problem with this is that I wouldn't get the standard error. One solution for this would be to use bootstrapping to compute the standard errors. This should work (I haven't tried it yet), but it seems to be too convoluted.</p>