Event study estimation for panel data

Estimates dynamic treatment effects by relative time (an event-study curve) on panel data. Supports the classic TWFE design (universal or staggered timing) plus modern staggered-robust estimators — Callaway-Sant'Anna, Sun-Abraham, Borusyak-Jaravel-Spiess, Wooldridge two-way Mundlak, and Deb et al. FLEX — selected via estimator=. Estimation runs on the package's internal fixed-effects OLS engine; results are returned as a tidy table with rich metadata.

Usage

event_study(
  formula,
  data,
  time,
  timing = NULL,
  unit = NULL,
  estimator = c("twfe", "cs", "sa", "bjs", "twm", "flex"),
  cluster = NULL,
  weights = NULL,
  baseline = -1L,
  lead_range = NULL,
  lag_range = NULL,
  conf_level = 0.95,
  vcov = "HC1",
  vcov_args = list(),
  control_group = c("nevertreated", "notyettreated"),
  anticipation = 0L,
  bootstrap = FALSE,
  boot_reps = 999L,
  boot_alpha = 0.05,
  boot_seed = NULL,
  group = NULL,
  trends = FALSE,
  interval = 1,
  time_transform = FALSE,
  rel_time = NULL,
  staggered = NULL,
  outcome,
  treatment = NULL,
  fe = NULL,
  covariates = NULL
)

Arguments

formula: A model formula outcome ~ treatment + covariates | fixed effects. The left-hand side may be an expression (e.g. log(y)); the part after | is optional. What the first right-hand-side term means depends on timing; see Specifying treatment. Fixed effects such as id^year require the optional {fixest} package. All variables must be columns of data.
data: A data.frame containing panel data. formula and data may come in either order, so df |> event_study(y ~ D | id + year, time = year) works.
time: Unquoted time variable (numeric or Date).
timing: Treatment timing; see Specifying treatment. Either an unquoted column giving each unit's first treatment period (NA = never treated) or a single treatment date (numeric or Date). Omit it when the formula holds a D column. A timing column must have the same value on every row of a unit; "cs", "sa", "bjs" and "twm" stop with an error naming the offending units otherwise. NA means never treated: this follows did::att_gt(), fixest::sunab(), and didimputation. If NA instead represents missing treatment timing for an otherwise treated unit, those observations are silently absorbed into the control group, which is almost certainly wrong. A D column avoids the coding question.
unit: Unquoted unit identifier. Taken from the fixed effects when possible; see Panel identifiers.
estimator: Estimation strategy: "twfe" (default) for the classic TWFE event study, "cs" (Callaway-Sant'Anna 2021), "sa" (Sun-Abraham 2021 interaction-weighted), "bjs" (Borusyak-Jaravel-Spiess 2024 imputation), "twm" (Wooldridge two-way Mundlak), or "flex" (Deb et al. 2024, repeated cross-sections).
cluster: Cluster specification (one-sided formula like ~ id + year, a single character column name, or a vector of length nrow(data)).
weights: Observation weights (a name/one-sided formula or a numeric vector of length nrow(data)).
baseline: Integer baseline period (default -1); the reference period excluded from the estimation results.
lead_range, lag_range: Integers for pre/post windows. If NULL, determined automatically.
conf_level: Numeric vector of confidence levels (default 0.95).
vcov: VCOV type: "HC1" (default), "iid", or "hetero" on the internal engine; other types are delegated to fixest::vcov() when the optional {fixest} package is installed.
vcov_args: List of additional arguments forwarded to fixest::vcov() (requires {fixest}).
control_group: For estimator = "cs": comparison group, either "nevertreated" (default) or "notyettreated".
anticipation: For estimator = "cs": number of anticipation periods before treatment (non-negative integer, default 0L).
bootstrap: Logical; if TRUE and estimator = "cs", compute simultaneous confidence bands via the multiplier bootstrap (Algorithm 1, Callaway and Sant'Anna 2021). Adds conf_low_sim and conf_high_sim columns to the result and stores the full (g,t)-level bootstrap object as attr(result, "bootstrap"). Default FALSE.
boot_reps: Integer number of bootstrap draws (default 999). Used only when bootstrap = TRUE and estimator = "cs".
boot_alpha: Significance level for the simultaneous band (default 0.05). Note: this is independent of conf_level, which governs the pointwise delta-method CIs.
boot_seed: Integer seed for the bootstrap RNG; NULL (default) does not set a seed. Pass an integer for reproducible results.
group: Unquoted group identifier for estimator = "flex" only. Identifies which treatment group (cohort) each observation belongs to in a repeated cross-section design (\(R_{ig}\) in Deb et al. 2024). Each group must map to exactly one value of timing (or NA for never-treated groups). With a D column, each group's start is the first period with D = 1. Not used by other estimators.
trends: Logical; for estimator = "twm" only. When TRUE, adds cohort-specific linear trend regressors \(d_g \cdot t\) to the Procedure 5.1 regression (Wooldridge 2025, Section 8), allowing each cohort's counterfactual trend to deviate linearly from the common time trend. Requires at least 2 pre-treatment periods per cohort. Default FALSE.
interval: Numeric spacing of the time variable (default 1; ignored internally for Dates).
time_transform: Logical; if TRUE, creates consecutive integer time within unit. A timing column or value is then read as an index; a D column needs no adjustment.
rel_time: Optional unquoted column holding a pre-built event time (time relative to treatment), e.g. the Time_to_Treatment produced by paneltools/fect's get.cohort(). When supplied, event_study() uses this column verbatim as the event-study factor — i(rel_time, treatment, ref = baseline) — instead of computing relative time internally from time and timing; the first right-hand-side term is then the treated-group dummy. Never-treated / control rows should be NA. Requires estimator = "twfe"; cannot be combined with timing, staggered = TRUE or time_transform = TRUE. See Details for how this differs from the calendar-difference and time_transform conventions.
staggered: Logical or NULL (default). Only consulted by estimator = "twfe": TRUE means a staggered design, FALSE a single universal treatment date. With the default NULL the design follows from how treatment is given (a timing column is staggered, a single date is universal, a D column is universal when every unit starts in the same period). Pass an explicit value to override the inference.
outcome, treatment, fe, covariates: The model arguments of the older interface: the unquoted outcome, the treated-group dummy (used only by estimator = "twfe"), a one-sided fixed-effects formula, and a one-sided covariate formula. They still work in fixes 1.x, with a deprecation warning that shows the equivalent formula call; fixes 2.0.0 is planned to accept only formula. Name them when you use them: the second argument is now data.

Value

A data.frame of class "es_result" with columns:

term, estimate, std.error, statistic, p.value
conf_low_XX, conf_high_XX (for each requested conf.level)
relative_time (integer; 0 = event), is_baseline (logical; marks the reference period)

Attributes include: lead_range, lag_range, baseline, interval, call, model_formula, conf.level, N, N_units, N_treated, N_nevertreated, fe, vcov_type, cluster_vars, staggered, sunab_used, unit_var, time_var.

Details

The model is a formula, outcome ~ treatment + covariates | fixed effects, as in fixest::feols(); the same formula works in did() and att().

event_study() supersedes run_es() as of fixes 1.0.0.

How relative (event) time is defined. event_study() can build the event-time axis in three different ways, and in sparse or unbalanced panels (e.g. monthly surveys observed only a few times a year) they do not coincide:

Method	Definition	How to invoke
Calendar difference (default)	`(time - timing) / interval`	a timing column, a scalar `timing`, or a D column
Within-unit rank	`dense_rank(time)` within each unit	`time_transform = TRUE`
Pre-built event time	supplied column, used verbatim	`rel_time = <col>`

In a balanced panel all three agree. When observation spacing is irregular they diverge: for a unit observed at periods {1, 5, 9, 20} treated at 20, the calendar difference puts period 9 at 9 - 20 = -11, whereas a global observed-period grid (as used by paneltools/fect's get.cohort()) places it two observed waves before treatment (-2). To reproduce an analysis built on such a grid, pass that grid's event-time column via rel_time so the design matches fixest::feols(y ~ i(rel_time, treatment, ref) | fe) exactly. att() uses the same calendar-difference convention as the default here. With a D column, a unit's start is the first period in which D = 1 is observed.

Missing values and unbalanced panels. With estimator = "cs", rows with a missing or non-finite outcome, time, or unit are dropped, and the panel is then balanced by dropping every unit not observed in all periods, as did::att_gt() does by default (allow_unbalanced_panel = FALSE). Both steps warn, and N, N_units, N_treated and N_nevertreated describe the balanced sample. The regression-based estimators drop rows with missing values from the fit as fixest::feols() does; with estimator = "sa" or "twm" the cohort shares used to aggregate cohort-specific effects count only units that remain in the estimation sample (Sun and Abraham 2021, Section 4.1).

Specifying treatment

Describe treatment in whichever way matches your data:

Your data has	Write
D (1 for treated units in post-treatment periods)	`y ~ D + x \| id + year, time = year`
each unit's first treated period (`NA` = never treated)	`y ~ x \| id + year, time = year, timing = g`
a treated-group dummy and one treatment date	`y ~ treated + x \| id + year, time = year, timing = 2005`

D column (no timing): the first right-hand-side term is D. A unit's treatment starts in the first period with D = 1 and must stay on from then on (staggered adoption). Rows with a missing D are dropped. A single start date across units gives the classic universal design, several give the staggered one. With time_transform = TRUE the start is read on the time index, so no index arithmetic is needed.
Timing column (timing = g): g defines treatment, so the right-hand side holds only covariates (write y ~ 1 | id + year when there are none). D or the treated-group dummy must not appear in the formula.
Treated-group dummy and a date (timing = 2005): the first term is the treated-group dummy (constant within units), and treatment starts at timing for every treated unit.

A pre-built event-time column goes in rel_time, with the treated-group dummy as the first term: y ~ treated | id + year, time = year, rel_time = ttt.

Panel identifiers

time is required. unit is taken from the fixed effects when they are plain columns made of time and exactly one other column (| id + year with time = year, in either order); a message then reports the unit and time used. Otherwise pass unit: it is needed for a D column and for estimator = "cs", "sa", "bjs" and "twm". print() shows the identifiers of every result.

Fixed effects and covariates by estimator

"twfe", "sa" and "twm" use the fixed effects as written. "cs", "bjs" and "flex" absorb the unit (group) and time effects themselves: those fixed effects are accepted silently, any others are ignored with a warning. Covariates are supported by "twfe", "twm" and "flex", and may be expressions such as log(x) or factor(region); "cs", "sa" and "bjs" stop with an error when the formula contains covariates.

Key Features

One formula for the model (outcome, treatment, covariates, fixed effects), shared with did() and att().
Six estimators behind one interface: "twfe", "cs", "sa", "bjs", "twm", "flex".
Flexible clustering, weights, and VCOV choices (e.g., vcov = "HC1" | "iid" | "hetero" ...).
Automatic lead/lag window detection and customizable baseline period.
Returns an "es_result" object compatible with plot(), print(), and autoplot().

Examples

if (FALSE) { # \dontrun{
# D column: 1 for treated units in post-treatment periods
res <- event_study(y ~ D | id + year, data = df, time = year)
plot(res)

# Timing column (NA = never treated), Callaway-Sant'Anna
res <- event_study(y ~ 1 | id + year, data = df, time = year, timing = g,
                   estimator = "cs")

# Treated-group dummy and a single treatment date, clustered by unit
res <- event_study(y ~ treated + x1 | id + year, data = df, time = year,
                   timing = 2005, cluster = ~ id)
} # }