Stata's reghdfe in R

If you are like me, you love Stata’s reghdfe command for linear regression. It offers a wide range of functionality desired in (financial) economics research, like multi dimensional fixed effects, instrumental variables and standard error clustering. Yet, for some modern data pipelines (in particular on distributed systems), it is not trivial to integrate Stata. R on the other hand has a lot of APIs that are useful in such a context (for example sparklyr package to write spark applications). In this post I am showing how you can use R’s linear fixed effects lfe packge and its felm command to replicate regression results from Stata’s reghdfe. If you are reading this post, chances are you are familiar with Stata’s reghdfe, so I won’t spend much time explaining what the below Stata code does. In a nutshell we have seven regression models here, that we are trying to replicate in R next:

cls
webuse nlswork, clear
xtset idcode year
gen cons = 1

* model 1
reg ln_w grade age ttl_exp tenure not_smsa south

* model 2
reghdfe ln_w grade age ttl_exp tenure not_smsa south, abs(idcode) 

* model 3
reghdfe ln_w grade age ttl_exp tenure not_smsa south, abs(idcode)  cluster(idcode)

* model 4
reghdfe ln_w grade age ttl_exp tenure not_smsa south, abs(year)  cluster(idcode)
xtreg ln_w grade age ttl_exp tenure not_smsa south i.year, re cluster(idcode)

* model 5
reghdfe ln_w  age ttl_exp tenure not_smsa south, abs(idcode year)  cluster(idcode)

* model 6
reghdfe ln_w  age ttl_exp tenure not_smsa south, abs(idcode year)  cluster(idcode year)

* model 7
reghdfe ln_w grade age ttl_exp tenure not_smsa south, abs(idcode year)  cluster(idcode wks_work)

I only print the results of the last model 7:

HDFE Linear regression                            Number of obs   =     26,834
Absorbing 2 HDFE groups                           F(   5,    104) =     107.26
Statistics robust to heteroskedasticity           Prob > F        =     0.0000
												  R-squared       =     0.6797
												  Adj R-squared   =     0.6215
Number of clusters (idcode)  =      4,107         Within R-sq.    =     0.0539
Number of clusters (wks_work) =        105        Root MSE        =     0.2922

  (Std. Err. adjusted for 105 clusters in idcode wks_work)
------------------------------------------------------------------------------
		 |               Robust
	 ln_wage |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
	   grade |          0  (omitted)
	     age |   .0122875    .011843     1.04   0.302    -.0111976    .0357725
	 ttl_exp |   .0329553   .0026807    12.29   0.000     .0276395    .0382712
	  tenure |   .0101173    .002155     4.69   0.000     .0058438    .0143907
	not_smsa |  -.0953025   .0139339    -6.84   0.000     -.122934   -.0676711
	   south |  -.0652122   .0168939    -3.86   0.000    -.0987133    -.031711
	   _cons |   1.137644   .3446994     3.30   0.001     .4540922    1.821196
------------------------------------------------------------------------------

Absorbed degrees of freedom:
-----------------------------------------------------+
 Absorbed FE | Categories  - Redundant  = Num. Coefs |
-------------+---------------------------------------|
      idcode |      4107        4107           0    *|
	year |        15           0          15     |
-----------------------------------------------------+
* = FE nested within cluster; treated as redundant for DoF computation

In this model we have two dimensional fixed effects (idcode and year) and two-way standard error clustering (idcode and wks_work).

First, prepare the workspace in R:

# clear work space
rm(list=ls()) # delete global environment
library(foreign)   
library(readstata13)
library(lfe)

# get sample data
df = read.dta13("http://www.stata-press.com/data/r14/nlswork.dta")

We can recreate the regression results from reghdfe using:

model7 = felm(
	  # define regression model
	  ln_wage ~  age + ttl_exp + tenure+  not_smsa  + south
	  # define fixed effects | instruments | standard errors
	  | idcode + year | 0 | idcode + wks_work, 
	  # the data
	  data = df,
	  # the method to compute SE
	  cmethod = 'cgm2',
	  # whether degree of freedom should be computed exactly
	  exactDOF=TRUE
	  )

summary(model7)

Note that you need the latest version of lfe in order to be able to use the cgm2 method. Otherwise standard errors will be different from Stata.

Call:
   felm(formula = ln_wage ~ age + ttl_exp + tenure + not_smsa +      south | idcode + year | 0 | idcode + wks_work, data = df,      exactDOF = TRUE, cmethod = "cgm2") 

Residuals:
	 Min       1Q   Median       3Q      Max 
-1.91479 -0.11707  0.00583  0.12838  2.96827 

Coefficients:
	  Estimate     Cluster s.e.   t value       Pr(>|t|)    
age       0.012287     0.011840       1.038     0.301767    
ttl_exp   0.032955     0.002680      12.297      < 2e-16 ***
tenure    0.010117     0.002154       4.696     8.13e-06 ***
not_smsa -0.095303     0.013930      -6.841     5.57e-10 ***
south    -0.065212     0.016889      -3.861     0.000196 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2921 on 22708 degrees of freedom
  (1133 observations deleted due to missingness)
Multiple R-squared(full model): 0.689   Adjusted R-squared: 0.6248 
Multiple R-squared(proj model): 0.05386   Adjusted R-squared: -0.1416 
F-statistic(full model, *iid*):10.72 on 4692 and 22708 DF, p-value: < 2.2e-16 
F-statistic(proj model): 107.3 on 5 and 104 DF, p-value: < 2.2e-16 

In the repo I stored the codes for model 1-6. The results are practically dentical (see the github issue page on remaining differences). At least they are close enough to win the star wars ;)