Two-way Fixed Effects and Event Studies

class: center, middle, inverse, title-slide

# Two-way Fixed Effects and Event Studies
## <html>
<div style="float:left">

</div>
<hr color='#EB811B' size=1px width=0px>
</html>
### Ian McCarthy | Emory University
### Workshop on Causal Inference with Panel Data

---

<style type="text/css">
.remark-slide-content {
    font-size: 30px;
    padding: 1em 2em 1em 2em;    
}
.remark-code {
  font-size: 15px;
}
.remark-inline-code { 
    font-size: 20px;
}
</style>

# Table of contents

1. [Two-way Fixed Effects](#twfe)
2. [Event Studies](#event)
3. [In Practice](#handson)
4. [What are we estimating?](#what)

---
class: inverse, center, middle
name: twfe

# The Idea of TWFE

---
# What is TWFE?

- Just a shorthand for a common regression specification
- Fixed effects for each unit and each time period, `$\lambda_{i}$` and `$\lambda_{t}$`
- More general than 2x2 DD but same result

---
count: false

# What is TWFE?
Want to estimate `$\delta$`:

`$$y_{it} = \alpha + \delta D_{it} + \gamma_{i} + \gamma_{t} + \varepsilon,$$`<br>

where `$\gamma_{i}$` and `$\gamma_{t}$` denote a set of unit `$i$` and time period `$t$` dummy variables (or fixed effects).

---
# TWFE in Practice

```r
library(tidyverse)
library(modelsummary)
mcaid.data <- read_tsv("https://raw.githubusercontent.com/imccart/202108-cdc-methodsworkshop/master/data/medicaid-expansion/mcaid-expand-data.txt")
reg.dat <- mcaid.data %>% filter(expand_year==2014 | is.na(expand_year), !is.na(expand_ever)) %>%
  mutate(perc_unins=uninsured/adult_pop,
         post = (year>=2014), 
         treat=post*expand_ever)
```

---
count: false

# TWFE in Practice

```r
summary(lm(perc_unins ~ post + expand_ever + post*expand_ever, data=reg.dat))
```

```
## 
## Call:
## lm(formula = perc_unins ~ post + expand_ever + post * expand_ever, 
##     data = reg.dat)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.115667 -0.027106 -0.006804  0.027765  0.117597 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               0.213965   0.007180  29.799  < 2e-16 ***
## postTRUE                 -0.054068   0.008496  -6.364 7.22e-10 ***
## expand_everTRUE          -0.046326   0.009166  -5.054 7.48e-07 ***
## postTRUE:expand_everTRUE -0.018403   0.010845  -1.697   0.0908 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04187 on 304 degrees of freedom
## Multiple R-squared:  0.4995,	Adjusted R-squared:  0.4946 
## F-statistic: 101.1 on 3 and 304 DF,  p-value: < 2.2e-16
```

---
count: false

# TWFE in Practice

```r
library(fixest)
twfe <- feols(perc_unins ~ expand | State + year, data=reg.dat)
twfe$coeftable
```

```
##               Estimate  Std. Error  t value     Pr(>|t|)
## expandTRUE -0.01840269 0.003702314 -4.97059 1.220461e-06
```

---
class: inverse, center, middle
name: event

# Event Studies

---
# What is an event study?

Event study is poorly named:
- In finance, even study is just an *interrupted time series*
- In econ and other areas, we usually have a treatment/control group *and* a break in time

---
count: false

# What is an event study?

- Allows for different effect estimates at each time period (maybe effects phase in over time or dissipate)
- Visually very appealing
- Offers easy "smell test" for parallel trends assumption

---
count: false

# What is an event study?

Estimate something akin to...
`$$y_{it} = \gamma_{i} + \gamma_{t} + \sum_{\tau = -q}^{-1}\delta_{\tau} D_{i \tau} + \sum_{\tau=0}^{m} \delta_{\tau}D_{i \tau} + x_{it} + \epsilon_{it},$$`

where `$q$` captures the number of periods before the treatment occurs and `$m$` captures periods after treatment occurs.

---
# How to do an event study?

First create all of the treatment/year interactions:

```r
event.dat <- reg.dat %>%
  mutate(expand_2012 = expand_ever*(year==2012),
         expand_2013 = expand_ever*(year==2013),
         expand_2014 = expand_ever*(year==2014),
         expand_2015 = expand_ever*(year==2015),
         expand_2016 = expand_ever*(year==2016),
         expand_2017 = expand_ever*(year==2017),
         expand_2018 = expand_ever*(year==2018))
```

---
count: false

# How to do an event study?

Second, run regression with full set of interactions and group/year dummies:

```r
event.ins.reg <- lm(perc_unins ~ expand_2012 + expand_2014 + 
                      expand_2015 + expand_2016 + expand_2017 + 
                      expand_2018 + factor(year) + factor(State), data=event.dat)
point.est <- as_tibble(c(event.ins.reg$coefficients[c("expand_2012","expand_2014","expand_2015",
                                            "expand_2016","expand_2017","expand_2018")]),
                       rownames = "term")
ci.est <- as_tibble(confint(event.ins.reg)[c("expand_2012","expand_2014","expand_2015",
                                   "expand_2016","expand_2017","expand_2018"),],
                    rownames = "term")
```

---
count: false

# How to do an event study?

Third, organize results into a new dataset:

```r
point.est <- point.est %>% rename(estimate = value)
ci.est <- ci.est %>% rename(conf.low = `2.5 %`, conf.high = `97.5 %`)
new.row <- tibble(
  term = "expand_2013",
  estimate = 0,
  conf.low = 0,
  conf.high = 0,
  year = 2013
)

event.plot.dat <- point.est %>%
  left_join(ci.est, by=c("term")) %>%
  mutate(year = c(2012, 2014, 2015, 2016, 2017, 2018)) %>%
  bind_rows(new.row) %>%
  arrange(year)
```

---
count: false

# How to do an event study?

Finally, plot coefficients and confidence intervals

.plot-callout[
<img src="twfe-cdc202108_files/figure-html/es-plot-small-1.png" style="display: block; margin: auto;" />
]

---
count: false

# How to do an event study?

---
class: inverse, center, middle
name: handson

# Seeing things in action

---
# Things to address

1. "Event time" vs calendar time
2. Define baseline period
3. Choose number of pre-treatment and post-treatment coefficients

---
# Event time vs calendar time

Essentially two "flavors" of event studies

1. Common treatment timing
2. Differential treatment timing

---
# Define baseline period

- Must choose an "excluded" time period (as in all cases of group dummy variables)
- Common choice is `$t=-1$` (period just before treatment)
- Easy to understand with calendar time
- For event time...manually set time to `$t=-1$` for all untreated units

---
# Number of pre-treatment and post-treatment periods

- On event time, sometimes very few observations for large lead or lag values
- Medicaid expansion example: Late adopting states have fewer post-treatment periods
- Norm is to group final lead/lag periods together

---
# Common treatment timing

.pull-left[
**Stata**<br>

```stata
ssc install reghdfe

insheet using "https://raw.githubusercontent.com/imccart/202108-cdc-methodsworkshop/master/data/medicaid-expansion/mcaid-expand-data.txt", clear
gen perc_unins=uninsured/adult_pop
keep if expand_year=="2014" | expand_year=="NA"
drop if expand_ever=="NA"
gen post=(year>=2014)
gen treat=(expand_ever=="TRUE")
gen treat_post=(expand=="TRUE")

reghdfe perc_unins treat##ib2013.year, absorb(state)
gen coef = .
gen se = .
forvalues i = 2012(1)2018 {
    replace coef = _b[1.treat#`i'.year] if year == `i'
    replace se = _se[1.treat#`i'.year] if year == `i'
}

* Make confidence intervals
gen ci_top = coef+1.96*se
gen ci_bottom = coef - 1.96*se

* Limit ourselves to one observation per year
keep year coef se ci_*
duplicates drop

* Create connected scatterplot of coefficients
* with CIs included with rcap 
* and a line at 0 from function
twoway (sc coef year, connect(line)) (rcap ci_top ci_bottom year) ///
    (function y = 0, range(2012 2018)), xtitle("Year") ///
    caption("Estimates and 95% CI from Event Study")
```
]

.pull-right[
**R**<br>

```r
library(tidyverse)
library(modelsummary)
library(data.table)
library(fixest)
mcaid.data <- read_tsv("https://raw.githubusercontent.com/imccart/202108-cdc-methodsworkshop/master/data/medicaid-expansion/mcaid-expand-data.txt")
reg.dat <- as.data.table(mcaid.data) %>% 
  filter(expand_year==2014 | is.na(expand_year), !is.na(expand_ever)) %>%
  mutate(perc_unins=uninsured/adult_pop,
         post = (year>=2014), 
         treat=post*expand_ever)

mod.twfe <- feols(perc_unins~i(year, expand_ever, ref=2013) | State + year,
                  cluster=~State,
                  data=reg.dat)
iplot(mod.twfe, 
      xlab = 'Time to treatment',
      main = 'Event study')
```
]

---
# Comparing results

.plot-callout[
<img src="twfe-cdc202108_files/figure-html/es-plot1-1.png" style="display: block; margin: auto;" />
]

---
# Differential treatment timing

- Now let's work with the full Medicaid expansion data
- Includes late adopters
- Requires putting observations on "event time"

---
count: false

# Differential treatment timing

.pull-left[
**Stata**<br>

```stata
ssc install reghdfe

insheet using "https://raw.githubusercontent.com/imccart/202108-cdc-methodsworkshop/master/data/medicaid-expansion/mcaid-expand-data.txt", clear
gen perc_unins=uninsured/adult_pop
drop if expand_ever=="NA"
replace expand_year="." if expand_year=="NA"
destring expand_year, replace
gen event_time=year-expand_year
replace event_time=-1 if event_time==.

forvalues l = 0/4 {
	gen L`l'event = (event_time==`l')
}
forvalues l = 1/2 {
	gen F`l'event = (event_time==-`l')
}
gen F3event=(event_time<=-3)

reghdfe perc_unins F3event F2event L0event L1event L2event L3event L4event, absorb(state year) cluster(state)
gen coef = .
gen se = .
forvalues i = 2(1)3 {
    replace coef = _b[F`i'event] if F`i'event==1
    replace se = _se[F`i'event] if F`i'event==1
}
forvalues i = 0(1)4 {
    replace coef = _b[L`i'event] if L`i'event==1
    replace se = _se[L`i'event] if L`i'event==1
}
replace coef = 0 if F1event==1
replace se=0 if F1event==1

* Make confidence intervals
gen ci_top = coef+1.96*se
gen ci_bottom = coef - 1.96*se

* Limit ourselves to one observation per year
keep if event_time>=-3 & event_time<=4
keep event_time coef se ci_*
duplicates drop

* Create connected scatterplot of coefficients
* with CIs included with rcap 
* and a line at 0 from function
sort event_time
twoway (sc coef event_time, connect(line)) (rcap ci_top ci_bottom event_time) ///
    (function y = 0, range(-3 4)), xtitle("Time") ///
    caption("Estimates and 95% CI from Event Study") xlabel(-3(1)4)
```
]

.pull-right[
**R**<br>

```r
library(tidyverse)
library(modelsummary)
library(data.table)
library(fixest)
mcaid.data <- read_tsv("https://raw.githubusercontent.com/imccart/202108-cdc-methodsworkshop/master/data/medicaid-expansion/mcaid-expand-data.txt")
reg.dat <- as.data.table(mcaid.data) %>% 
  filter(!is.na(expand_ever)) %>%
  mutate(perc_unins=uninsured/adult_pop,
         post = (year>=2014), 
         treat=post*expand_ever,
         time_to_treat = ifelse(expand_ever==FALSE, 0, year-expand_year),
         time_to_treat = ifelse(time_to_treat < -3, -3, time_to_treat))

mod.twfe <- feols(perc_unins~i(time_to_treat, expand_ever, ref=-1) | State + year,
                  cluster=~State,
                  data=reg.dat)
iplot(mod.twfe, 
      xlab = 'Time to treatment',
      main = 'Event study')
```
]

---
class: inverse, center, middle
name: what

# What are we estimating?

---
# Problems with TWFE

- Recall goal of estimating ATE or ATT
- TWFE and 2x2 DD identical with homogeneoues effects and common treatment timing
- Otherwise...TWFE is biased and inconsistent for ATT

---
# Inutition

- OLS is a weighted average of all 2x2 DD groups
- Weights are function of size of subsamples, size of treatment/control units, and timing of treatment
- Units treated in middle of sample receive larger weights
- Prior-treated units act as controls for late-treated units

--
Just the length of the panel will change the estimate!

---
# Does it really matter?

- Definitely! But how much?
- Large treatment effects for early treated units could reverse the sign of final estimate
- Let's explore this nice Shiny app from Kyle Butts: [Bacon-Decomposition Shiny App](https://kyle-butts.shinyapps.io/did_bacon/).