Skip to contents

Introduction

The fixes package provides an easy-to-use toolkit for estimating and visualizing difference-in-differences designs on panel data: event-study curves via event_study(), aggregated ATT via att(), and basic two-way fixed-effects DiD via did(). Every result has a base plot() method, so a complete analysis is a two-line pipeline.

Estimation runs on the package’s internal C++ fixed-effects OLS engine; the fixest package is an optional dependency (used only in this vignette for its example datasets and a few advanced specifications).

Upgrading from fixes < 1.0.0? The old verb-style functions (run_es(), calc_att(), run_did(), plot_es(), …) still work and return identical results, but are deprecated. This vignette uses the new noun-style API introduced in 1.0.0.

Installation

Install the released version from CRAN:

Or with pak (recommended for fast install):

pak::pak("fixes")

To install the latest development version from GitHub:

pak::pak("yo5uke/fixes")

Minimal Example

Below is a basic example using the built-in fixest::base_did dataset, running an event study, and visualizing the results.

# Load example data
df <- fixest::base_did

# Run the event study (supports multiple confidence levels)
es <- event_study(
  data       = df,
  outcome    = y,
  treatment  = treat,
  time       = period,
  timing     = 5,  # Treatment occurs at period 5 (universal timing)
  fe         = ~ id + period,
  cluster    = ~ id,
  baseline   = -1,
  interval   = 1,
  lead_range = 3,
  lag_range  = 3,
  conf_level = c(0.90, 0.95, 0.99)  # Multiple CIs supported!
)

# View results
head(es)
#> Event Study Result (fixes)
#>   N: 1080  | Units: NA  | Treated units: 1080  | Never-treated: NA 
#>   FE: id + period
#>   VCOV: cluster  | Cluster: id 
#>   Estimator: twfe 
#>   Method: classic  | lead_range: 3  lag_range: 3  baseline: -1

Because timing = 5 is a scalar, event_study() recognizes a universal (non-staggered) design automatically; when timing names a column of adoption times, a staggered design is inferred instead.

Visualizing Event Study Results

Results plot directly with plot(). You can switch between ribbon-style or error bar CIs, select the displayed CI level, and customize appearance.

# Basic plot (default: ribbon, 95% CI)
plot(es)

# Plot with error bars and 99% CI
plot(es, type = "errorbar", ci_level = 0.99)

# Customize further with ggplot2
plot(es, type = "errorbar", ci_level = 0.9, theme_style = "classic") +
  scale_x_continuous(breaks = seq(-3, 3, by = 1)) +
  ggtitle("Event Study with 90% CI and Classic Theme")

For interactive exploration with hover tooltips (estimates, CIs, standard errors, p-values), pass interactive = TRUE (requires the suggested plotly package):

plot(es, interactive = TRUE)

Package Highlights

  • event_study():

    • Fast, one-step event study for panel data on an internal C++ fixed-effects engine.
    • Six estimators behind one interface (estimator =): classic TWFE, Callaway & Sant’Anna, Sun & Abraham, Borusyak-Jaravel-Spiess, Wooldridge two-way Mundlak, and Deb et al. FLEX.
    • Automatic creation of lead/lag dummies relative to treatment, with staggered/universal timing inferred from timing.
    • Covariates, clustering, weights, and flexible baseline normalization.
    • Multiple confidence interval levels (e.g., 90%, 95%, 99%).
    • Handles irregular time panels via time_transform.
  • att(): aggregated ATT — overall, by cohort, or by calendar time — for the CS and BJS estimators, with broom::tidy()/glance() methods.

  • did(): single-coefficient TWFE DiD compatible with modelsummary::modelsummary().

  • plot(): every result class has a plot method — event-study curves, ATT point-range charts, ATT(g,t) heatmaps (plot(att_gt(x))), honest sensitivity plots, and contamination-weight heatmaps.

Staggered adoption estimators

Standard TWFE event-study regressions can produce biased and sign-reversed estimates under heterogeneous treatment effects when units adopt treatment at different times (Callaway & Sant’Anna 2021; Sun & Abraham 2021; Borusyak, Jaravel & Spiess 2024). fixes provides robust alternatives, all accessible through the same event_study() interface via the estimator argument.

We demonstrate three of them on fixest::base_stagg, a simulated staggered-adoption panel with true ATT = 1 for all cohorts and horizons.

df_stagg <- fixest::base_stagg
# Mark never-treated units with NA (convention for all staggered estimators)
df_stagg$timing <- df_stagg$year_treated
df_stagg$timing[df_stagg$year_treated == 10000] <- NA

Callaway & Sant’Anna (2021) — estimator = "cs"

Estimates a separate ATT(g,t) for every cohort-by-period cell, then aggregates to an event-study curve using cohort-size weights.

res_cs <- event_study(
  data          = df_stagg,
  outcome       = y,
  time          = year,
  timing        = timing,
  unit          = id,
  estimator     = "cs",
  control_group = "nevertreated"
)
plot(res_cs) + ggplot2::ggtitle("Callaway & Sant'Anna (2021)")

Sun & Abraham (2021) — estimator = "sa"

Builds cohort x relative-time interactions and aggregates with cohort-share weights. Numerically identical to fixest::sunab() (which the deprecated run_es(method = "sunab") used directly).

res_sa <- event_study(
  data      = df_stagg,
  outcome   = y,
  treatment = treated,
  time      = year,
  timing    = timing,
  unit      = id,
  fe        = ~ id + year,
  estimator = "sa",
  cluster   = ~ id
)
plot(res_sa) + ggplot2::ggtitle("Sun & Abraham (2021)")

Borusyak, Jaravel & Spiess (2024) — estimator = "bjs"

Fits TWFE on untreated observations, imputes counterfactuals for treated units, then averages by horizon.

res_bjs <- event_study(
  data      = df_stagg,
  outcome   = y,
  time      = year,
  timing    = timing,
  unit      = id,
  estimator = "bjs"
)
plot(res_bjs) + ggplot2::ggtitle("Borusyak, Jaravel & Spiess (2024)")

Under homogeneous treatment effects (as in this DGP), the estimators give similar results, each recovering a post-treatment ATT close to the true value of 1. Wooldridge’s two-way Mundlak ("twm") and the Deb et al. FLEX estimator for repeated cross-sections ("flex") round out the set.

Aggregated ATT — att()

When you want a single number (or one per cohort / per period) instead of the full curve:

att(df_stagg, outcome = y, time = year, timing = timing, unit = id,
    estimator = "cs", aggregation = "simple")
#> ATT Estimation  [estimator: CS | aggregation: Simple (overall)]
#> N = 950 obs | 95 units | 45 treated
#> 
#>   group estimate std.error statistic  p.value conf_low_95 conf_high_95
#> 1    NA   -0.755     0.226     -3.35 0.000813        -1.2       -0.313

Bootstrap simultaneous confidence bands

Pointwise vs. simultaneous CIs

Standard confidence intervals are pointwise: the 95% CI at each horizon covers the true value with 95% probability, but the probability that all intervals simultaneously cover their true values can be much lower when many periods are plotted.

Simultaneous confidence bands (Callaway & Sant’Anna 2021, Corollary 1) control the joint coverage probability. With probability at least 1 − α, the entire event-study curve is contained within the simultaneous band. This is especially important for parallel-trends pre-testing: a pre-trend test based on simultaneous bands controls the family-wise error rate.

Usage

Pass bootstrap = TRUE to event_study() together with the CS estimator. The multiplier bootstrap (Algorithm 1 of Callaway & Sant’Anna 2021) is used.

# NOTE: boot_reps = 199 shown here for brevity; use 999 in practice
res_cs_boot <- event_study(
  data          = df_stagg,
  outcome       = y,
  time          = year,
  timing        = timing,
  unit          = id,
  estimator     = "cs",
  control_group = "nevertreated",
  bootstrap     = TRUE,
  boot_reps     = 199,
  boot_seed     = 42
)
# The lighter outer band is the simultaneous CI; the darker inner band is
# the standard pointwise CI.
plot(res_cs_boot, show_simultaneous = TRUE)

The simultaneous critical value ĉ and per-period simultaneous CI bounds are stored in attr(res_cs_boot, "bootstrap"). The simultaneous band is always at least as wide as the pointwise band, since the critical value ĉ_{1-α} >= z_{1-α/2}.


ATT(g,t) visualization

The CS estimator produces a separate ATT estimate for every (cohort g, calendar period t) pair. Extract the full matrix with att_gt() and plot it.

Heatmap

Tiles are filled by the ATT(g,t) estimate. Cells whose pointwise CI excludes zero are marked with a filled dot (●); when bootstrap data are available, simultaneously significant cells also receive an open diamond (◇).

plot(att_gt(res_cs), type = "heatmap")

The vertical dashed lines mark each cohort’s treatment onset (t = g).

Facet plot

One panel per cohort showing ATT over calendar time, with a pointwise CI ribbon. Useful for inspecting heterogeneous dynamics across cohorts.

plot(att_gt(res_cs), type = "facet")


Conclusion

The fixes package streamlines event study estimation and visualization for panel data researchers. With a minimal API, multiple CI support, and robust visualization, it accelerates the workflow for dynamic treatment effect analysis.

For further details and full argument documentation, see:

?event_study
?att
?did
?plot.es_result

Happy analyzing!