Background

This vignette is written under the assumption that you’re at least vaguely familiar with developments in difference-in-differences with staggered adoptions since about 2018. Just to make sure we’re on the same page, the brief recap is:

• Historically, under staggered adoptions researchers used the standard two-way fixed effects estimator and interpreted the coefficient on the treatment dummy as an average treatment effect on the treated units.
• In the late 2010s, econometricians formally checked what this estimator was doing and found that in fact, this coefficient was not any kind of reasonable average treatment effect estimator.
• Since then, a number of new difference-in-differences estimators that are unbiased under staggered adoptions have been developed.

The estimator in this package, fused extended two-way fixed effects (FETWFE), is one of those unbiased estimators. Of course, I made this estimator because I think FETWFE brings something to the table that the others don’t. Here’s a brief summary on that:

One issue with these estimators has been that they’ve worked so hard to be unbiased that they are inefficient (in the language of econometrics), or high-variance (in the language of machine learning). These estimators add extra parameters in order to remove bias, but estimating extra parameters means you have less data per parameter and your estimates are noisier.

In machine learning, creating a more flexible estimator with lots of parameters and then finding that it is too high variance (that is, it overfits) is a familiar issue. The most common solution has been regularization.

You could just add \(\ell_2\) or \(\ell_1\) regularization to a difference-in-differences regression estimator and probably see an improvement in your efficiency, but FETWFE does something more sophisticated than that. (Plus, that approach wouldn’t allow you to get valid standard errors for your treatment effect estimates, but FETWFE does.)

That’s all the description I’ll give you in this vignette. You can learn all of the details in the paper on arXiv:

Fused Extended Two-Way Fixed Effects for Difference-in-Differences With Staggered Adoptions

If you want to learn a little more before you dive into the full paper, here are some other resources with descriptions of the methodology that provide a little more detail than this vignette:

• My blog post announcing the paper.
• Some slides I made for a presentation on FETWFE.
• Another blog post focused on what this package is doing under the hood.

But the headline summary of what fused extended two-way fixed effects brings to the table in a crowded field of estimators is: fused extended two-way fixed effects is not only unbiased, it also uses machine learning to maximize efficiency (minimize variance). Further, unlike many machine learning estimators, fused extended two-way fixed effects gives you valid standard errors for the treatment effect estimates.

Simulated Data Example

I’ll start illustrating how to use fetwfe() by using a simulated data set. The example below simulates a balanced panel with 60 time periods, 30 individuals, and 5 waves of treatment.

In the simulation, each individual is assigned a random cohort (which determines the timing of treatment) and three time-invariant covariates are generated. The response variable is constructed so that, after treatment, its evolution depends on a treatment effect (which varies by cohort) and a linear trend, plus the covariates and some random noise.

Below is the complete code for simulating the data, converting it into the required pdata format, and running the fetwfe() function.

I borrowed some of the below code from Asjad Naqvi’s helpful website for DiD estimators. Thanks for sharing the code publicly!

# Set seed for reproducibility
set.seed(123456L)

# 60 time periods, 30 individuals, and 5 waves of treatment
tmax = 60; imax = 30; nlvls = 5

dat = 
  expand.grid(time = 1:tmax, id = 1:imax) |>
  within({
    
    # Generate time-invariant covariates
    cov1 = rep(runif(imax, 0, 1), each = tmax)  # Random uniform values (0, 1) per individual
    cov2 = rep(sample(1:5, imax, replace = TRUE), each = tmax)  # Random categorical values (1-5) per individual
    cov3 = rep(rnorm(imax, mean = 0, sd = 1), each = tmax)  # Random Gaussian values (mean=0, sd=1) per individual
    
    # Initialize columns
    cohort      = NA
    effect      = NA
    first_treat = NA
    
    for (chrt in 1:imax) {
      cohort = ifelse(id==chrt, sample.int(nlvls, 1), cohort)
    }
    
    for (lvls in 1:nlvls) {
      effect      = ifelse(cohort==lvls, sample(2:10, 1), effect)
      first_treat = ifelse(cohort==lvls, sample(1:(tmax+20), 1), first_treat)
    }
    
    first_treat = ifelse(first_treat>tmax, Inf, first_treat)
    treat       = time >= first_treat
    rel_time    = time - first_treat
    y           = id + time + ifelse(treat, effect*rel_time, 0) + cov1 + cov2 + cov3 + rnorm(imax*tmax)
    
    rm(chrt, lvls, cohort, effect)
  })

head(dat)

##   time id        y rel_time treat first_treat      cov3 cov2      cov1
## 1    1  1 3.995163     -Inf FALSE         Inf 0.1582893    2 0.7977843
## 2    2  1 7.215283     -Inf FALSE         Inf 0.1582893    2 0.7977843
## 3    3  1 6.089467     -Inf FALSE         Inf 0.1582893    2 0.7977843
## 4    4  1 8.205118     -Inf FALSE         Inf 0.1582893    2 0.7977843
## 5    5  1 8.869905     -Inf FALSE         Inf 0.1582893    2 0.7977843
## 6    6  1 9.463189     -Inf FALSE         Inf 0.1582893    2 0.7977843

The simulated data (dat) now has columns for time, id, covariates (cov1, cov2, cov3), a treatment indicator (treat), and a response variable (y). Next, we convert this data into the panel data format required by fetwfe().

library(dplyr)

# Specify column names for the pdata format
time_var <- "time"       # Column for the time period
unit_var <- "unit"       # Column for the unit identifier
treatment <- "treated"   # Column for the treatment dummy indicator
covs <- c("cov1", "cov2", "cov3")  # Columns for covariates
response <- "response"   # Column for the response variable

# Convert the dataset
pdata <- dat |>
  mutate(
    # Rename id to unit and convert to character
    {{ unit_var }} := as.character(id),
    # Ensure treatment dummy is 0/1
    {{ treatment }} := as.integer(treat),
    # Rename y to response
    {{ response }} := y
  ) |>
  select(
    {{ time_var }}, {{ unit_var }}, {{ treatment }}, all_of(covs), {{ response }}
  ) 

# Preview the resulting pdata dataframe
head(pdata)

##   time unit treated      cov1 cov2      cov3 response
## 1    1    1       0 0.7977843    2 0.1582893 3.995163
## 2    2    1       0 0.7977843    2 0.1582893 7.215283
## 3    3    1       0 0.7977843    2 0.1582893 6.089467
## 4    4    1       0 0.7977843    2 0.1582893 8.205118
## 5    5    1       0 0.7977843    2 0.1582893 8.869905
## 6    6    1       0 0.7977843    2 0.1582893 9.463189

Now that pdata is properly formatted, we run the FETWFE estimator on the simulated data.

library(fetwfe)

# Run the FETWFE estimator on the simulated data
result <- fetwfe(
  pdata = pdata,              # The panel dataset
  time_var = "time",          # The time variable
  unit_var = "unit",          # The unit identifier
  treatment = "treated",      # The treatment dummy indicator
  covs = c("cov1", "cov2", "cov3"),  # Covariates
  response = "response"      # The response variable
)

# Display the overall average treatment effect estimate
cat("Estimated Overall ATT:", result$att_hat, "\n")

## Estimated Overall ATT: 32.82035

When you run this code, the function internally performs all the necessary data preparation, applies the fusion penalty via a bridge regression (using the grpreg package), and returns a list with overall and cohort-specific treatment effect estimates, standard errors (if available), and additional diagnostics.

A “Real Data” Example

Next I illustrate FETWFE in an empirical context. I’ll use data from Stevenson and Wolfers (2006), via the divorce data set from the bacondecomp package, on the impact of no-fault divorce laws on women’s suicide rates. (See also Goodman-Bacon (2021) for an alternative analysis.) The below is identical to the data application from my paper.

In this application, the panel data consist of state-level observations over 33 years. After removing states that received treatment in the first period, we are left with 42 states, of which 5 are never treated and 12 cohorts adopt treatment at various times.

Time-varying covariates (such as the state homicide rate, logged personal income, and welfare participation) are included as controls (using the pre-treatment values). In this example, I use FETWFE to estimate the marginal average treatment effect of no-fault divorce laws. Note: the below code takes about 15 seconds to run on my laptop.

library(bacondecomp)  # for the example data

# Load the example data
data(divorce)

set.seed(23451)

# Suppose we wish to estimate the effect of a policy (here represented by the variable "changed")
# on the response "suiciderate_elast_jag" using covariates "murderrate", "lnpersinc", and "afdcrolls".
# Here
# - 'year' is the time period variable (as an integer),
# - 'st' is the unit identifier,
# - 'changed' is the treatment indicator (with 0 = untreated, 1 = treated),
# 
# The `fetwfe()` function will automatically take care of removing units that were treated in the
# first time period.

# Call the estimator
res <- fetwfe(
    pdata=divorce[divorce$sex == 2, ],
    time_var="year",
    unit_var="st",
    treatment="changed",
    covs=c("murderrate", "lnpersinc", "afdcrolls"),
    response="suiciderate_elast_jag"
    )

## Warning in idCohorts(df = data, time_var = time_var, unit_var = unit_var, : 9
## units were removed because they were treated in the first time period

## Warning in processCovs(df = data, units = units, unit_var = unit_var, times =
## times, : 1 covariates were removed because they contained missing values.

# Average treatment effect on the treated units (in percentage point
# units)
100 * res$att_hat

## [1] -3.76012

# Conservative 95% confidence interval for ATT (in percentage point units)

low_att <- 100 * (res$att_hat - qnorm(1 - 0.05 / 2) * res$att_se)
high_att <- 100 * (res$att_hat + qnorm(1 - 0.05 / 2) * res$att_se)

c(low_att, high_att)

## [1] -12.73424   5.21400

# Cohort average treatment effects and confidence intervals (in percentage
# point units)

catt_df_pct <- res$catt_df
catt_df_pct[["Estimated TE"]] <- 100 * catt_df_pct[["Estimated TE"]]
catt_df_pct[["SE"]] <- 100 * catt_df_pct[["SE"]]
catt_df_pct[["ConfIntLow"]] <- 100 * catt_df_pct[["ConfIntLow"]]
catt_df_pct[["ConfIntHigh"]] <- 100 * catt_df_pct[["ConfIntHigh"]]

catt_df_pct

##    Cohort Estimated TE       SE ConfIntLow ConfIntHigh
## 1    1969     0.000000 0.000000   0.000000   0.0000000
## 2    1970   -40.142218 7.813430 -55.456259 -24.8281776
## 3    1971     0.000000 0.000000   0.000000   0.0000000
## 4    1972     0.000000 0.000000   0.000000   0.0000000
## 5    1973    -3.465602 2.169996  -7.718716   0.7875132
## 6    1974     0.000000 0.000000   0.000000   0.0000000
## 7    1975     0.000000 0.000000   0.000000   0.0000000
## 8    1976    -4.702911 3.807125 -12.164740   2.7589171
## 9    1977    -5.338497 3.503595 -12.205417   1.5284229
## 10   1980     0.000000 0.000000   0.000000   0.0000000
## 11   1984     0.000000 0.000000   0.000000   0.0000000
## 12   1985     0.000000 0.000000   0.000000   0.0000000

For the data application, FETWFE yielded an overall ATT of approximately –3.76% change in the female suicide rate, similar to other estimates in the literature. In addition, the output table (stored in result_emp$catt_df) displays the cohort-specific estimates. (Note that standard errors for the individual cohort estimates are less reliable when the number of units per cohort is small.)