# 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