---
title: "Extended TWFE and other DiD estimators"
date: "2024-12-09"
---
## Staggered DiD without treatment reversal
This blog is adopted from the tutorial of paneltools package.
https://yiqingxu.org/tutorials/panel.html#Installing_Packages
I added Wooldridge's ETWFE, which is my favorate estimator with staggered DiD. It is easy to understand and implement. It also works with nonlinear model such as logit and poisson.
```{r}
#| message: false
#| warning: false
#| cache: true
#devtools:: install_github("xuyiqing/paneltools")
library(paneltools)
data(paneltools)
library(tidyverse)
hh2019 <- as.tibble(hh2019)
hh2019
```
This data is from Hainmueller and Hopkins (2019), in which the authors investigate the effects of indirect democracy (versus direct democracy) on naturalization rates in Switzerland using municipality-year level panel data from 1991 to 2009. The study shows that a switch from direct to indirect democracy increased naturalization rates by an average of 1.22 percentage points (Model 1 of Table 1 in the paper).
```{r}
#| message: false
#| warning: false
#| cache: true
#devtools:: install_github("xuyiqing/panelView")
library(panelView)
panelview(nat_rate_ord ~ indirect, data = hh2019, index = c("bfs","year"),
xlab = "Year", ylab = "Unit", display.all = T,
gridOff = TRUE, by.timing = TRUE)
```
### TWFE
First we do a TWFE model on it.
```{r}
#| message: false
#| warning: false
#| cache: true
library(fixest)
model_twfe <- feols(nat_rate_ord~indirect|bfs+year,
data=hh2019, cluster = "bfs")
summary( model_twfe)
```
TWFE has issues when there is staggered adoption and heterogeneous treatment effects.
Goodman-Bacon (2021) demonstrates that the two-way fixed-effects (TWFE) estimator in a staggered adoption setting can be represented as a weighted average of all possible 2x2 difference-in-differences (DID) estimates between different cohorts. However, when treatment effects change over time heterogeneously across cohorts, the “forbidden” comparisons that use post-treatment data from early adopters as controls for late adopters may introduce bias in the TWFE estimator. We employ the procedure outlined in Goodman-Bacon (2021) and decompose the TWFE estimate using the command bacon in the bacondecomp package.
"Negative weights" can be the other issue. Basically it is a weighted average of comparison between different cohorts. The weights can be negative in a TWFE model.
```{r}
#| message: false
#| warning: false
#| cache: true
library(bacondecomp)
df_bacon <- bacon(nat_rate_ord~indirect,
data = hh2019,
id_var = "bfs",
time_var = "year")
coef_bacon <- sum(df_bacon$estimate * df_bacon$weight)
coef_bacon
ggplot(df_bacon) +
aes(x = weight, y = estimate, shape = factor(type), color = factor(type)) +
labs(x = "Weight", y = "Estimate", shape = "Type", color = 'Type') +
geom_point()
```
The problem shown in Bacon decomposion here is the comparison between "Later" and "Earlier Treated". It's a "forbiden" comparison.
TWFE with event history type of graph.
```{r}
#| message: false
#| warning: false
#| cache: true
df.use <- hh2019 |>
na.omit()
df <- as.data.frame(
hh2019 |>
group_by(bfs) |>
mutate(treatment_mean = mean(indirect,na.rm = TRUE))
)
df.use <- df[which(df$treatment_mean<1),]
df.twfe <- as.tibble(get.cohort(df.use,D="indirect",index=c("bfs","year"),start0 = TRUE))
# time_to_treatment is NA when the unit is always treated, or never treated. It needs to be set to an arbitrary value for the twfe model.
df.twfe$treat <- as.numeric(df.twfe$treatment_mean>0) # drop always treated units
df.twfe[which(is.na(df.twfe$Time_to_Treatment)),'Time_to_Treatment'] <- 0 # can be an arbitrary value
twfe.est <- feols(nat_rate_ord ~ i(Time_to_Treatment, treat, ref = -1)| bfs + year,
data = df.twfe, cluster = "bfs")
summary(twfe.est)
twfe.output <- as.matrix(twfe.est$coeftable)
twfe.output <- as.data.frame(twfe.output)
twfe.output$Time <- c(c(-18:-2),c(0:17))+1
p.twfe <- esplot(twfe.output,Period = 'Time',Estimate = 'Estimate',
SE = 'Std. Error', xlim = c(-12,10))
p.twfe
```
### Sun and Abraham
Sun and Abraham (2021) estimator is a weighted average of the average treatment effect (ATT) estimates for each cohort, obtained from a TWFE regression that includes cohort dummies fully interacted with indicators of relative time to the treatment’s onset. The iw estimator is robust to heterogeneous treatment effects (HTE) and can be implemented using the command sunab in the package fixest.
```{r}
#| message: false
#| warning: false
#| cache: true
df.sa <- df.twfe
df.sa[which(is.na(df.sa$FirstTreat)),"FirstTreat"] <- 1000 # replace NA with an arbitrary number
model.sa.1 <- feols(nat_rate_ord~sunab(FirstTreat,year)|bfs+year,
data = df.sa, cluster = "bfs")
summary(model.sa.1,agg = "ATT")
sa.output <- as.matrix(model.sa.1$coeftable)
sa.output <- as.data.frame(sa.output)
sa.output$Time <- c(c(-18:-2),c(0:17))+1
p.sa <- esplot(sa.output,Period = 'Time',Estimate = 'Estimate',
SE = 'Std. Error', xlim = c(-12,10))
p.sa
```
## Callaway and Sant’Anna
Callaway and Sant’Anna (2021) estimator is a weighted average of the average treatment effect (ATT) estimates for each cohort, obtained from a TWFE regression that includes cohort dummies fully interacted with indicators of relative time to the treatment’s onset. The iw estimator is robust to heterogeneous treatment effects (HTE) and can be implemented using the command callaway in the package fixest.
```{r}
#| message: false
#| warning: false
#| cache: true
library(did)
df.cs <- df.twfe
df.cs[which(is.na(df.cs$FirstTreat)),"FirstTreat"] <- 0 # replace NA with 0
cs.est.1 <- att_gt(yname = "nat_rate_ord",
gname = "FirstTreat",
idname = "bfs",
tname = "year",
xformla = ~1,
control_group = "nevertreated",
allow_unbalanced_panel = TRUE,
data = df.cs,
est_method = "reg")
cs.est.att.1 <- aggte(cs.est.1, type = "simple", na.rm=T, bstrap = F)
print(cs.est.att.1)
cs.att.1 <- aggte(cs.est.1, type = "dynamic",
bstrap=FALSE, cband=FALSE, na.rm=T)
print(cs.att.1)
cs.output <- cbind.data.frame(Estimate = cs.att.1$att.egt,
SE = cs.att.1$se.egt,
time = cs.att.1$egt + 1)
p.cs.1 <- esplot(cs.output,Period = 'time',Estimate = 'Estimate',
SE = 'SE', xlim = c(-12,10))
p.cs.1
```
### Imai, Kim, and Wang
PanelMatch is an R package implementing a set of methodological tools proposed by Imai, Kim, and Wang (2021) that enables researchers to apply matching methods for causal inference on time-series cross-sectional data with binary treatments. The package includes implementations of matching methods based on propensity scores and Mahalanobis distance, as well as weighting methods. PanelMatch enables users to easily calculate a variety of possible quantities of interest, along with standard errors. The software is flexible, allowing users to tune the matching, refinement, and estimation procedures with a large number of parameters. The package also offers a variety of visualization and diagnostic tools for researchers to better understand their data and assess their results.
```{r}
#| message: false
#| warning: false
#| cache: true
library(PanelMatch)
df.pm <- as.data.frame(df.twfe)
# we need to convert the unit and time indicator to integer
df.pm[,"bfs"] <- as.integer(as.factor(df.pm[,"bfs"]))
df.pm[,"year"] <- as.integer(as.factor(df.pm[,"year"]))
df.pm <- df.pm[,c("bfs","year","nat_rate_ord","indirect")]
# Create PanelData object (new API)
pd <- PanelData(df.pm, unit.id = "bfs", time.id = "year",
treatment = "indirect", outcome = "nat_rate_ord")
PM.results <- PanelMatch(lag=3,
panel.data = pd,
refinement.method = "none",
qoi = "att",
lead = c(0:3),
match.missing = TRUE)
## For pre-treatment dynamic effects
PM.results.placebo <- PanelMatch(lag=3,
panel.data = pd,
refinement.method = "none",
qoi = "att",
lead = c(0:3),
match.missing = TRUE,
placebo.test = TRUE)
PE.results.pool <- PanelEstimate(PM.results, panel.data = pd, pooled = TRUE)
summary(PE.results.pool)
PE.results <- PanelEstimate(PM.results, panel.data = pd)
PE.results.placebo <- placebo_test(PM.results.placebo, panel.data = pd, plot = F)
est_lead <- as.vector(estimates(PE.results))
est_lag <- tryCatch(as.vector(estimates(PE.results.placebo)),
error = function(e) as.vector(PE.results.placebo$estimates))
sd_lead <- tryCatch(apply(PE.results$bootstrapped.estimates,2,sd),
error = function(e) rep(NA, length(est_lead)))
sd_lag <- tryCatch(apply(PE.results.placebo$bootstrapped.estimates,2,sd),
error = function(e) rep(NA, length(est_lag)))
coef <- c(est_lag, 0, est_lead)
sd <- c(sd_lag, 0, sd_lead)
pm.output <- cbind.data.frame(ATT=coef, se=sd, t=c(-2:4))
p.pm <- esplot(data = pm.output,Period = 't',
Estimate = 'ATT',SE = 'se')
p.pm
```
### Liu, Wang and Xu
Liu, Wang, and Xu (2022) proposes a method for estimating the fixed effect counterfactual model by utilizing the untreated observations. It is essentially an imputation method that uses the untreated observations to estimate the counterfactual outcome for the treated.
```{r}
#| message: false
#| warning: false
#| cache: true
library(fect)
out.fect <- fect(nat_rate_ord~indirect, data = df,
index = c("bfs","year"),
method = 'fe', se = TRUE)
print(out.fect$est.avg)
fect.output <- as.matrix(out.fect$est.att)
print(fect.output)
fect.output <- as.data.frame(fect.output)
fect.output$Time <- c(-17:18)
p.fect <- esplot(fect.output,Period = 'Time',Estimate = 'ATT',
SE = 'S.E.',CI.lower = "CI.lower",
CI.upper = 'CI.upper',xlim = c(-12,10))
p.fect
```
Borusyak, Jaravel, and Spiess (2021) with package didimputation almost gives the same result as Liu, Wang, and Xu (2022).
### Wooldridge
Wooldridge 2021, also called extended TWFE (ETWFE), is to allow more flexibiltiy to make TWFE work properly. The problems with TWFE is that we are not specify it flexibly enough. If we allow heterogeneous treatment effect, then we cannot expect the canonical TWFE to work properly.
```{r}
#| message: false
#| warning: false
#| cache: true
df.wool <- df.twfe
df.wool[which(is.na(df.wool$FirstTreat)),"FirstTreat"] <- 0 # replace NA with 0
library(etwfe)
model.wool = etwfe(
fml = nat_rate_ord ~ 1, # outcome ~ controls
tvar = year, # time variable
gvar = FirstTreat, # group variable
data = df.wool, # dataset
vcov = ~bfs, # vcov adjustment (here: clustered)
cgroup = "never"
)
model.wool
emfx(model.wool)
mod_es = emfx(model.wool, type = "event", post_only=FALSE)
mod_es
wool.output <- mod_es |>
filter(event != -1) |>
mutate(Time = event + 1) |>
as.data.frame()
p.wool <- esplot(wool.output,Period = 'Time',Estimate = 'estimate',
SE = 'std.error', xlim = c(-12,10))
p.wool
```
Note that in order to get pre-event estimates, we need to set "cgroup=never" in the etwfe function. Otherwise, ETWFE will only give us post-event estimates. Under the default setting, the control group comprises the “not-yet” treated units, therefore all pre-treatment effects are mechanistically set to zero. By setting "cgroup=never", we only use never treated units as the control group. This way we can get pre-treatment estimates, but we are not using all available information for the control group, therefore the estimates may not be as efficient as the default setting.
Also we don't have control variable here. If we have, it is very easy to incorporate into ETWFE. Wooldridge used the Mundlak device to incorporate control variables, which is a very simple and intuitive way to do it.
If we do the default setting and plot only post-event estimates:
```{r}
#| message: false
#| warning: false
#| cache: true
model.wool2 = etwfe(
fml = nat_rate_ord ~ 1, # outcome ~ controls
tvar = year, # time variable
gvar = FirstTreat, # group variable
data = df.wool, # dataset
vcov = ~bfs # vcov adjustment (here: clustered)
)
model.wool2
emfx(model.wool2)
mod_es2 = emfx(model.wool, type = "event")
mod_es2
wool.output2 <- mod_es2 |>
filter(event != -1) |>
mutate(Time = event + 1) |>
as.data.frame()
p.wool2 <- esplot(wool.output2, Period = 'Time',Estimate = 'estimate',
SE = 'std.error', xlim = c(0,10))
p.wool2
```
## Staggered DiD with treatment reversal
With treatment reversal (ie., the treated units go back to control), Wooldridge's method still works. However, I don't think package "etwfe" has it implemented.
So far, packages "fect", "PanelMatch" can handle treatment reversal.
We can probably manually do ETWFE with treatment reversal. See Wooldridge's dropbox for the code in Stata. The basic idea is to treat cohorts as different treatment beginning time and exit time combinations. Then do the same regression with interactions of cohort dummies with the time dummies. It is very flexible. The tricky part is to get the marginal effect and its standard error right.
