Difference-in-Differences

Module 3.2: Staggered Treatment Timing

Author

Affiliation

Alex Cardazzi

Old Dominion University

All materials can be found at alexcardazzi.github.io.

Staggered Treatment Timing

Difference-in-Differences is one of the most widely used research designs in the social sciences. This is likely due to how (relatively) easy it is to set up and to understand. While economists tend to prefer designs like Regression Discontinuity (Module 7) in terms of credibility, DiD remains a much more popular method. There have been many econometric innovations in the past decade or so, but none as big as the DiD one. Certainly, no other method has had a “revolution” in the way DiD has.

In these notes, we are going to discuss the innovations that have come to light since ~2020.

When treatment timing is staggered (for example: policies are rolled out in different places at different times), and treatment effects are heterogeneous (for example: different across time/units), the DiD coefficient occasionally fails to represents what we think it should represent. In short, this is due to improper comparisons of groups at different points in time.

In a TWFE design, OLS¹ compares the following:

• Newly treated units to never treated units (good!)
• Newly treated units to not-yet treated units (good!)
• Newly treated units to already treated units (bad!)

This last group, where newly treated units are compared to already treated units, can result in negative estimates of positive treatment effects and vice versa. This issue puzzled social scientists and ubiquitously stalled research. For a few years, economists everywhere were scrambling to figure out what was wrong with two-way fixed effects models, who Andrew Goodman-Bacon was, and how to do his difference-in-differences thing. This was such a FAQ that in late 2019, Goodman-Bacon wrote a paper So You’ve Been Told to Do My Difference-in-Differences Thing: A Guide.

Let’s spend some time motivating and understanding the problem with TWFE before moving onto the shortcuts some solutions to the problem.

TWFE

Consider the following example that we will simulate data for. We have three groups over ten time periods. The first group is never treated throughout the sample. The second and third groups are treated at \(t = 4.5\) and \(t = 7.5\), but we only observe the units at \(t \in \{1, 2, ..., 10\}\). Put more simply: at \(t = 4\), the second group is 100% untreated and at \(t = 5\), the group is 100% treated. The treatment is either going to change the trajectory of growth for each group, or simply be a static level shift, depending on the arguments you choose for the simulation. You should experiment with both, however.

Below are the simulation functions that create the data and then plot the data. I am providing the code so you can see what it looks like under the hood.²

create_data <- function(dynamic = FALSE, rnd = FALSE){
  
  set.seed(757)
  df <- data.frame(id = rep(1:3, times = 10),
                   tt = rep(1:10, each = 3))
  df$col <- ifelse(df$id == 1, "black",
                   ifelse(df$id == 2, "tomato",
                          "dodgerblue"))
  
  # Setting the binary treatment variable
  df$D <- ifelse((df$id == 2 & df$tt >= 5) | (df$id == 3 & df$tt >= 8), 1, 0)
  
  df$btrue <- ifelse(df$D == 1 & df$id == 3, 5,
                     ifelse(df$D == 1 & df$id == 2, 3, 0))
  # Decide if treatment is dynamic or static.
  if(dynamic){
    
    df$rel_time <- df$tt - c(NA, 5, 8)
    df$rel_time <- ifelse(is.na(df$rel_time), -100, df$rel_time)
  } else {
    
    df$rel_time <- 1
  }
  
  # Setting the "true" treatment effect for each time period
  #   This will be nice for later comparison
  df$treatment_effect <- df$btrue*df$D*df$rel_time
  
  # Create outcome variable
  df$y <- df$id + df$tt + df$treatment_effect + if(rnd) rnorm(nrow(df), 0, 0.01) else 0
  
  return(df)
}
plot_data <- function(df){
  
  plot(df$tt, df$y, col = df$col, pch = 19,
       xlab = "Time", ylab = "Y", las = 1)
  for(i in 1:3){
    
    lines(df$tt[df$id == i],
          df$y[df$id == i],
          col = df$col[df$id == i])
  }
  abline(v = c(4.5, 7.5), lty = 2)
}

Compare the static and dynamic data by switching between TRUE and FALSE below:

Since we simulated these data, and we know the true treatment effect for each observation (via treatment_effect), we should be able to average them across treated units. Of course, since the treatments are inherently different for groups 2 and 3, there are many ways we can summarize these treatment effects. For example, since the two treatment effects are 3 and 5 in the static case, we can use the average (4) for our treatment effect. We could also weigh 3 and 5 by how many time periods they’ve been observed for (6 and 3, respectively). This would give us an average of 3.67. We could also find an average treatment by time-since-treatment like in an event study.

Let’s think about what a regression would give us. Suppose we have a TWFE model of the following form³:

\[Y_{it} = \alpha_i + \lambda_t + \delta^{DD}D_{it} + \epsilon_{it}\]

Intuitively, OLS should return an average of all the treatment effects in the data as \(\delta^{DD}\). This would be consistent with a time-weighted treatment effect of 3.67 in static case (or 6.67 in the dynamic case).

Let’s compare our time-weighted treatment effect with what TWFE returns.

Spoiler: our coefficient is (seemingly) too large in the static case and too small in the dynamic case.⁴

Goodman-Bacon

Wait, so what gives? Well, it’s complicated and we (they: the econometricians) just figured this all out.

As simple as this example is, this was a non-trivial observation to econometricians. Andrew Goodman-Bacon, the most recent scapegoat for grumbling economists, demonstrated in an immediately seminal work that the DiD coefficient in a TWFE framework (\(\delta^{DD}\)) is a (sometimes strangely) weighted average of all possible 2x2 DiD comparisons.

What does “all possible” 2x2 DiD comparisons mean?

In this three-group case with one untreated group, there are four possible 2x2 DiD comparisons.

1. The first two comparisons are the two treated groups each separately compared to the never-treated group.
2. The third comparison is the earlier treated group compared to the later treated group while the later treated group remains untreated.
3. The final comparison is the later treated group compared to the earlier treated group while the earlier treated group is always treated.

It’s easier to visualize this as the following:

df <- create_data(TRUE)
par(mfrow = c(2, 2), mar = c(2.1, 2.1, 2.1, 0.6))
## Comp. 1 -- Early vs Never
lim <- df$id %in% c(1, 2)
r1 <- fixest::feols(y ~ D | id + tt, df[lim,])
plot(df$tt[lim], df$y[lim], pch = 19,
     ylim = range(df$y),
     xlim = range(df$tt),
     xlab = "", ylab = "", las = 1,
     main = "I. Early vs Never",
     col = df$col[lim])
lines(df$tt[df$id == 1], df$y[df$id == 1], col = df$col[df$id == 1], lty = 2)
lines(df$tt[df$id == 2], df$y[df$id == 2], col = df$col[df$id == 2])
abline(v = c(4.5, 7.5), lty = 2)
legend("topleft", legend = c("Treatment", "Control"),
       lty = 1:2, bty = "n")

## Comp. 2 -- Late vs Never
lim <- df$id %in% c(1, 3)
r2 <- fixest::feols(y ~ D | id + tt, df[lim,])
plot(df$tt[lim], df$y[lim], pch = 19,
     ylim = range(df$y),
     xlim = range(df$tt),
     xlab = "", ylab = "", las = 1,
     main = "II. Late vs Never",
     col = df$col[lim])
lines(df$tt[df$id == 1], df$y[df$id == 1], col = df$col[df$id == 1], lty = 2)
lines(df$tt[df$id == 3], df$y[df$id == 3], col = df$col[df$id == 3])
abline(v = c(4.5, 7.5), lty = 2)

## Comp. 3 -- Early vs Late
lim <- df$id %in% c(2, 3) & df$tt < min(df$tt[df$id == 3 & df$D == 1])
r3 <- fixest::feols(y ~ D | id + tt, df[lim,])
plot(df$tt[lim], df$y[lim], pch = 19,
     ylim = range(df$y),
     xlim = range(df$tt),
     xlab = "", ylab = "", las = 1,
     main = "III. Early vs Late",
     col = df$col[lim])
lim2 <- lim & df$id == 2
lines(df$tt[lim2], df$y[lim2], col = df$col[lim2])
lim2 <- lim & df$id == 3
lines(df$tt[lim2], df$y[lim2], col = df$col[lim2], lty = 2)
abline(v = c(4.5, 7.5), lty = 2)


## Comp. 4 -- Late vs Early
lim <- df$id %in% c(2, 3) & df$tt >= min(df$tt[df$id == 2 & df$D == 1])
r4 <- fixest::feols(y ~ D | id + tt, df[lim,])
plot(df$tt[lim], df$y[lim], pch = 19,
     ylim = range(df$y),
     xlim = range(df$tt),
     xlab = "", ylab = "", las = 1,
     main = "IV. Late vs Early",
     col = df$col[lim])
lim2 <- lim & df$id == 2
lines(df$tt[lim2], df$y[lim2], col = df$col[lim2], lty = 2)
lim2 <- lim & df$id == 3
lines(df$tt[lim2], df$y[lim2], col = df$col[lim2])
abline(v = c(4.5, 7.5), lty = 2)

par(mfrow = c(1, 1))

Display Output

Each of these 2x2 comparisons will be given positive weights when making up \(\delta^{DD}\). However, we only necessarily want to consider comparisons I-III. Comparison IV is flawed because the “control” group is already tainted by treatment dynamics.

One of the innovations in Goodman-Bacon was the ability to decompose \(\delta^{DD}\) into it’s parts. Using the bacondecomp package, we can apply Goodman-Bacon’s diagnostic methodology to split up \(\delta^{DD}\) into its different 2x2 comparisons and their resulting estimates and weights.

library("bacondecomp")
bacon(formula = y ~ D, data = df,
      id_var = "id", time_var = "tt",
      quietly = FALSE) -> our_decomp

Standard Output

                      type  weight  avg_est
1 Earlier vs Later Treated 0.18182  3.00000
2 Later vs Earlier Treated 0.13636 -4.00000
3     Treated vs Untreated 0.68182  6.33333

This output summarizes the three possible comparisons: treated vs never-treated (good!), early vs late treated (good!), and late vs early treated (maybe bad?). What should raise an eyebrow, though, is the late vs early comparison group with an average estimate of -4! This is, of course, not representative of what we want OLS to reflect. Looking at the plot of our dynamic simulation data, it’s crystal clear that our treatment effects are positive.⁵

An aside about bacondecomp::bacon()

If we print the result of our decomposition, we can observe each 2x2. We can also always use this to re-aggregate to the summary bacon() output for us before.

our_decomp
cat("\n########## Summary ##########\n")
aggregate(list(weight = our_decomp$weight,
               avg_est = our_decomp$estimate * our_decomp$weight),
          list(type = our_decomp$type), sum) -> agg
agg$avg_est <- agg$avg_est / agg$weight
agg

Standard Output

  treated untreated estimate    weight                     type
2       5     99999      7.5 0.3636364     Treated vs Untreated
3       8     99999      5.0 0.3181818     Treated vs Untreated
6       8         5     -4.0 0.1363636 Later vs Earlier Treated
8       5         8      3.0 0.1818182 Earlier vs Later Treated

########## Summary ##########
                      type    weight   avg_est
1 Earlier vs Later Treated 0.1818182  3.000000
2 Later vs Earlier Treated 0.1363636 -4.000000
3     Treated vs Untreated 0.6818182  6.333333

This is no mistake. When I plotted the four comparisons, I also calculated the TWFE regressions for each. Displaying these coefficients, we can see how Goodman-Bacon’s decomposition gets it’s estimates.⁶

library("modelsummary")
options(modelsummary_factory_default = 'tinytable')
options("modelsummary_format_numeric_latex" = "plain")

modelsummary(list("E vs U" = r1, "L vs U" = r2,
                  "E vs L" = r3, "L vs E" = r4),
             gof_map = "", estimate = "{estimate}",
             title = "All 2x2 DiD Coefficients")

Display Output

All 2x2 DiD Coefficients

	E vs U	L vs U	E vs L	L vs E
D	7.500	5.000	3.000	-4.000
	(2.864)	(1.725)	(1.449)	(2.082)

With this very simple simulation, we’ve found a case where all treatment effects are positive (visually), but part of the decomposed \(\delta^{DD}\) is negative.⁷ In fact, it’s possible for this comparison to turn the entire treatment effect negative. See below for an example.⁸

baker <- haven::read_dta('https://github.com/scunning1975/mixtape/raw/master/baker.dta')
baker$treated <- ifelse(baker$treat_date > 0, 1, 0)

tmp <- aggregate(list(y = baker$y),
                 list(year = baker$year,
                      group = baker$group),
                 mean)
plot(tmp$year, tmp$y, pch = 19,
     col = scales::alpha(tmp$group, .4),
     xlab = "Year", ylab = "Outcome")
lines(tmp$year[tmp$group == unique(tmp$group)[1]],
      tmp$y[tmp$group == unique(tmp$group)[1]],
      col = tmp$group[tmp$group == unique(tmp$group)[1]])
lines(tmp$year[tmp$group == unique(tmp$group)[2]],
      tmp$y[tmp$group == unique(tmp$group)[2]],
      col = tmp$group[tmp$group == unique(tmp$group)[2]])
lines(tmp$year[tmp$group == unique(tmp$group)[3]],
      tmp$y[tmp$group == unique(tmp$group)[3]],
      col = tmp$group[tmp$group == unique(tmp$group)[3]])
lines(tmp$year[tmp$group == unique(tmp$group)[4]],
      tmp$y[tmp$group == unique(tmp$group)[4]],
      col = tmp$group[tmp$group == unique(tmp$group)[4]])
abline(v = unique(baker$treat_date)-0.5, col = 1:4)
cat("\n##### Goodman-Bacon Decomp #####\n")
bacon(formula = y ~ treat, data = baker,
      id_var = "id", time_var = "year",
      quietly = FALSE) -> this_decomp
cat("\n##### TWFE Results #####\n")
fixest::feols(y ~ treat | id + year, data = baker)

Standard Output

##### Goodman-Bacon Decomp #####
                      type weight   avg_est
1 Earlier vs Later Treated    0.5  51.80009
2 Later vs Earlier Treated    0.5 -65.18046

##### TWFE Results #####
OLS estimation, Dep. Var.: y
Observations: 30,000
Fixed-effects: id: 1,000,  year: 30
Standard-errors: IID 
      Estimate Std. Error  t value  Pr(>|t|)    
treat -6.69019   0.623867 -10.7237 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 29.1     Adj. R2: 0.803436
             Within R2: 0.003954

Display Output

As a note: unfortunately, the weights are not computed in any intuitive way. According to Andrew Baker⁹, “we can think of the TWFE estimator as a weighted average of individual 2x2 DiD estimators, with the weights proportional to group sizes and the variance of the treatment indicators in each pair, which is highest for units treated in the middle of the panel.” Therefore, \(\delta^{DD}\) is sensitive to both panel length and which groups appear towards the center. These are not necessarily things we want to influence our estimates. If you remain interested in how weights are derived, I encourage you to read Goodman-Bacon (2019).

Solutions

So, now what?

There have been many estimators proposed that solve this problem. In a blog post, Scott Cunningham notes three categories of these estimators: weighted group-time ATT, stacking in relative time, and imputation. While this is a helpful dichotomy to have, these estimators are all quite similar in reality. Of course, each has its own strengths and weaknesses, so it’s important to be familiar with more than just one, but the main thread through all of them is to fix the issues associated with wonky comparisons.

Some of the proposed estimators, just to name a few, are as follows:

• Gardner (2022)
• Dube et al. (2023)
• Callaway and Sant’Anna (2021)
• Sun and Abraham (2021)
• de Chaisemartin and D’Haultfouille (2020)
• Borusyak et al. (2023)
• Wooldridge (2021, 2022)

We do not have enough time to cover each method, so instead we are going to focus on two: Gardner (2022) and Sun and Abraham (2021). Gardner (2022) proposes something called “Two-Stage Difference-in-Differences”, while Sun and Abraham (2021) aggregate group-specific effects into a single quantity. While I would imagine that Callaway & Sant’Anna (2021) is the most “popular” estimator, I believe Gardner (2022) and Sun and Abraham (2021) are two of the more intuitive.

Gardner (2022)

The first solution we will discuss is, in my opinion, the most intuitive. This is one of those papers that make you think, “I wish I thought of that!” Gardner (2022) takes a clever approach towards estimating the TWFE model. The idea starts like this:

If the parallel trends assumption holds (is true), then the potential outcome for \(Y_{it}^0\) is simply a sum of unit \(i\)’s unit-specific effect and time period \(t\)’s common time-specific effect.¹⁰ Perhaps we should re-frame the problem and switch our priority to accurately measure these parameters, since they provide information about the unobservable counterfactual.

Gardner goes through quite a bit of math to effectively (albeit more elegantly) say the following: since we only have a single, static \(\delta^{DD}\) in our equation, any heterogeneity in treatment effects gets absorbed by the fixed effects \(\alpha_i\) and \(\lambda_t\). Informally, OLS is trying to fit a square peg in a round hole. Forcing the peg and biasedly estimating \(\alpha_i\) and \(\lambda_t\) in turn biases \(\widehat{\delta^{DD}}\)!

The solution Gardner proposes is rather simple: if estimating \(\alpha_i\) and \(\lambda_t\) causes us problems, let’s remove them. Wait… how can we just remove them?! If we can somehow get unbiased estimates of our fixed effects, we can simply subtract them from our outcome before estimating \(\delta^{DD}\)!

Before moving on, let’s think about what fixed effects are. Suppose the data we have contain Virginia from 2010 to 2019. A Virginia-fixed-effect is a time-invariant effect. So, if we estimate a Virginia-fixed-effect, it shouldn’t matter whether we use data from 2010-2014, 2015-2019, strictly odd years, or any other combination of years. In theory, all estimations should yield the same (or similar) result. The same logic holds for time-period-fixed-effects: it shouldn’t matter which units you use to estimate them since the effect is invariant across units.

Since we know treatment is messing up our fixed effects, we can estimate \(\alpha_i\) and \(\lambda_t\) using only data from untreated time periods. We will call this our first-stage. The resulting fixed effects from this estimating will necessarily be unbiased (assuming parallel trends) since there is no treatment involved. Then, in a second stage, we can fit the following model using the entire dataset:

\[ (Y_{it} - \widehat{\alpha}_i - \widehat{\lambda}_t) = \delta^{DD}D_{it} + \epsilon_{it} \]

Let’s take a look at the mechanics of this. We are going to start with the first stage, and regress our outcome variable on only fixed effects only using observations from untreated units.

r1 <- lm(y ~ as.factor(id) + as.factor(tt) - 1, df[df$D == 0,])
coef(r1)

Standard Output

 as.factor(id)1  as.factor(id)2  as.factor(id)3  as.factor(tt)2  as.factor(tt)3 
              2               3               4               1               2 
 as.factor(tt)4  as.factor(tt)5  as.factor(tt)6  as.factor(tt)7  as.factor(tt)8 
              3               4               5               6               7 
 as.factor(tt)9 as.factor(tt)10 
              8               9 

Next, we are going to collect the estimated fixed effects and subtract them from the outcome variable.

fe_id <- coef(r1)[grepl("id", names(coef(r1)))]
fe_tt <- coef(r1)[grepl("tt", names(coef(r1)))]
df$fe_id <- fe_id[match(df$id, gsub("as.factor\\(id|\\)", "", names(fe_id)))]
df$fe_tt <- fe_tt[match(df$tt, gsub("as.factor\\(tt|\\)", "", names(fe_tt)))]
df$fe_tt <- ifelse(is.na(df$fe_tt), 0, df$fe_tt)
df$y_2sdid <- df$y - df$fe_id - df$fe_tt

Finally, estimate the second-stage regression without fixed effects.

lm(y_2sdid ~ D, df)

Standard Output

Call:
lm(formula = y_2sdid ~ D, data = df)

Coefficients:
(Intercept)            D  
 -3.243e-16    6.667e+00  

Notice how this is exactly equivalent to our average treatment effect of 6.67 for the dynamic group! The only issue with doing this two-step procedure is that the standard errors can get weird. This is true for all two-stage procedures. There’s an R package called did2s that can help implement two-stage DiD for you. I highly recommend this package.

library("did2s")
did2s(df, yname = "y",
      first_stage = ~ 0 | id + tt,
      second_stage = ~ D, treatment = "D",
      cluster_var = "id")

Standard Output

OLS estimation, Dep. Var.: y
Observations: 30
Standard-errors: Corrected Clustered (id) 
  Estimate Std. Error t value  Pr(>|t|)    
D  6.66667   0.785674 8.48528 2.383e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 2.70801   Adj. R2: 0.544828

For more information about this package, see the vignette and/or the R Journal article.

In addition to estimating \(\delta^{DD}\), we can also estimate an event study using 2SDiD. To do this, simply replace the treatment variable with the fixest notation for an event study. See Module 3.1 for this information.

As a final note, keep in mind that in order to estimate the unit and/or time fixed effect, you need to observe some untreated units in both. What I mean is that there cannot be any units that are always-treated and there cannot be any times that only contain treated units. Their fixed effects will, of course, not be estimated in the first stage, so they automatically get removed from the second stage.

Sun and Abraham (2021)

The standard TWFE event study framework is as follows:

\[Y_{it} = \alpha_i + \lambda_t + \sum_{k\neq-1}\delta_k D_{it}^{k} + u_{it}\]

where \(D_{it}^{k}\) is a dummy equal to 1 if unit \(i\) is \(k\) periods from treatment at time \(t\). [When treatment effects are heterogeneous across cohorts, that is, units treated in different periods have different dynamic paths, \(\widehat{\delta}_k\) is not a clean average of cohort-specific effects at horizon \(k\). It gets contaminated by treatment effects from other relative periods, because already-treated cohorts end up serving as part of the comparison group for later cohorts through the two-way fixed effects.

As mentioned, there are three main ways to classify the “solutions” to the problems introduced by staggered treatment timing. Gardner’s fix was a two-stage procedure, imputing untreated potential outcomes from a first-stage model, then comparing observed to imputed outcomes in a second stage. Sun and Abraham (2021) takes a different route entirely, and stays inside a single regression.

The problem is that \(D_{it}^{k}\) pools observations across cohorts \(g\) as if one \(\delta_k\) could describe every cohort’s effect at horizon \(k\). Sun and Abraham (2021)’s fix is to stop pooling, let each cohort have its own relative-time path, and then average those paths after the fact.

Instead of a single dummy per relative period, Sun and Abraham (2021) interacts the relative-time dummies with cohort indicators \(\mathbb{1}\{G_i = g\}\), where \(G_{i}\) is the period unit \(i\) is first treated:

\[Y_{it} = \alpha_i + \lambda_t + \sum_{g}\sum_{k\neq-1}\delta_{g,k} (\mathbb{1}\{G_i = g\} \times D_{it}^{k}) + u_{it}\]

Now, \(\delta_{g,k}\) is cohort \(g\)’s treatment effect at relative time \(k\), estimated only off variation within that cohort’s own path relative to never-treated (or last-treated) units. In this set up, there are no forbidden comparisons, since the comparison group is fixed and explicit rather than baked into a shared TWFE.

The issue with \(\delta_{g,k}\), though, is that it is specific to each cohort, whereas what we’d really like is a single estimate or event-study (like 2S-DiD). To get our target parameters, we can aggregate the \(\delta_{g,k}\) parameters using cohort shares as weights. Luckily, fixest can do this for us with sunab().

# I need to merge two datasets because
#   one obs per group-time is breaking the SEs. 
df0 <- create_data(TRUE, TRUE)
df1 <- create_data(TRUE, TRUE)
df <- rbind(df0, df1)

# Construct group indicators by treatment time.
agg <- aggregate(df$tt[df$D == 1], list(df$id[df$D == 1]), min)
df$g <- agg$x[match(df$id, agg$Group.1)]
df$g <- ifelse(is.na(df$g), Inf, df$g)
# Untreated units are given 'Inf'

mod_sa <- feols(y ~ sunab(g, tt) | id + tt, data = df)
# or: feols(y ~ sunab(g, tt, att = T) | id + tt, data = df)
summary(mod_sa, agg = "ATT")

Standard Output

OLS estimation, Dep. Var.: y
Observations: 60
Fixed-effects: id: 3,  tt: 10
Standard-errors: IID 
    Estimate Std. Error      t value  Pr(>|t|)    
ATT   6.6636   6.55e-15 1.017208e+15 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 5.673e-15     Adj. R2: 1
                  Within R2: 1

To get the event-study version, we can just use iplot().

Stacked Difference-in-Differences

The next approach we’ll discuss is stacked difference-in-differences. There is not a specifically associated paper with this estimator, but the intuition has been used in multiple studies. One of the most famous studies that uses a variation of this method is Cengiz et at. (2019), but Wing et al. (2024) is another good resource.

First, let’s define some terms:

1. Event: an event is going to be any time some policy is changed or some shock is experienced within a unit (e.g. state, etc.) For example, New Jersey changing its minimum wage would constitute an event.
2. Sub-experiment: a sub-experiment is a collection of treatment and clean control groups. We’ll define “clean” in a moment. An event might be New Jersey changing its minimum wage. However, if New Jersey and Michigan change their minimum wages at the same time, then their events, in addition to their controls, would constitute a single sub-experiment. The way to refer to New Jersey’s and Michigan’s events would be a cohort.
3. Clean control: a clean control group is made up of units that do not experience any treatment, or residual treatment dynamics, during a sub-experiment’s treatment periods.

Let’s start with a fresh df so we can experiment with stacking.

df <- create_data(TRUE)
# Since this dataset is so small, I need to increase
#   the number of observations for this demo.
#   You would **not** do this in real life.  
df <- rbind(df, df)
df <- df[order(df$id, df$tt),]

Since there are two treated units in this example, we would have two sub-experiments (one for each event). To focus in on the treatment, we can create a window of \(\pm2\) time period relative to treatment for each sub-experiment.

par(mfrow = c(1, 2), mar = c(2.1, 2.1, 2.1, 0.6))
lim <- df$tt %in% (5 + (-2:2))
df1 <- df[lim,]
plot(df$tt[lim], df$y[lim],
     main = "Sub-Experiment 1",
     col = df$col[lim], pch = 19,
     xlim = range(df$tt), ylim = range(df$y),
     xlab = "Time", ylab = "Y", las = 1)
for(i in 1:3){
    
  lines(df$tt[df$id == i & lim],
        df$y[df$id == i & lim],
        col = df$col[df$id == i & lim])
}
abline(v = c(4.5, 7.5), lty = 2)

lim <- df$tt %in% (8 + (-2:2)) & df$id != 2
df2 <- df[lim,]
plot(df$tt[lim], df$y[lim],
     main = "Sub-Experiment 2",
     col = df$col[lim], pch = 19,
     xlim = range(df$tt), ylim = range(df$y),
     xlab = "Time", ylab = "Y", las = 1)
for(i in 1:3){
    
  lines(df$tt[df$id == i & lim],
        df$y[df$id == i & lim],
        col = df$col[df$id == i & lim])
}
abline(v = c(4.5, 7.5), lty = 2)
par(mfrow = c(1, 1))

Display Output

Now that we have our two separately defined sub-experiments, with clean control units, we can create variables to enumerate each experiment and denote relative time within each experiment. Finally, we’ll stack the two subsets of data on top of each other in relative time. By “stack”, I mean that we will literally be vertically concatenating the data. In R, this is done via rbind().

# Enumerate the sub-experiment
df1$experiment <- 1
df2$experiment <- 2
# Add in relative time
df1$rel_time <- df1$tt - 5
df2$rel_time <- df2$tt - 8
# Stack the data with rbind()
stacks <- rbind(df1, df2)
# Creat a treatment variable
stacks$ever_treat <- ave(stacks$D, paste0(stacks$id, stacks$experiment), FUN = max)

Using these stacked data, we can now estimate a simple 2x2 TWFE specification. Note: the fixed effects should be interacted with the sub-experiment identifier. With fixest, we can do this by using FE1^subexperiment.

feols(y ~ D | id^experiment + tt^experiment, stacks)

Standard Output

OLS estimation, Dep. Var.: y
Observations: 50
Fixed-effects: id^experiment: 5,  tt^experiment: 10
Standard-errors: IID 
  Estimate Std. Error t value  Pr(>|t|)    
D  3.85714   0.884924 4.35873 0.0001093 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 1.23889     Adj. R2: 0.870041
                Within R2: 0.351834

As I see it, there are two main advantages of stacking. First, after stacking, we’re back to the familiar 2x2 case that we’re already comfortable with. Moreover, and this is the second benefit, we can also estimate \(\delta^{DD}\) for each sub-experiment and summarize/combine the results however we’d like. See below for an example of how to do this with fixest.

# Estimate an event study by sub-experiment using split
feols(y ~ i(rel_time, ever_treat, -1)
      | id^experiment + tt^experiment,
      stacks,
      split = ~experiment) -> r2
iplot(r2, grid = F)

# Estimate the aggregate event study
feols(y ~ i(rel_time, ever_treat, -1)
      | id^experiment + tt^experiment, stacks) -> r1
iplot(r1, add = TRUE, col = "dodgerblue")

Display Output

While this is seems easy enough with just three groups, the complexity increases when the number of sub-experiments and/or the size of the data becomes large. This becomes even more difficult when treatment is non-absorbing. Treatments that are absorbing forever alter the trajectory of the unit. It is the opposite for non-absorbing treatments, which “wears off”. In other words, after some time, the effect of the treatment levels off, a new equilibrium is reached, and the unit effectively returns back to being untreated.

The type of treatment, absorbing or non-absorbing, does not imply the magnitude of the effect, but rather whether the treated units can eventually return to the control group after some time. This is an assumption that necessarily depends on the institutional setting.

An example of non-absorbing treatment would be a pipe in my apartment bursting. I would likely have to move out for a few days or so while the pipe is fixed, but eventually I would be able to return to my normal life. An example of absorbing treatment would be my landlord evicting me from my apartment. I would need to find a new place to live and (not to be dramatic) my life would be changed forever (in a causal inference sense).

More concrete examples would be minimum wage laws compared to texting-and-driving laws. Minimum wages are increased every few years to account for inflation, which necessarily makes the changes non-absorbing. On the other hand, once texting-and-driving has been made illegal it will always be illegal.

Unfortunately, there are not any (well developed) R packages that address stacked difference-in-differences. However, that last sentence may be out-dated by the time you read it with how quickly theory and software are developing around this topic. Again, at the time of writing, I am working on developing a package to do this.

Method Comparison

Finally, let’s revisit that dataset where all the treatment effects were clearly positive but the TWFE estimate turned out negative. We can compare the results the 2SDiD, SA and TWFE. Let’s try all of these with a pooled effect:

baker <- as.data.frame(baker)
r_twfe <- feols(y ~ treat | id + year, data = baker,
                cluster = ~ state)
r_did2s <- did2s(baker, "y",
                 first_stage = ~ 0 | id + year,
                 second_stage = ~ treat,
                 treatment = "treat",
                 cluster_var = "state")
r_sunab <- feols(y ~ sunab(treat_date, year) | id + year, 
                 baker[baker$year<2004,],
                 cluster = ~state)

cat("TWFE DiD:\n")
r_twfe$coeftable; cat("\n\nTwo-Stage DiD:\n")
r_did2s$coeftable; cat("\n\n\nSun & Abraham:\n")
summary(r_sunab, agg = "ATT")$coeftable

Standard Output

TWFE DiD:
       Estimate Std. Error   t value   Pr(>|t|)
treat -6.690186   3.391857 -1.972426 0.05568123
attr(,"vcov_type")
[1] "Clustered (state)"


Two-Stage DiD:
      Estimate Std. Error  t value   Pr(>|t|)
treat 68.32906   5.191276 13.16229 1.9815e-39
attr(,"vcov_type")
[1] "Corrected Clustered (state)"



Sun & Abraham:
    Estimate Std. Error  t value      Pr(>|t|)
ATT 68.33072 0.01417939 4819.018 3.118791e-114
attr(,"vcov_type")
[1] "Clustered (state)"

Last but not least, let’s visualize the event study versions.

es_twfe <- feols(y ~ i(time_til, treated, -1) | id + year, data = baker)
es_did2s <- did2s(baker, "y",
                  first_stage = ~ 0 | id + year,
                  second_stage = ~ i(time_til, treated, -1),
                  treatment = "treat",
                  cluster_var = "id")
iplot(es_twfe, grid = F, xlab = "Time until Treatment", main = "")
iplot(es_did2s, add = T, x.shift = 0.1, col = "tomato")
iplot(r_sunab, add = T, x.shift = -0.1, col = "dodgerblue")
legend("topright", c("TWFE", "2SDiD", "SA"),
       col = c("black", "tomato", "dodgerblue"),
       bty = "n", pch = 19, ncol = 2)

Display Output

Before wrapping up, it is important to make the following note. Figures like the above are bad and you should not be doing this. Although each estimator is similar, there are subtle differences. While figures like the above are tempting, it is more likely to signal a lack of understanding than robustness. See Roth (2026).

Footnotes

1. The beautiful, sentient being that it is

2. That said, if you want to change how create_data works, you can copy, paste, and modify it in one of the following webr code chunks.

3. Notation taken from Baker et al. (2022)

4. Note: it’s not the case that static \(\rightarrow\) too small and dynamic \(\rightarrow\) too large. This is just an artifact of our simulated data. My point is that the coefficient estimate can be either too large or too small and it’s not necessarily clear which will happen.

5. Note: in real life, things are hardly ever as clear!

6. Note: “U” represents never-treated units, “E” represents early-treated units, and “L” represents later-treated units.

7. This is not to say that some parts cannot be negative. In real life, this will likely be the case. However, in this illustration, it’s clear that the negative is unintended.

8. Note: in these data, there is no never-treated group. Therefore, the decomposition only contains two groups: early vs late and late vs early.

9. This is an excellent blog post, and would be worth a read.

10. I am using \(Y^0\) to denote \(Y\) as untreated.