Event Study

Event Study Skill

This skill covers everything needed to conduct event studies in R - from the traditional two-way fixed effects (TWFE) approach through modern robust estimators that handle staggered treatment timing and heterogeneous treatment effects.

When to Use This Skill

Use this skill when the user wants to:

• Create event study plots (coefficient plots around a treatment event)
• Test for parallel pre-trends in a difference-in-differences design
• Estimate dynamic treatment effects
• Work with staggered treatment adoption (units treated at different times)
• Implement any of the modern DiD estimators (Callaway & Sant'Anna, Sun & Abraham, etc.)
• Diagnose problems with TWFE under treatment effect heterogeneity

Quick Decision Tree

Before writing any code, determine the user's situation:

1. Is treatment timing the same for all treated units?

• Yes → Traditional TWFE event study is fine. See traditional-event-study.md
• No (staggered) → Go to step 2.

2. Is there reason to expect heterogeneous treatment effects across cohorts or over time?

• No strong reason → Traditional TWFE may still be OK, but consider robust alternatives.
• Yes or unsure → Use a robust estimator. See modern-extensions.md

3. What is the user's goal?

• Quick visualization of pre-trends → Traditional approach, fast and familiar.
• Publication-quality estimation → Recommend a robust estimator + traditional as comparison.
• Diagnostic check → Bacon decomposition to understand TWFE weights.

Core Concepts (Brief)

An event study plot displays estimated coefficients for leads (pre-treatment periods) and lags (post-treatment periods) relative to a baseline period (typically one period before treatment). The key elements are:

• Relative time: Time reindexed so that treatment occurs at period 0 for each unit.
• Baseline/reference period: Usually t = -1 (omitted from regression, normalized to zero).
• Pre-trend coefficients: Periods before treatment. If these are close to zero, it supports

the parallel trends assumption.

• Post-treatment coefficients: These capture the dynamic treatment effect over time.

The parallel trends assumption states that, absent treatment, treated and control units would have followed the same trajectory. Pre-trend coefficients near zero are necessary (but not sufficient) evidence for this.

Workflow

For any event study request, follow this general workflow:

1. Understand the data: Identify the panel structure (unit ID, time variable, treatment indicator, outcome). Ask the user if not clear. 2. Determine treatment timing: Is it uniform or staggered? 3. Choose the estimator: Use the decision tree above. 4. Estimate the model: Follow the relevant reference file. 5. Create the plot: Always produce a clean, publication-ready event study plot. 6. Interpret results: Discuss pre-trends, post-treatment dynamics, and any concerns.

Reference Files

Read these as needed based on the user's situation:

• traditional-event-study.md - Traditional TWFE event study with fixest.

Read this for any event study request. It covers the baseline approach, plotting, and the assumptions behind it.

• modern-extensions.md - Modern robust estimators for staggered designs.

Read this when treatment timing varies across units or when the user asks about heterogeneous treatment effects, or any of the newer DiD methods.

• diagnostics-and-testing.md - Pre-trend testing, placebo checks,

Bacon decomposition, sensitivity analysis. Read this when the user wants to validate their design or when you spot potential issues.

Key R Packages

Plotting Standards

All event study plots should follow these defaults (user can override):

• Use ggplot2 with theme_minimal() or a clean custom theme
• Include a vertical dashed line at t = 0 (treatment onset)
• Include a horizontal dashed line at y = 0 (null effect)
• **IMPORTANT: Use discrete point estimates with vertical error bars (geom_point + geom_errorbar),

NOT connected lines with shaded ribbons (geom_line + geom_ribbon).** Each period t should show an individual point with its own error bar. This is the standard format in economics journals.

• Show 95% confidence intervals as error bars (vertical lines at each point)
• Label axes clearly: "Periods Relative to Treatment" (x) and the outcome name (y)
• Normalize the reference period coefficient to 0 explicitly
• Use colorblind-friendly palettes when comparing multiple groups
• If comparing estimators, use facets or distinct colors/shapes with a legend

Common Pitfalls to Watch For

• Forgetting to set the reference period: Always drop one pre-treatment period.
• Binning endpoint periods: With limited pre/post periods, bin the endpoints

(e.g., "-5+" and "5+") to avoid small-sample noise.

• Not-yet-treated as controls: In staggered settings, TWFE uses not-yet-treated

units as controls, which can introduce bias if treatment effects are heterogeneous.

• Interpreting pre-trends as proof of parallel trends: Pre-trends being zero is

necessary but not sufficient. Absence of evidence is not evidence of absence.

• Unbalanced panels: Missing observations can distort event study estimates,

especially at the endpoints.

Diagnostics, Testing, and Sensitivity Analysis

Table of Contents

1. Pre-Trend Testing 2. Bacon Decomposition 3. Placebo and Falsification Tests 4. Sensitivity Analysis with HonestDiD 5. Balance and Covariate Checks 6. Aggregation and Power Considerations

---

1. Pre-Trend Testing

Visual Inspection

The most important diagnostic is the event study plot itself. Look at pre-treatment coefficients: are they close to zero? Do they show a trend? Visual inspection is more informative than formal tests because it reveals patterns (drifts, level shifts) that a single test statistic misses.

Formal Joint Test of Pre-Trends

Test whether all pre-treatment coefficients are jointly zero:

library(fixest)
library(car)

# Estimate the event study
es_model <- feols(
  outcome ~ i(rel_time, ref = -1) | unit_id + year,
  data = df,
  cluster = ~unit_id
)

# Joint F-test of pre-treatment coefficients
# Identify pre-treatment coefficient names
pre_coefs <- grep("rel_time::-[2-9]|rel_time::-[1-9][0-9]", names(coef(es_model)), value = TRUE)

# Using fixest's built-in wald test
wald(es_model, keep = pre_coefs)

Roth (2022) - Pre-Test Bias Warning

Jonathan Roth (2022) "Pre-test with Caution" shows that conditioning on passing a pre-trends test introduces bias. The intuition:

• If you only proceed when pre-trends look flat, you're selecting on a noisy signal.
• This selection can bias your post-treatment estimates, especially when pre-trends

have low power to detect violations.

• The bias is toward finding significant post-treatment effects even when the true

effect is zero.

Implications:

• Don't treat a non-significant pre-trends test as proof of parallel trends.
• Consider the HonestDiD sensitivity analysis (Section 4) as a complement.
• Report the event study plot regardless of whether pre-trends "pass."
• Be especially cautious when the confidence intervals on pre-trends are wide.

---

2. Bacon Decomposition

Goodman-Bacon (2021) showed that TWFE DiD with staggered treatment is a weighted average of all possible 2x2 DiD comparisons. The decomposition reveals which comparisons drive your estimate and whether problematic comparisons (using already-treated as controls) get large weights.

# install.packages("bacondecomp")
library(bacondecomp)

# Requires a balanced panel with a binary treatment indicator
# df needs: unit_id, time, treatment (0/1), outcome

bacon_out <- bacon(
  outcome ~ treatment,
  data = df,
  id_var = "unit_id",
  time_var = "year"
)

# View the decomposition
print(bacon_out)

# The output shows:
# - Type: "Earlier vs Later Treated", "Later vs Earlier Treated", "Treated vs Untreated"
# - Weight: how much each comparison contributes to TWFE
# - Estimate: the 2x2 DiD estimate for each comparison

# Plot the decomposition
ggplot(bacon_out, aes(x = weight, y = estimate, color = type, shape = type)) +
  geom_point(size = 3) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  labs(
    x = "Weight in TWFE Estimate",
    y = "2x2 DiD Estimate",
    title = "Bacon Decomposition of TWFE Estimate",
    color = "Comparison Type",
    shape = "Comparison Type"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "bottom",
    panel.grid.minor = element_blank()
  )

Interpreting the Decomposition

• Treated vs Untreated: Clean comparisons. These are fine.
• Earlier vs Later Treated: Uses not-yet-treated as controls. Fine if no anticipation.
• Later vs Earlier Treated: Uses already-treated as controls. Problematic if

treatment effects change over time.

If "Later vs Earlier Treated" comparisons have large weights and different estimates from other types, your TWFE is likely contaminated. Switch to a robust estimator.

---

3. Placebo and Falsification Tests

Placebo Treatment Timing

Assign a fake treatment date before the actual treatment and test for effects:

# Create a placebo treatment, e.g., 3 years before actual treatment
df_placebo <- df |>
  mutate(
    placebo_treat_year = treat_year - 3,
    placebo_rel_time = year - placebo_treat_year
  ) |>
  # Only use pre-treatment data for this test
  filter(year < treat_year | !treated)

es_placebo <- feols(
  outcome ~ i(placebo_rel_time, ref = -1) | unit_id + year,
  data = df_placebo,
  cluster = ~unit_id
)

iplot(es_placebo,
      main = "Placebo Test: Fake Treatment 3 Years Early",
      xlab = "Periods Relative to Placebo Treatment")
abline(h = 0, lty = 2)

Placebo Outcome

Test whether treatment affects an outcome it theoretically should not affect:

es_placebo_outcome <- feols(
  placebo_outcome ~ i(rel_time, ref = -1) | unit_id + year,
  data = df,
  cluster = ~unit_id
)

iplot(es_placebo_outcome,
      main = "Placebo Test: Effect on Unrelated Outcome")
abline(h = 0, lty = 2)

Randomization / Permutation Inference

Randomly reassign treatment across units and re-estimate. Repeat many times to build a distribution under the null of no effect:

set.seed(58104)
n_perms <- 500
perm_estimates <- numeric(n_perms)

# Get the true estimate (e.g., simple post-treatment ATT)
true_model <- feols(outcome ~ treated_post | unit_id + year, data = df, cluster = ~unit_id)
true_est <- coef(true_model)["treated_post"]

for (i in 1:n_perms) {
  # Randomly reassign treatment groups
  df_perm <- df
  unit_ids <- unique(df_perm$unit_id)
  perm_treated <- sample(unit_ids, sum(df$treated[df$year == min(df$year)]))
  df_perm$treated_perm <- df_perm$unit_id %in% perm_treated
  df_perm$treated_post_perm <- df_perm$treated_perm & (df_perm$year >= df_perm$treat_year)

  perm_model <- feols(outcome ~ treated_post_perm | unit_id + year, data = df_perm)
  perm_estimates[i] <- coef(perm_model)["treated_post_perm"]
}

# Permutation p-value
p_perm <- mean(abs(perm_estimates) >= abs(true_est))
cat("Permutation p-value:", p_perm, "\n")

# Plot the distribution
ggplot(data.frame(est = perm_estimates), aes(x = est)) +
  geom_histogram(bins = 50, fill = "gray70", color = "white") +
  geom_vline(xintercept = true_est, color = "red", linewidth = 1) +
  labs(
    x = "Estimated Effect",
    y = "Count",
    title = "Permutation Distribution",
    caption = paste0("Red line = true estimate. Permutation p-value = ", round(p_perm, 3))
  ) +
  theme_minimal(base_size = 13)

---

4. Sensitivity Analysis with HonestDiD

Rambachan & Roth (2023) "A More Credible Approach to Parallel Trends" provides tools to assess how robust your conclusions are to violations of parallel trends.

The key question: how much would parallel trends need to be violated to explain away your treatment effect?

# install.packages("HonestDiD")
library(HonestDiD)

# First, estimate the event study with fixest
es_model <- feols(
  outcome ~ i(rel_time, ref = -1) | unit_id + year,
  data = df,
  cluster = ~unit_id
)

# Extract results for HonestDiD
# Need: coefficient vector and variance-covariance matrix
# for both pre-treatment and post-treatment coefficients

betahat <- coef(es_model)
sigma <- vcov(es_model)

# Identify pre and post coefficients
# Pre-treatment: negative rel_time (excluding -1 which is the reference)
# Post-treatment: rel_time >= 0

# Get coefficient names and indices
coef_names <- names(betahat)
pre_idx <- grep("rel_time::-[2-9]|rel_time::-[1-9][0-9]", coef_names)
post_idx <- grep("rel_time::[0-9]", coef_names)

# Sort by relative time
pre_times <- as.numeric(stringr::str_extract(coef_names[pre_idx], "-?\\d+"))
post_times <- as.numeric(stringr::str_extract(coef_names[post_idx], "-?\\d+"))

pre_idx <- pre_idx[order(pre_times)]
post_idx <- post_idx[order(post_times)]

# Relative magnitudes approach:
# How large could the post-treatment violation be relative to the max pre-trend?
honest_rm <- createSensitivityResults_relativeMagnitudes(
  betahat = betahat[c(pre_idx, post_idx)],
  sigma = sigma[c(pre_idx, post_idx), c(pre_idx, post_idx)],
  numPrePeriods = length(pre_idx),
  numPostPeriods = length(post_idx),
  Mbarvec = seq(0, 2, by = 0.5)  # Grid of M-bar values
)

# Smoothness-based approach:
# Allows for smooth violations of parallel trends
honest_smooth <- createSensitivityResults(
  betahat = betahat[c(pre_idx, post_idx)],
  sigma = sigma[c(pre_idx, post_idx), c(pre_idx, post_idx)],
  numPrePeriods = length(pre_idx),
  numPostPeriods = length(post_idx),
  Mvec = seq(0, 0.05, by = 0.01)  # Grid of M values (slope change bound)
)

# Original event study plot with sensitivity
createSensitivityPlot_relativeMagnitudes(honest_rm,
  rescaleFactor = 1,
  gridPoints = 100
)

Interpreting HonestDiD

• M-bar = 0: Assumes parallel trends hold exactly. This is the standard result.
• M-bar = 1: Allows post-treatment violations up to the size of the largest pre-trend.
• M-bar = 2: Allows violations twice as large as the largest pre-trend.

If your confidence intervals exclude zero even at M-bar = 1 or 2, your results are robust to meaningful violations of parallel trends. If they don't, your conclusions are more fragile.

---

5. Balance and Covariate Checks

Pre-Treatment Covariate Balance

Check whether treated and control groups are balanced on observables:

library(dplyr)

# Pre-treatment period only
pre_data <- df |> filter(year < treat_year | !treated)

# Balance table
balance <- pre_data |>
  group_by(treated) |>
  summarise(
    across(c(covariate1, covariate2, covariate3),
           list(mean = mean, sd = sd),
           .names = "{.col}_{.fn}")
  )

print(balance)

# Normalized differences (preferred over t-tests for balance)
norm_diff <- function(x_treat, x_control, sd_treat, sd_control) {
  (x_treat - x_control) / sqrt((sd_treat^2 + sd_control^2) / 2)
}

Event Study on Covariates

Run event studies on pre-treatment covariates. If treatment "affects" a covariate that it shouldn't, something is wrong:

# Event study on a covariate
es_covariate <- feols(
  covariate1 ~ i(rel_time, ref = -1) | unit_id + year,
  data = df,
  cluster = ~unit_id
)

iplot(es_covariate,
      main = "Covariate Balance: Event Study on Covariate 1")
abline(h = 0, lty = 2)

---

6. Aggregation and Power Considerations

Power Analysis for Event Studies

Burlig, Preonas, and Woerman (2020) discuss power in event study settings:

# Rough power calculation for event study
# Key parameters:
# - N: number of units
# - T_pre: pre-treatment periods
# - T_post: post-treatment periods
# - prop_treated: share of units treated
# - sigma_unit: unit-level SD
# - rho: within-unit correlation
# - mde: minimum detectable effect

power_es <- function(N, T_pre, T_post, prop_treated, sigma_unit, rho, alpha = 0.05) {
  n_treated <- N * prop_treated
  n_control <- N * (1 - prop_treated)

  # Effective sample size accounting for clustering
  T_total <- T_pre + T_post
  deff <- 1 + (T_total - 1) * rho  # Design effect

  # SE of the DiD estimate (approximate)
  se <- sigma_unit * sqrt(deff * (1/n_treated + 1/n_control) / T_post)

  # MDE for 80% power
  mde <- (qnorm(1 - alpha/2) + qnorm(0.8)) * se

  return(list(se = se, mde = mde))
}

# Example
pow <- power_es(N = 500, T_pre = 5, T_post = 5,
                prop_treated = 0.5, sigma_unit = 1, rho = 0.5)
cat("MDE (80% power):", pow$mde, "\n")

Trimming Extreme Relative Time Periods

When some relative time periods have very few observations, estimates are noisy. Consider trimming or binning:

# Count observations per relative time period
obs_per_period <- df |>
  filter(is.finite(rel_time)) |>
  count(rel_time)

# Only keep periods with sufficient observations
min_obs <- 50  # Adjust based on your data
valid_periods <- obs_per_period |> filter(n >= min_obs) |> pull(rel_time)

es_trimmed <- feols(
  outcome ~ i(rel_time, ref = -1) | unit_id + year,
  data = df |> filter(rel_time %in% valid_periods | !treated),
  cluster = ~unit_id
)

Aggregating Post-Treatment Effects

Often you want a single summary measure of the treatment effect:

# Average post-treatment effect from the event study
post_coefs <- grep("rel_time::[0-9]", names(coef(es_model)), value = TRUE)

# Using fixest's aggregate function
# For overall ATT in sunab models:
# summary(sa_model, agg = "ATT")

# Manual aggregation with lincom:
lincom_formula <- paste0("(", paste(post_coefs, collapse = " + "), ") / ", length(post_coefs))
# Note: exact syntax depends on fixest version

# Or extract and average:
post_est <- coef(es_model)[post_coefs]
post_vcov <- vcov(es_model)[post_coefs, post_coefs]

avg_effect <- mean(post_est)
avg_se <- sqrt(sum(post_vcov)) / length(post_coefs)
cat("Average post-treatment effect:", avg_effect, "(SE:", avg_se, ")\n")

Modern Extensions to Event Studies

Table of Contents

1. Why TWFE Fails with Staggered Treatment 2. Callaway & Sant'Anna (2021) - The did Package 3. Sun & Abraham (2021) - Interaction-Weighted Estimator via fixest 4. Gardner (2022) - Two-Stage DiD via did2s 5. de Chaisemartin & D'Haultfoeuille (2020) - DIDmultiplegt 6. Borusyak, Jaravel & Spiess (2024) - Imputation Estimator via didimputation 7. Roth et al. (2023) Overview and Choosing Between Estimators 8. Comparing Estimators Side by Side

---

1. Why TWFE Fails with Staggered Treatment

The core problem: when treatment rolls out at different times and treatment effects are heterogeneous (vary across cohorts or over time), TWFE can produce biased estimates.

The intuition:

• TWFE implicitly uses already-treated units as controls for newly-treated units.
• This means it's comparing newly-treated units to units who have already been experiencing

treatment for a while.

• If treatment effects change over time (e.g., grow or shrink), this comparison is contaminated.
• The resulting estimates can be a weighted average of treatment effects where some weights

are negative, meaning the overall estimate can even have the wrong sign.

Goodman-Bacon (2021) formalized this by showing that the TWFE estimator is a weighted average of all possible 2x2 DiD comparisons in the data. Some of these comparisons use already-treated units as controls, and the weights on different comparisons can be negative.

The solution: modern estimators that only use valid comparisons - typically comparing newly-treated units to not-yet-treated or never-treated units.

---

2. Callaway & Sant'Anna (2021) - The did Package

Key idea: Estimate group-time average treatment effects ATT(g,t) separately for each treatment cohort g at each time period t, using only clean comparison groups (never-treated or not-yet-treated). Then aggregate these into summary measures.

When to use: This is the most popular and well-tested robust estimator. Good default choice for staggered DiD with heterogeneous effects.

Installation and Setup

# install.packages("did")
library(did)
library(ggplot2)
library(dplyr)

Basic Usage

# Data requirements:
# - yname: outcome variable name (character)
# - tname: time period variable name (character)
# - idname: unit ID variable name (character)
# - gname: group variable = the time period when unit is first treated
#          Use 0 for never-treated units

# Estimate group-time ATTs
cs_out <- att_gt(
  yname = "outcome",
  tname = "year",
  idname = "unit_id",
  gname = "treat_year",  # 0 for never-treated
  data = df,
  control_group = "notyettreated",  # or "nevertreated"
  est_method = "dr",  # doubly robust (default); also "ipw" or "reg"
  clustervars = "unit_id",
  base_period = "varying"  # "varying" (default) or "universal"
)

summary(cs_out)

Aggregation

# Dynamic/event study aggregation (the event study plot)
cs_es <- aggte(cs_out, type = "dynamic")
summary(cs_es)

# Plot the event study
ggdid(cs_es,
      title = "Event Study (Callaway & Sant'Anna)",
      xlab = "Periods Relative to Treatment",
      ylab = "ATT")

# Simple overall ATT aggregation
cs_simple <- aggte(cs_out, type = "simple")
summary(cs_simple)

# Group-specific aggregation (ATT by treatment cohort)
cs_group <- aggte(cs_out, type = "group")
summary(cs_group)

# Calendar time aggregation
cs_calendar <- aggte(cs_out, type = "calendar")
summary(cs_calendar)

Custom ggplot2 Event Study Plot

# Extract event study coefficients for custom plotting
cs_df <- data.frame(
  rel_time = cs_es$egt,
  estimate = cs_es$att.egt,
  se = cs_es$se.egt
) |>
  mutate(
    conf.low = estimate - 1.96 * se,
    conf.high = estimate + 1.96 * se
  )

ggplot(cs_df, aes(x = rel_time, y = estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_vline(xintercept = -0.5, linetype = "dashed", color = "gray50") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2, color = "#E69F00") +
  geom_point(color = "#E69F00", size = 2) +
  labs(
    x = "Periods Relative to Treatment",
    y = "ATT",
    title = "Event Study (Callaway & Sant'Anna 2021)",
    caption = "95% pointwise confidence intervals. Control group: not-yet-treated."
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )

Options and Variants

# With covariates (must be time-invariant or pre-treatment values)
cs_out_cov <- att_gt(
  yname = "outcome",
  tname = "year",
  idname = "unit_id",
  gname = "treat_year",
  data = df,
  xformla = ~ covariate1 + covariate2,  # Covariate formula
  control_group = "notyettreated",
  est_method = "dr"
)

# Using never-treated only as controls
cs_out_nt <- att_gt(
  yname = "outcome",
  tname = "year",
  idname = "unit_id",
  gname = "treat_year",
  data = df,
  control_group = "nevertreated"
)

# Simultaneous confidence bands (uniform inference, more conservative)
cs_es_unif <- aggte(cs_out, type = "dynamic")
# ggdid shows both pointwise and simultaneous CIs by default
ggdid(cs_es_unif)

---

3. Sun & Abraham (2021) - Interaction-Weighted Estimator

Key idea: Decompose the TWFE event study into cohort-specific estimates, then aggregate using the right weights (cohort shares). Implemented directly in fixest via sunab().

When to use: When you like the fixest workflow and want a drop-in robust replacement for the traditional event study. Requires a never-treated or last-treated group.

Implementation

library(fixest)

# sunab(treat_year, year) handles everything:
# - cohort_var: the period each unit was first treated (Inf or NA for never-treated)
# - period_var: the calendar time
sa_model <- feols(
  outcome ~ sunab(treat_year, year) | unit_id + year,
  data = df,
  cluster = ~unit_id
)

summary(sa_model)

# Event study plot (fixest native)
iplot(sa_model,
      xlab = "Periods Relative to Treatment",
      ylab = "Estimated Effect",
      main = "Event Study (Sun & Abraham 2021)")
abline(h = 0, lty = 2, col = "gray50")
abline(v = -0.5, lty = 2, col = "gray50")

# Get the aggregated ATT
summary(sa_model, agg = "ATT")

# Aggregate by relative period (cohort-averaged dynamic effects)
summary(sa_model, agg = "period")

ggplot2 Custom Plot

# Extract coefficients
sa_coefs <- broom::tidy(sa_model, conf.int = TRUE) |>
  mutate(
    rel_time = as.numeric(stringr::str_extract(term, "-?\\d+"))
  ) |>
  filter(!is.na(rel_time))

# Add reference period
sa_coefs <- bind_rows(
  sa_coefs,
  data.frame(term = "ref", estimate = 0, std.error = 0,
             conf.low = 0, conf.high = 0, rel_time = -1)
)

ggplot(sa_coefs, aes(x = rel_time, y = estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_vline(xintercept = -0.5, linetype = "dashed", color = "gray50") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2, color = "#56B4E9") +
  geom_point(color = "#56B4E9", size = 2) +
  labs(
    x = "Periods Relative to Treatment",
    y = "Estimated Effect",
    title = "Event Study (Sun & Abraham 2021)"
  ) +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"), panel.grid.minor = element_blank())

---

4. Gardner (2022) - Two-Stage DiD via did2s

Key idea: In stage 1, estimate unit and time fixed effects using only untreated observations (not-yet-treated and never-treated in their pre-treatment periods). In stage 2, residualize the outcome for all observations, then regress the residualized outcome on treatment indicators. This "cleans out" the fixed effects without contamination from treatment effects.

When to use: Conceptually simple, handles staggered treatment well. Good complement to Callaway & Sant'Anna.

Implementation

# install.packages("did2s")
library(did2s)

# Basic usage
g_model <- did2s(
  data = df,
  yname = "outcome",
  first_stage = ~ 0 | unit_id + year,  # Fixed effects (first stage)
  second_stage = ~ i(rel_time, ref = -1),  # Event study (second stage)
  treatment = "treated_post",  # Binary: 1 if treated and in post-treatment period
  cluster_var = "unit_id"
)

# fixest-style plotting works directly
iplot(g_model,
      xlab = "Periods Relative to Treatment",
      ylab = "Estimated Effect",
      main = "Event Study (Gardner 2022, did2s)")

# Or with ggplot2 (same extraction method as fixest models)

With Controls

g_model_cov <- did2s(
  data = df,
  yname = "outcome",
  first_stage = ~ covariate1 + covariate2 | unit_id + year,
  second_stage = ~ i(rel_time, ref = -1),
  treatment = "treated_post",
  cluster_var = "unit_id"
)

---

5. de Chaisemartin & D'Haultfoeuille (2020) - DIDmultiplegt

Key idea: Estimate treatment effects period by period, comparing newly-treated units (switchers) to units whose treatment status doesn't change (stable units). Uses only "clean" comparisons.

When to use: When you have a binary treatment that switches on (and possibly off). Especially useful when treatment can turn on and off.

Implementation

# install.packages("DIDmultiplegt")
library(DIDmultiplegt)

# Basic estimation
dcdh_out <- did_multiplegt(
  df = df,
  Y = "outcome",
  G = "unit_id",
  T = "year",
  D = "treatment",  # Binary treatment indicator (0/1)
  placebo = 4,  # Number of pre-treatment placebo periods to estimate
  dynamic = 4,  # Number of post-treatment dynamic effects to estimate
  brep = 100,   # Bootstrap replications for inference
  cluster = "unit_id"
)

# The output is a list with effects, placebos, and SEs
# Create a manual event study plot:
n_placebo <- 4
n_dynamic <- 4

dcdh_df <- data.frame(
  rel_time = c(-(n_placebo:1), 0:n_dynamic),
  estimate = c(
    rev(dcdh_out$placebo),   # Placebo estimates (reverse order)
    dcdh_out$effect,          # Dynamic effects (including contemporaneous)
    dcdh_out$dynamic
  ),
  se = c(
    rev(dcdh_out$se_placebo),
    dcdh_out$se_effect,
    dcdh_out$se_dynamic
  )
) |>
  mutate(
    conf.low = estimate - 1.96 * se,
    conf.high = estimate + 1.96 * se
  )

ggplot(dcdh_df, aes(x = rel_time, y = estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_vline(xintercept = -0.5, linetype = "dashed", color = "gray50") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2, color = "#009E73") +
  geom_point(color = "#009E73", size = 2) +
  labs(
    x = "Periods Relative to Treatment",
    y = "Estimated Effect",
    title = "Event Study (de Chaisemartin & D'Haultfoeuille 2020)"
  ) +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"), panel.grid.minor = element_blank())

Note: The DIDmultiplegt package has been updated in recent versions. Check the package documentation for the latest syntax. A newer version DIDmultiplegtDYN is also available with additional features.

---

6. Borusyak, Jaravel & Spiess (2024) - Imputation Estimator

Key idea: Impute untreated potential outcomes for treated units using a model estimated on untreated observations (similar to Gardner's first stage), then compute treatment effects as the difference between observed and imputed outcomes. Efficient and handles heterogeneity well.

When to use: Good all-around choice, especially for large panels. Computationally efficient.

Implementation

# install.packages("didimputation")
library(didimputation)

# Basic event study
bjs_model <- did_imputation(
  data = df,
  yname = "outcome",
  gname = "treat_year",  # First treatment period (0 for never-treated)
  tname = "year",
  idname = "unit_id",
  horizon = c(-5:5),  # Relative time periods to estimate
  pretrends = TRUE     # Include pre-treatment periods
)

# Plot
bjs_model |>
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_vline(xintercept = -0.5, linetype = "dashed", color = "gray50") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2, color = "#CC79A7") +
  geom_point(color = "#CC79A7", size = 2) +
  labs(
    x = "Periods Relative to Treatment",
    y = "Estimated Effect",
    title = "Event Study (Borusyak, Jaravel & Spiess 2024)"
  ) +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"), panel.grid.minor = element_blank())

---

7. Roth et al. (2023) - Choosing Between Estimators

Roth, Sant'Anna, Bilinski, and Poe (2023) "What's Trending in Difference-in-Differences?" provide an excellent overview. Key takeaways for choosing:

Practical advice:

• For most applied work, start with Callaway & Sant'Anna (most flexible, well-tested).
• If you're already using fixest, Sun & Abraham via sunab() is a natural choice.
• Run two or three estimators and check that results are qualitatively similar.
• If results differ substantially across estimators, investigate why - it may reveal

important features of your setting.

---

8. Comparing Estimators Side by Side

Here's a template for plotting multiple estimators together:

library(ggplot2)
library(dplyr)

# Assume you have data frames from each estimator with columns:
# rel_time, estimate, conf.low, conf.high, estimator

all_estimates <- bind_rows(
  twfe_df |> mutate(estimator = "TWFE"),
  cs_df |> mutate(estimator = "Callaway & Sant'Anna"),
  sa_df |> mutate(estimator = "Sun & Abraham")
)

# Color palette (colorblind-friendly)
colors <- c(
  "TWFE" = "#0072B2",
  "Callaway & Sant'Anna" = "#E69F00",
  "Sun & Abraham" = "#56B4E9"
)

ggplot(all_estimates, aes(x = rel_time, y = estimate, color = estimator)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_vline(xintercept = -0.5, linetype = "dashed", color = "gray50") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2,
                position = position_dodge(width = 0.3)) +
  geom_point(size = 2, position = position_dodge(width = 0.3)) +
  scale_color_manual(values = colors) +
  labs(
    x = "Periods Relative to Treatment",
    y = "Estimated Effect",
    title = "Event Study: Comparing Estimators",
    color = "Estimator"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

# Alternative: faceted
ggplot(all_estimates, aes(x = rel_time, y = estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_vline(xintercept = -0.5, linetype = "dashed", color = "gray50") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2, color = "steelblue") +
  geom_point(color = "steelblue", size = 1.5) +
  facet_wrap(~estimator, ncol = 1) +
  labs(
    x = "Periods Relative to Treatment",
    y = "Estimated Effect",
    title = "Event Study: Comparing Estimators"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.minor = element_blank(),
    strip.text = element_text(face = "bold")
  )

Traditional Event Study with TWFE

The Setup

A traditional event study uses a two-way fixed effects (TWFE) regression with leads and lags of a treatment indicator. The idea: if you can show that outcomes evolved similarly between treated and control groups before treatment (the "pre-trends"), and then diverged after, you have compelling evidence that the treatment caused the divergence.

The Model

The canonical specification is:

Y_it = α_i + λ_t + Σ_{k ≠ -1} β_k × D_it^k + X_it'γ + ε_it

Where:

• Y_it is the outcome for unit i at time t
• α_i are unit fixed effects
• λ_t are time fixed effects
• D_it^k is an indicator for unit i being k periods from treatment at time t
• k = -1 is omitted (the reference period)
• X_it are optional time-varying controls
• β_k are the event study coefficients: pre-trends for k < -1, treatment effects for k ≥ 0

Assumptions

1. Parallel trends: Absent treatment, treated and control units would have followed the same outcome trajectory. This is the key identifying assumption and is untestable in the post-treatment period. Pre-trends provide indirect evidence.

2. No anticipation: Units don't change behavior before treatment actually occurs. If violated, you can adjust the reference period or allow for anticipation effects.

3. SUTVA (Stable Unit Treatment Value Assumption): One unit's treatment doesn't affect another unit's outcome (no spillovers), and treatment is well-defined.

4. Homogeneous treatment effects (for TWFE validity in staggered designs): When treatment timing varies, TWFE assumes treatment effects are the same across cohorts and over time since treatment. If this fails, TWFE estimates can be biased - this is the key motivation for modern robust estimators.

Implementation in R with fixest

fixest is the recommended package for traditional event studies. It's fast, has excellent syntax, and handles event studies natively with i() for interaction terms and sunab() for Sun & Abraham.

Basic Event Study

library(fixest)
library(ggplot2)

# ---- Data Preparation ----
# Assume: df has columns: unit_id, year, treat_year (year of treatment, Inf or NA if never treated), outcome
# Create relative time variable
df$rel_time <- df$year - df$treat_year

# For never-treated units, set rel_time to a large negative or handle via NA
# fixest handles never-treated units well if treat_year is set appropriately

# ---- Estimation ----
# Using fixest::i() to create the event study dummies
# ref = -1 sets the reference period
# bin argument bins extreme values (e.g., bin = c(-5, 5) bins ≤-5 and ≥5)

es_model <- feols(
  outcome ~ i(rel_time, ref = -1) | unit_id + year,
  data = df,
  cluster = ~unit_id
)

summary(es_model)

Plotting with fixest's Built-in Tools

# fixest has a native event study plotting function
iplot(es_model,
      xlab = "Periods Relative to Treatment",
      ylab = "Estimated Effect",
      main = "Event Study: Effect on Outcome")

# Add reference lines manually if needed
abline(h = 0, lty = 2, col = "gray50")
abline(v = -0.5, lty = 2, col = "gray50")

Publication-Quality Plot with ggplot2

For more control over the plot, extract coefficients and use ggplot2:

library(ggplot2)
library(broom)

# Extract coefficients
es_coefs <- broom::tidy(es_model, conf.int = TRUE) |>
  # Parse the relative time from coefficient names
  # fixest names them like "rel_time::-5", "rel_time::0", etc.
  dplyr::mutate(
    rel_time = as.numeric(stringr::str_extract(term, "-?\\d+"))
  ) |>
  dplyr::filter(!is.na(rel_time))

# Add the reference period (coefficient = 0 by construction)
ref_row <- data.frame(
  term = "ref",
  estimate = 0,
  std.error = 0,
  conf.low = 0,
  conf.high = 0,
  rel_time = -1
)
es_coefs <- dplyr::bind_rows(es_coefs, ref_row)

# Plot
ggplot(es_coefs, aes(x = rel_time, y = estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_vline(xintercept = -0.5, linetype = "dashed", color = "gray50") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2, color = "steelblue") +
  geom_point(color = "steelblue", size = 2) +
  labs(
    x = "Periods Relative to Treatment",
    y = "Estimated Effect on Outcome",
    title = "Event Study",
    caption = "Notes: Reference period is t = -1. 95% confidence intervals shown.\nStandard errors clustered at the unit level."
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )

With Binned Endpoints

When you have very long pre or post periods, bin the endpoints to avoid noisy estimates from periods with few observations:

# Bin at -5 and +5
es_model_binned <- feols(
  outcome ~ i(rel_time, ref = -1, bin = c(-5, 5)) | unit_id + year,
  data = df,
  cluster = ~unit_id
)

iplot(es_model_binned)

With Controls

# Add time-varying controls
es_model_controls <- feols(
  outcome ~ i(rel_time, ref = -1) + control_var1 + control_var2 | unit_id + year,
  data = df,
  cluster = ~unit_id
)

With Heterogeneity by Group

# Estimate separate event studies by group (e.g., by region)
es_by_group <- feols(
  outcome ~ i(rel_time, ref = -1) | unit_id + year,
  data = df,
  cluster = ~unit_id,
  split = ~group_var
)

# Plot all together
iplot(es_by_group)

When Traditional TWFE Breaks Down

The traditional event study with TWFE works well when:

• All treated units are treated at the same time, OR
• Treatment effects are homogeneous across cohorts and over time

It can produce misleading estimates when:

• Treatment is staggered AND effects vary across cohorts or evolve differently over time
• Already-treated units serve as "controls" for later-treated units
• Negative weights emerge in the TWFE estimand (see Bacon decomposition in diagnostics)

When these concerns apply, move to the modern extensions (see modern-extensions.md).

Simulated Example for Testing

Here's a complete simulated example you can use to verify code works:

library(fixest)
library(ggplot2)
library(dplyr)

set.seed(90841)

# Simulate panel data
n_units <- 200
n_periods <- 20
treat_period <- 11  # Treatment happens at period 11 for treated units

df_sim <- expand.grid(
  unit_id = 1:n_units,
  time = 1:n_periods
) |>
  mutate(
    treated = unit_id <= 100,  # First 100 units are treated
    treat_year = ifelse(treated, treat_period, Inf),
    rel_time = time - treat_year,
    # True effect: grows linearly after treatment
    true_effect = ifelse(treated & time >= treat_period, (time - treat_period + 1) * 0.5, 0),
    # Unit and time fixed effects + noise
    unit_fe = rep(rnorm(n_units, 0, 2), each = n_periods),
    time_fe = rep(rnorm(n_periods, 0, 1), times = n_units),
    outcome = unit_fe + time_fe + true_effect + rnorm(n_units * n_periods, 0, 1)
  )

# Estimate event study
es <- feols(
  outcome ~ i(rel_time, ref = -1) | unit_id + time,
  data = df_sim |> filter(is.finite(rel_time) | !treated),
  cluster = ~unit_id
)

# Plot
iplot(es,
      xlab = "Periods Relative to Treatment",
      ylab = "Estimated Effect",
      main = "Simulated Event Study")
abline(h = 0, lty = 2, col = "gray50")
abline(v = -0.5, lty = 2, col = "gray50")