CS estimator: character / factor unit IDs produced
silently wrong results (run_es(estimator = "cs")
and calc_att(estimator = "cs")): unit IDs were coerced with
as.integer(), which turns every character ID into
NA — all units collapsed into one inside the Rcpp kernel
and the ATT(g,t) table was meaningless, with no error raised. Unit IDs
are now densely re-encoded with match(); character, factor,
and arbitrary numeric IDs give results identical to integer
IDs.
CS estimator: a cohort with timing == 0
leaked into the control group: 0L was the internal
never-treated sentinel, so a genuine cohort first treated at time 0
(e.g., event-time-centred data) was silently used as never-treated
controls. time/timing are now shifted by a
common internal offset so the sentinel can never collide; results are
invariant to common shifts of the time axis.
Classic staggered TWFE dropped the entire never-treated
control group (estimator = "twfe",
staggered = TRUE): the relative-time variable is
NA for never-treated units and fixest::feols()
removes rows with NA regressors, so every control row was excluded from
the estimation sample (visible only as a fixest NOTE). Those rows are
now assigned the (excluded) baseline level — their event dummies are
identically zero because the i() interaction is multiplied
by treatment = 0 — so the control group is correctly
retained. N_treated / N_nevertreated metadata
still reflect the original NA pattern.
method = "sunab" dropped the entire
never-treated control group: fixest::sunab() has
no NA-cohort convention and drops those rows. Never-treated units are
now recoded internally to a cohort far beyond the sample period (the
fixest::base_stagg convention) and retained as controls.
Verified numerically identical to a direct
fixest::feols(y ~ sunab(...)) call using the same
recoding.
Rcpp (unit, time) lookup key: replaced the
multiplicative hash (which could collide for time values above
1,000,003, e.g. Unix timestamps) with a collision-free 64-bit
(unit << 32) | time key in
src/compute_att_gt.cpp.
g - 1 base-period assumption), instead of
the generic “No ATT(g,t) estimates were computed”.g - 1 - anticipation) is not observed.time /
timing values now raise an error instead of being silently
truncated by as.integer() (new shared helper
.assert_integerish()).baseline
reference period is not feasible for any cohort, in which case no period
would be excluded and the regression is unnormalised.match() on keys instead of
per-cohort scans).tests/testthat/test-robustness.R (19 assertions)
covering all of the above: ID-type invariance, time-shift invariance,
control-group retention (classic TWFE and sunab vs. direct fixest
reference), grid-spacing errors, fractional-time rejection, and BJS
singleton imputation.run_es() — diagnostic warning for ambiguous
timing = NA (classic TWFE,
staggered = TRUE only): a warning() is now
emitted when any unit has treatment = 1 in at least one row
but timing = NA. Such units are silently absorbed as
never-treated controls — correct when NA is intentional,
but dangerous when NA represents genuinely missing
treatment timing. The warning names up to 5 affected unit IDs (when
unit is supplied) for quick diagnosis.
@param timing documentation
clarified: all staggered estimators (cs, sa, bjs, twm, flex,
classic with staggered = TRUE) share the same
NA = never treated convention as
did::att_gt() and fixest::sunab(). The roxygen
documentation now states this explicitly and notes the risk of silent
misclassification when NA represents missing data rather
than “never treated”.
honest_sensitivity() — honest robust
inference (Rambachan & Roth 2023) (R/honest_did.R): sensitivity analysis for
event-study / DiD designs when parallel trends may be violated. Instead
of assuming parallel trends holds exactly, it reports confidence sets
for a post-treatment effect under progressively weaker restrictions on
the difference in trends, and a breakdown value (the largest
restriction at which the effect is still significant):
res <- run_es(df, outcome = y, treatment = treat, time = year, timing = 6,
fe = ~ id + year)
h <- honest_sensitivity(res, type = "relative_magnitude",
Mvec = c(0, 0.5, 1, 1.5, 2))
plot_honest(h)"relative_magnitude" () and
"smoothness" ().run_es() es_result
(estimator "twfe": classic or
method = "sunab", which now carry the event-study
coefficient covariance via the new es_vcov attribute), or
from a raw betahat / sigma pair for any
estimator.plot_honest() / autoplot() draw the
“top-down” sensitivity plot; print.honest_result() reports
the breakdown value."Conditional"): (single post period) matches to
machine precision; matches to grid resolution. Heavy numeric helpers
(lpSolveAPI, Rglpk,
TruncatedNormal, Matrix, pracma)
are in Suggests and gated with
requireNamespace().Estimator-core consolidation
(behaviour-preserving): factored the boilerplate duplicated across
run_es(), calc_att(), run_did(),
and the five estimator backends into shared helpers in R/utils-internal.R —
.resolve_col() (NSE column resolution),
.add_ci_columns() (CI construction),
.validate_panel_cols(),
.compute_cohort_sizes(), .model_vcov_full(),
and .aggregate_iw(). The calc_att()
aggregation loops were vectorised. All existing tests pass
unchanged.
run_did() — basic TWFE DiD
estimator (R/run_did.R): New function
for classic two-way fixed-effects difference-in-differences. Accepts a
pre-built binary treatment indicator D_it and returns a
did_result S3 object with full modelsummary /
tinytable compatibility:
df$D <- as.integer(df$treated & df$year >= 2006)
result <- run_did(df, outcome = y, treatment = D, fe = ~ id + year)
print(result)
modelsummary::modelsummary(result) # tinytable output (modelsummary >= 2.0)run_es() and
calc_att() — specify variables by name, not by formula
construction.fe auto-inferred as ~ unit + time when
both unit and time are supplied and
fe = NULL.tidy.did_result() delegates to
broom::tidy.fixest() for full regression table (all
regressors); glance.did_result() delegates to
broom::glance.fixest() for model-level stats
(within.r.squared, nobs, AIC,
etc.).print.did_result() shows a clean header: estimator,
sample sizes, FE spec, VCOV type.log(y) supported.modelsummary and tinytable added to
Suggests.calc_att() — ATT aggregation
function (R/calc_att.R): New
function that separates “calculating ATT” from “running an event study”.
While run_es() produces a full dynamic event-study curve
(effects by relative time), calc_att() yields aggregated
ATT estimates with a clear, minimal interface:
calc_att(data, outcome = y, time = year, timing = g, unit = id,
estimator = "cs",
aggregation = "simple") # overall ATT
calc_att(..., aggregation = "by_cohort") # ATT per cohort
calc_att(..., aggregation = "by_time") # ATT per calendar period"cs" (Callaway-Sant’Anna
2021) and "bjs" (Borusyak et al. 2024). SA, TWM, FLEX are
event-study-only estimators and are not supported (use
run_es()).did::aggte() for point estimates:
"simple": exposure-weighted average over all
post-treatment (g,t) pairs."by_cohort": mean ATT(g,t) over post-treatment periods
per cohort."by_time": cohort-size-weighted ATT per calendar
period.did::aggte SEs,
which use influence-function variance accounting for within-cohort time
correlation. Cluster-robust SE for BJS aggregations planned for
v0.12.att_result S3 object
(data.frame subclass) with columns group,
estimate, std.error, statistic,
p.value,
conf_low_XX/conf_high_XX. Attributes:
aggregation, estimator, N,
N_units, N_treated, control_group
(CS only), att_gt (raw CS ATT(g,t) table),
tau_it (BJS unit-time effects table).print.att_result() S3 method provides
a concise summary.tests/testthat/test-att.R:
did::aggte() at tolerance
1e-4mean(tau_it$tau_hat) at
tolerance 1e-12conf.level support verifiedrun_es.R refactoring: Extracted
.make_tidy() and .es_finalize() private
helpers that centralise the repeated tidy-construction / baseline-row /
range-filter / attr-stamping pattern that was previously duplicated
across the cs, sa, bjs,
twm, and flex branches (~360 lines removed).
No API change; all 593 tests pass.inst/profile/benchmark.R with
profile_estimator() and benchmark_all()
utilities. Run interactively to generate flamegraph HTML and wall-clock
summary across estimators before starting Phase 3 Rcpp work.src/compute_att_gt.cpp (shipped in v0.8.1 but missing from
the Rcpp Acceleration Roadmap), and updated Package Structure to list
all 13 R source files, 5 Rcpp files, and 16 test files that now
exist.test-twm.R Test 16 — trends=TRUE with
exactly 2 pre-treatment periods succeeds without warning (boundary
condition of the ≥ 2 requirement).test-twm.R Test 17 — explicit
lead_range/lag_range correctly trims TWM
output and preserves estimates in the shared window (exercises
.es_finalize() filter).test-flex.R Test 12 — explicit
lead_range/lag_range trims FLEX output
correctly.test-flex.R Test 13 — FLEX succeeds with an all-treated
panel (no never-treated groups; N_nevertreated == 0).test-sa.R — explicit
lead_range/lag_range trims SA output and
preserves estimates within the shared window.run_es(estimator = "twm", trends = TRUE):
Cohort-specific linear trend detrending (Wooldridge 2025, Section 8).
Adds d_g * t regressors so each cohort’s counterfactual
trend can deviate linearly from the common time trend. Implemented as
post-treatment-only cells + trend columns; requires >= 2
pre-treatment periods per cohort. Output shows
relative_time >= 0 only (pre-trend testing is deferred
to the no-trend model).run_es(estimator = "twm", covariates = ~ x1 + x2): Full
Wooldridge (2025) Procedure 5.1 covariate interactions. Adds
treatment-cell × cohort-demeaned covariate terms
(ẋ_{ig} = x_i - x̄_g) and i(time, x_j)
conditional parallel-trends controls (Eqs. 5.2-5.3).run_es(estimator = "flex", covariates = ~ x1 + x2):
Full Deb et al. (2024) Eq. 3.1 covariate interactions for RCS data.
Cell-level centering X_{i,t} - X̄_{g,t} (Eq. 2.11) with
i(time, x_j) and i(group, x_j) conditional PT
controls.src/indicator_matrix.cpp —
build_indicator_matrix_cpp(): Replaces the pure-R nested
for-loop that fills the 0/1 indicator matrix in the SA, TWM, and FLEX
estimators. Single shared Rcpp utility, O(N×K) integer fill with no R
allocation overhead per column.src/iw_aggregation.cpp —
aggregate_iw_cpp(): Replaces the R event-time aggregation
loop and the t(w) %*% V_sub %*% w quadratic-form VCOV step
in SA, TWM, and FLEX. Uses RcppArmadillo
arma::mat::submat() for O(1) VCOV block extraction and
arma::as_scalar(w.t() * V_sub * w) for the quadratic
form.src/cov_demeaning.cpp —
build_cov_interactions_cpp(): Replaces the pure-R covariate
demeaning loops in TWM (cohort-level centering) and FLEX (group×time
cell-level centering) with a single
std::unordered_map-based O(N) grouping pass followed by
column-wise mean subtraction, then builds the N × (K×p) treatment-cell ×
centred-covariate interaction matrix in C++.src/bootstrap_cs.cpp — bootstrap_cs_cpp():
Replaces the R-level apply/sweep calls in the
CS multiplier bootstrap with RcppArmadillo DGEMM for both the pilot and
main stages. Uses R::runif() for Mammen weight generation
(column-major fill order matches R’s matrix()) so results
are numerically identical to the pure-R path under the same
set.seed().src/Makevars and src/Makevars.win
linking $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) for
cross-platform RcppArmadillo support.RcppArmadillo added to LinkingTo in
DESCRIPTION.run_es(estimator = "twm"): Wooldridge (2025) Two-Way
Mundlak (TWM) estimator for panel data. Implements Procedure 5.1 via
POLS on cohort x calendar-time treatment cells with two-way FE.
Algebraically identical to estimator = "sa" in the
no-covariate base case (verified numerically at tolerance 1e-10). The
trends = TRUE extension (cohort-specific linear trend
detrending, Wooldridge 2025 Section 8) is reserved for v0.9.1.run_es(estimator = "flex", group = <group_id>):
Deb, Norton, Wooldridge & Zabel (2024) FLEX estimator for
repeated cross-section (RCS) data. Uses group x
calendar-time OLS with group+time FE (no unit FE). Algebraically
equivalent to the multi-step imputation estimator (Proposition 2.1). New
group argument identifies the treatment group (R_{ig}) each
individual belongs to. Default clustering is by group. Covariate
interactions deferred to v0.9.1.papers/Deb et al. (2024).txt,
papers/Woodridge (2025).txt.run_es(estimator = "cs"): Callaway-Sant’Anna (2021)
estimator. Computes group-time ATT(g,t) via the unconditional DiD
estimand (eq. 2.8) with never-treated or not-yet-treated control groups.
Supports control_group = "nevertreated" (default) or
"notyettreated". The full ATT(g,t) matrix is stored as the
att_gt attribute of the result.run_es(estimator = "sa"): Sun-Abraham (2021)
interaction-weighted estimator. Aggregates cohort x relative-time
interactions by cohort share weights. Numerically identical to
fixest::sunab() to machine precision.run_es(estimator = "bjs"): Borusyak, Jaravel &
Spiess (2024) imputation estimator. Fits TWFE on untreated observations
only, imputes counterfactuals for treated observations, and averages
treatment effects by horizon. Handles singleton unit fixed effects via
closed-form recovery.run_es(bootstrap = TRUE): Multiplier bootstrap for
simultaneous confidence bands (Algorithm 1, Callaway & Sant’Anna
2021). Adds conf_low_sim / conf_high_sim
columns to the result and stores the (g,t)-level bootstrap object as
attr(result, "bootstrap"). Controlled by B
(draws, default 999) and boot_seed arguments. Available
only when estimator = "cs".plot_att_gt(): Visualize the full ATT(g,t) matrix from
CS results as a heatmap (type = "heatmap") or
cohort-faceted time series (type = "facet"). Requires
estimator = "cs" result as input. When
bootstrap = TRUE was used, the heatmap adds a subtitle with
the simultaneous critical value and open-diamond markers for
simultaneously significant cells; the facet adds a lighter simultaneous
CI ribbon.plot_es(show_simultaneous = TRUE): Overlays the
simultaneous bootstrap CI (lighter band, alpha 0.15) alongside the
pointwise CI (alpha 0.3) with a two-entry legend. Errors informatively
if bootstrap was not run.plot_es_interactive(show_simultaneous = TRUE): Adds a
second lighter ribbon trace for the simultaneous CI and extends the
hover tooltip with simultaneous CI bounds.run_es() bootstrap path: base::merge()
silently dropped all custom attributes (e.g. att_gt,
N_units, lead_range) from the
es_result object. Attributes are now saved and restored
around the merge, so plot_att_gt() works correctly after a
bootstrap = TRUE run.did and didimputation to Suggests
for numerical agreement tests against reference implementations.cluster in
run_es() did not produce clustered standard errorsvcov = "HC1" was being applied after model
estimation, overriding the clustered SE from feols()cluster is specified and vcov is
left at its default ("HC1"), the model’s clustered standard
errors are usedcluster is set,
pass vcov = "HC1" together with
cluster = NULLclassic and sunab
methodsrun_es() results
and equivalent direct fixest::feols() callslead_range and lag_range
parameters were ignored, causing all estimated coefficients to be
returned regardless of specified ranges[-lead_range, lag_range]classic and sunab
methods'sunab' is not an exported object from 'namespace:fixest'
error when using method = "sunab"sunab function from fixest packagebaseline parameter (default: -1) now applies to
both classic and sunab methodsestimate = 0,
std.error = 0, and is_baseline = TRUEbaseline is used
for both methods"-9", "0", "3")"year::-9") for better readability%>%) with base R pipe
(|>) throughout the package|>)relative_time was NA for all
coefficients except baseline when using non-staggered treatment
timingfixest::i() now correctly
corresponds to the baseline parametertiming = 5 and
baseline = -1, period 4 (not period 5) is now used as the
referenceplot_es_interactive() for creating
interactive event study plotses_result objects from
run_es()sapply()
approach for significant performance gainsplotly package (optional,
suggested dependency)run_es() and plot_es() now support returning
and visualizing multiple confidence levels (e.g., 90%, 95%, 99%) in a
single analysis.plot_es() adds ci_level selection, more theme
options (theme_style), and improved ribbon/error bar
display.NA in
timing):
timing are now retained as
never-treated controls (staggered only).weights input:
~ popwt), bare names
(popwt), or character strings ("popwt").cluster input handling:
unit is supplied without
time_transform = TRUE.staggered = TRUE option allows timing
to vary by unit (e.g., treatment year column).NA in timing are safely
retained as untreated.weights argument (e.g., ~ popwt)
to run weighted regressions.lead_range or lag_range is
NULL, the function computes the maximum feasible range from
the data.unit is specified without
time_transform = TRUE.plot_es():
ggplot2::scale_x_continuous() with
integer breaks spaced by 1, aligned to the relative_time
range.Date:
time_transform = TRUE to automatically
convert the time variable into a unit-level sequential
index (1, 2, 3, …) for event study estimation.unit argument to specify the panel unit
identifier required when time_transform = TRUE.Date class in the time
variable and converts it automatically to numeric if
time_transform = FALSE.unit is missing or time is of unsupported
type.@examples in the function documentation to
include Date-based examples.time_transform
usage.README.md to describe irregular time handling
and demonstrate new use cases.time_transform,
unit handling, and Date conversion edge
cases.lead1, lag0) already exist in the dataset to
prevent accidental overwriting.lead_range, lag_range, and
interval) has fewer than 10 rows, helping users identify
overly narrow estimation windows.treatment variable: it is now
coerced to logical using as.logical() to support both
binary numeric (0/1) and logical (TRUE/FALSE)
formats.fe argument
(e.g., ~ id + year) were combined using
model_formula | fe_text, which caused evaluation errors
during tests.as.formula() to ensure compatibility with
fixest::feols().run_es():
~ x1 + x2).fe and cluster arguments must now be
specified using a one-sided formula (e.g.,
~ id + year).cluster is still
accepted.fe_var argument now supports additive notation
(firm_id + year) instead of character vectors.plot_es() efficiency and documentation.fe notation.This version introduced several enhancements and refinements to improve usability and maintainability.
outcome_var, treated_var, and
time_var are now processed using
rlang::ensym() for better robustness.fe_var and cluster_var handling improved
for more reliable column referencing.plot_es() function:
relative_time, estimate, etc.) are
present.baseline handling could lead
to incorrect sorting of lead/lag terms.This is the first release of the fixes package,
providing tools for estimating and visualizing event study models with
fixed effects.
run_es(): A function to estimate event
study models using fixest::feols(), generating lead and lag
variables automatically.
fe_var as
character vector).cluster_var.interval argument.plot_es(): A function to visualize
event study results with ggplot2.
type = "ribbon", default).type = "errorbar").fixest::feols().interval
argument).c("firm_id", "year"))."state_id").fe_var.