Two Ways to fit Two-Way Fixed Effects in R

1. Introduction

Recently, a friend asked me how to fit a two-way fixed effects model in R. A fixed effects model is a regression model in which the intercept of the model is allowed to move across individuals and groups. We most often see it in panel data contexts. Two-way fixed effects have seen massive interest from the methodological community. Some recent papers of interest are Imai and Kim 2019, Goodman-Bacon 2019, and Abraham and Sun 2018.

In this post, I show two ways to fit two-way fixed effects models. The first is the least squares with dummy variables approach, and the second is a fixed effects approach.

2. Packages

To show that both methods work, use a reproducible example from a dataset located in the package bacondecomp.

# for a dataset
library(bacondecomp) 
# for robust standard error estimation
library(lmtest) 
# To calculate correct vcov matrix with 2WFE
library(multiwayvcov) 
# For a package way to do FE
library(plm)

3. Estimation

Least Squares Dummy Variables

Here is the way to fit with lm(). Because we only have one independent variable of interest, I subset coeftest() to just that variable. The rest of the variables will be all the factor dummies.

df <- bacondecomp::castle
# The way with lm. 
fit_tw <- lm(l_homicide ~ post + factor(state) + factor(year), 
             data = df)

# The coefficient of interest is the second in this object. 
vcov_tw <- multiwayvcov::cluster.vcov(fit_tw, 
                                      cbind(df$state, df$year),
                                      use_white = F, 
                                      df_correction = F)
# Just get coefficient of interest
# Here it's the second row from coeftest
coeftest(fit_tw, vcov_tw)[2,] 

##   Estimate Std. Error    t value   Pr(>|t|) 
## 0.08181162 0.05671410 1.44252696 0.14979397

Alternatively, we can do it with plm(). Under this method, we are calculating a fixed effects estimator. We will get the same result, but this way is more computationally efficient, especially as our models become more complex.

fit_plm <- plm(l_homicide ~ post, 
               data = df, 
               index = c("state", "year"), 
               model = "within", 
               effect = "twoways")

# Note how this is functionally identical to the lm() way 
coeftest(fit_plm, vcov = vcovHC, type = "HC1")

## 
## t test of coefficients:
## 
##      Estimate Std. Error t value Pr(>|t|)
## post 0.081812   0.057748  1.4167   0.1572

Note that both methods return functionally identical answers. For other ways, check out LOST Stats