tutorial120 minutes

Lab: Staggered Difference-in-Differences

Tackle the staggered DiD problem step by step. Learn why two-way fixed effects fails with heterogeneous treatment effects, decompose TWFE with Bacon, and estimate robust alternatives using Callaway-Sant'Anna and Sun-Abraham.

Overview

In this lab you will estimate the effect of a state-level policy that is adopted at different times across states. You will discover that the standard two-way fixed effects (TWFE) estimator can be severely biased when treatment effects vary over time, and you will implement modern robust alternatives. This development represents a significant recent methodological advance in applied microeconometrics.

What you will learn:

Why standard TWFE can give negative estimates even when all treatment effects are positive
How to decompose the TWFE estimator using the Bacon decomposition
How to estimate robust group-time average treatment effects with Callaway and Sant'Anna (2021)
How to estimate interaction-weighted estimators with Sun and Abraham (2021)
How to compare and interpret results across methods

Prerequisites: Familiarity with difference-in-differences and two-way fixed effects (see the DiD and FE labs).

Step 1: Simulate Staggered Adoption Data

We create a panel of 40 states over 20 years. States adopt a policy in different years, and the treatment effect grows over time since adoption (dynamic treatment effects). A control group of states never adopts.

1library(fixest)
2library(did)
3library(bacondecomp)
4library(modelsummary)
5
6set.seed(42)
7n_states <- 40
8n_years <- 20
9years <- 2000:(2000 + n_years - 1)
10
11adoption <- c(rep(0, 10), rep(2005, 8), rep(2010, 8),
12            rep(2015, 7), rep(2017, 7))
13
14state_fe <- rnorm(n_states, sd = 2)
15year_fe <- seq(0, 3, length.out = n_years)
16
17df <- expand.grid(state = 1:n_states, year = years)
18df$state_fe <- state_fe[df$state]
19df$year_fe <- year_fe[df$year - 1999]
20df$adoption_year <- adoption[df$state]
21df$treat <- as.integer(df$year >= df$adoption_year & df$adoption_year > 0)
22df$years_since <- ifelse(df$treat == 1, df$year - df$adoption_year, 0)
23df$tau_true <- ifelse(df$treat == 1, 0.5 + 0.1 * df$years_since, 0)
24df$y <- df$state_fe + df$year_fe + df$tau_true + rnorm(nrow(df))
25df$cohort <- ifelse(df$adoption_year > 0,
26                   as.character(df$adoption_year), "Never")
27df$state <- factor(df$state)
28df$year_f <- factor(df$year)
29
30cat("True ATT:", mean(df$tau_true[df$treat == 1]), "\n")

Requiresfixest did bacondecomp modelsummary

Expected output:

state	year	y	treat	adoption_year	cohort	tau_true
0	2000	-1.5	0	0	Never	0.0
0	2001	-1.3	0	0	Never	0.0
10	2005	2.1	1	2005	2005	0.5
10	2006	2.9	1	2005	2005	0.6
10	2010	3.8	1	2005	2005	1.0

Panel: 40 states x 20 years = 800 obs

Cohort sizes:
  Never: 10 states
  2005:   8 states (early adopters)
  2010:   8 states (middle adopters)
  2015:   7 states (late adopters)
  2017:   7 states (very late adopters)

True ATT (all treated): ~1.05

The true ATT varies because treatment effects grow over time (tau = 0.5 + 0.1 * years_since_treatment). Early adopters have accumulated larger effects by the end of the panel.

Step 2: The TWFE Regression (and Its Problems)

1# Standard TWFE
2m_twfe <- feols(y ~ treat | state + year_f, data = df, vcov = ~state)
3
4true_att <- mean(df$tau_true[df$treat == 1])
5cat("=== Standard TWFE ===\n")
6cat("TWFE estimate:", coef(m_twfe)["treat"], "\n")
7cat("True ATT:", true_att, "\n")
8cat("Bias:", coef(m_twfe)["treat"] - true_att, "\n")

Expected output:

=== Standard TWFE ===
TWFE estimate:    ~0.65
True ATT:         ~1.05
Bias:             ~-0.40
SE (clustered):   ~0.10

The TWFE estimate is biased because treatment effects
are heterogeneous across time. TWFE puts negative weight
on some group-time treatment effects.

The TWFE estimate (~0.65) substantially underestimates the true ATT (~1.05). This underestimation occurs because TWFE implicitly uses early-adopter states (whose effects have grown large) as controls for late-adopter states, generating "forbidden comparisons" that pull the estimate downward.

Concept Check

The true average treatment effect on the treated is about 1.0, but your TWFE estimate is only 0.6. What is causing this underestimation?

The sample size is too small for the fixed effects to work properly.TWFE implicitly uses already-treated units (whose effects have grown larger) as controls for newly-treated units. These 'bad comparisons' assign negative weights to some treatment effects, biasing the overall estimate.The state fixed effects absorb too much variation.You need to add state-specific time trends.

Step 3: Bacon Decomposition

The Bacon decomposition ((Goodman-Bacon, 2021)) shows exactly which 2x2 comparisons make up the TWFE estimate and what weights they receive.

1# Bacon decomposition using bacondecomp package
2# Requires numeric treatment timing variable
3df_bacon <- df[, c("state", "year", "y", "adoption_year")]
4df_bacon$state <- as.numeric(df_bacon$state)
5df_bacon$D <- ifelse(df_bacon$adoption_year > 0, df_bacon$adoption_year, 10000)
6
7bacon_out <- bacon(y ~ treat,
8                  data = df[df$adoption_year > 0 | df$adoption_year == 0, ],
9                  id_var = "state",
10                  time_var = "year")
11
12# If bacon() requires specific format:
13# Show the decomposition
14print(bacon_out)
15
16# Weighted average should equal TWFE
17cat("\nWeighted average:", sum(bacon_out$estimate * bacon_out$weight), "\n")
18cat("TWFE estimate:", coef(m_twfe)["treat"], "\n")

Requiresbacondecomp

Expected output:

Comparison	2x2 DiD Estimate	True ATT for Cohort
Cohort 2005 vs Never	~1.25	~1.25
Cohort 2010 vs Never	~0.75	~0.75
Cohort 2015 vs Never	~0.60	~0.60
Cohort 2017 vs Never	~0.50	~0.50

2x2 DiD estimates by comparison type:
  Cohort 2005 vs Never: est = ~1.25, true ATT = ~1.25
  Cohort 2010 vs Never: est = ~0.75, true ATT = ~0.75
  Cohort 2015 vs Never: est = ~0.60, true ATT = ~0.60
  Cohort 2017 vs Never: est = ~0.50, true ATT = ~0.50

TWFE combines all these comparisons — including using
already-treated cohorts as controls for late adopters.

Clean comparisons (treated vs. never-treated) recover accurate estimates. The bias in TWFE arises from the additional comparisons where early adopters serve as controls for late adopters: their growing treatment effects create a "moving baseline" that contaminates the estimate.

Step 4: Callaway and Sant'Anna (2021)

The Callaway-Sant'Anna estimator computes group-time ATTs — the treatment effect for each cohort at each post-treatment period — using only clean comparisons (never-treated or not-yet-treated as controls).

1# Callaway-Sant'Anna using the did package
2df$id <- as.numeric(df$state)
3df$G <- ifelse(df$adoption_year > 0, df$adoption_year, 0)
4
5cs_out <- att_gt(
6yname = "y",
7tname = "year",
8idname = "id",
9gname = "G",
10data = as.data.frame(df),
11control_group = "nevertreated"
12)
13
14# Summary
15summary(cs_out)
16
17# Aggregate to simple ATT
18cs_agg <- aggte(cs_out, type = "simple")
19summary(cs_agg)
20
21cat("\nCS simple ATT:", cs_agg$overall.att, "\n")
22cat("TWFE:", coef(m_twfe)["treat"], "\n")
23cat("True ATT:", true_att, "\n")
24
25# Event study aggregation
26cs_es <- aggte(cs_out, type = "dynamic")
27ggdid(cs_es)

Requiresdid

Expected output:

Callaway-Sant'Anna (simple aggregation):
  Estimated ATT: ~1.05
  True ATT:      ~1.05
  TWFE estimate: ~0.65

CS is much closer to the truth than TWFE.

Expected output: Group-time ATTs (selected)

Cohort	Year	ATT(g,t)	True Effect
2005	2005	~0.50	0.50
2005	2008	~0.80	0.80
2005	2015	~1.50	1.50
2010	2010	~0.50	0.50
2010	2015	~1.00	1.00
2015	2015	~0.50	0.50
2017	2017	~0.50	0.50

Step 5: Sun and Abraham (2021)

Sun and Abraham (2021) propose an interaction-weighted estimator that corrects TWFE by interacting cohort indicators with relative time indicators.

1# Sun and Abraham using fixest::sunab()
2df$G_sa <- ifelse(df$adoption_year > 0, df$adoption_year, 10000)
3
4m_sa <- feols(y ~ sunab(G_sa, year) | state + year_f,
5            data = df, vcov = ~state)
6summary(m_sa)
7
8# Event study plot
9iplot(m_sa, main = "Sun-Abraham Event Study")
10
11# Aggregate ATT
12cat("\nSun-Abraham aggregate ATT:", summary(m_sa, agg = "ATT")$coeftable[1], "\n")
13cat("Callaway-Sant'Anna ATT:", cs_agg$overall.att, "\n")
14cat("TWFE:", coef(m_twfe)["treat"], "\n")
15cat("True:", true_att, "\n")

Requiresfixest

Expected output: Sun-Abraham event study

Relative Time (e)	SA Estimate	True Effect
e = 0	~0.50	0.50
e = +1	~0.60	0.60
e = +2	~0.70	0.70
e = +3	~0.80	0.80
e = +5	~1.00	1.00
e = +7	~1.20	1.20

=== Sun-Abraham Event Study ===
  e = +0: SA estimate = ~0.50, True = 0.50
  e = +1: SA estimate = ~0.60, True = 0.60
  e = +2: SA estimate = ~0.70, True = 0.70
  e = +3: SA estimate = ~0.80, True = 0.80
  ...

The Sun-Abraham estimates correctly recover the dynamic treatment effects, showing the linear growth pattern from the DGP. Each estimate is a cohort-size-weighted average of the cohort-specific effects at that relative time.

Concept Check

Both Callaway-Sant'Anna and Sun-Abraham produce estimates close to the true ATT, while TWFE is biased. What is the key difference in how these estimators handle staggered treatment?

They use different fixed effects.They avoid using already-treated units as controls. Instead, they either use never-treated (or not-yet-treated) units as the comparison group, ensuring that treatment effect dynamics in already-treated units do not contaminate the estimates.They use a different outcome variable.They require larger sample sizes.

Step 6: Compare All Estimates

1# Summary comparison
2cat("=" , rep("=", 50), "\n")
3cat("Method", "\t\t\t", "Estimate", "\t", "True\n")
4cat("TWFE:", "\t\t\t", round(coef(m_twfe)["treat"], 3), "\t",
5  round(true_att, 3), "\n")
6cat("Callaway-Sant'Anna:", "\t", round(cs_agg$overall.att, 3), "\t",
7  round(true_att, 3), "\n")
8cat("Sun-Abraham:", "\t\t",
9  round(summary(m_sa, agg = "ATT")$coeftable[1], 3), "\t",
10  round(true_att, 3), "\n")
11
12# Side-by-side event study
13# CS dynamic
14cs_es <- aggte(cs_out, type = "dynamic")
15ggdid(cs_es) + ggtitle("Callaway-Sant'Anna Event Study")

Expected output:

Method	Estimate	True ATT	Bias	Bias (%)
TWFE	~0.65	~1.05	~-0.40	~-38%
Callaway-Sant'Anna	~1.05	~1.05	~0.00	~0%
Sun-Abraham	~1.05	~1.05	~0.00	~0%

Step 7: Pre-Trends and Diagnostics

1# Pre-trends from CS event study
2cs_es <- aggte(cs_out, type = "dynamic", min_e = -5, max_e = 10)
3summary(cs_es)
4
5# Visual check
6ggdid(cs_es) +
7ggtitle("Event Study with Pre-Trends") +
8geom_hline(yintercept = 0, linetype = "dashed")

Expected output: Pre-trends check

Relative Time	Pre-Treatment ATT	Significant?
e = -3	~0.02	No
e = -2	~-0.01	No

=== Pre-Trends Check ===
  e = -3: ATT = ~0.02
  e = -2: ATT = ~-0.01

Pre-treatment estimates should be close to zero.

All pre-treatment estimates are statistically indistinguishable from zero, supporting the parallel trends assumption. This pattern is expected because treatment is not anticipated in the DGP.

Step 8: Exercises

Try these on your own:

Constant treatment effects. Modify the simulation so that tau = 0.5 for all cohorts and all post-periods (no dynamics). Show that TWFE recovers the correct estimate in this case.
Not-yet-treated controls. Re-estimate the CS model using control_group = "notyettreated". Compare with the never-treated control group. When might this matter?
Rambachan-Roth sensitivity. Install HonestDiD (R) and assess how sensitive your results are to possible violations of parallel trends.
de Chaisemartin and D'Haultfoeuille. Estimate the treatment effect using their did_multiplegt estimator. Compare with CS and SA.
Heterogeneity analysis. Modify the simulation so that early adopters have larger treatment effects than late adopters. How does this additional heterogeneity affect the TWFE bias?

Expected Output

If your code runs correctly, you should see:

True ATT (all treated observations): Around 0.8–1.5, depending on the mix of cohorts and time since treatment (DGP: 0.5 + 0.1 * years_since)
TWFE static estimate: A single coefficient that may underestimate (or even be negative) the true ATT due to "forbidden comparisons" using already-treated units as controls
Bacon decomposition: Reveals negative weights on some 2x2 comparisons, particularly "later vs. earlier" comparisons where already-treated units serve as controls
Callaway-Sant'Anna group-time ATTs: Positive and growing over time for each cohort (e.g., 0.5 at k=0, 0.6 at k=1, etc.), matching the true DGP
Sun-Abraham interaction-weighted estimates: Similar to CS, correctly recovering the dynamic treatment effects
CS aggregated ATT: Close to the true average ATT, unlike the TWFE estimate
Event study plot (robust estimators): Flat pre-trends and growing post-treatment effects
Panel dimensions: 40 states x 20 years = 800 observations, with 4 adoption cohorts (2005, 2010, 2015, 2017) and 10 never-treated

Summary

In this lab you learned:

Standard TWFE can be severely biased in staggered designs when treatment effects vary across cohorts or over time
The Bacon decomposition reveals that TWFE is a weighted average of all 2x2 DiD comparisons, some of which use already-treated units as controls
Callaway and Sant'Anna (2021) estimate clean group-time ATTs using never-treated or not-yet-treated units as controls
Sun and Abraham (2021) correct TWFE by using interaction-weighted estimation with cohort-specific relative time effects
Both robust estimators recover the true dynamic treatment effects, while TWFE produces a misleading single number
Reporting event study plots with pre-trends is increasingly common in practice, and being transparent about the comparison between TWFE and robust estimators strengthens credibility