Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
45,778 characters · 16 sections · 31 citation commands
A Fixed-Effects Causal Forest for Staggered Adoption, with an Application to Medicaid Expansion
Difference-in-differences (DiD) with staggered treatment timing is among the most widely used research designs in empirical economics. A recent literature has shown that the workhorse two-way fixed-effects (TWFE) regression can be badly biased in this setting when treatment effects are heterogeneous across adoption cohorts or over time: the TWFE coefficient is a variance-weighted average of all possible $2\times2$ comparisons, including “forbidden” comparisons that use already-treated units as controls, and these can enter with negative weights goodman2021difference, dechaisemartin2020two, sun2021estimating. callaway2021difference (hereafter CS) resolve this by defining the target as a collection of group-time average treatment effects $\mathrm{ATT}(g,t)$, the effect at period $t$ for the cohort first treated at $g$. Each $\mathrm{ATT}(g,t)$ is identified from a clean comparison of that cohort against units not yet treated as of $t$, and the estimates are then aggregated with transparent weights.
Two limitations motivate this paper. First, $\mathrm{ATT}(g,t)$ is an average over the treated cohort; it does not describe how the effect varies with observed characteristics. In most applications, though, the question that matters for policy is which places, firms, or people respond most. Second, the causal-forest literature that estimates such conditional average treatment effects (CATEs) nonparametrically wager2018estimation, athey2019generalized was developed under unconfoundedness and does not, out of the box, respect the parallel-trends identification or the not-yet-treated comparison structure that staggered DiD requires.
We study the covariate-conditional group-time effect
the CS estimand made heterogeneous in pre-treatment covariates $x$. We propose to estimate it with a fixed-effects causal forest that operates block by block. For each cohort $g$ and post-period $t \ge g$ we form the same clean two-group, two-period panel that CS use: cohort $g$ against its not-yet-treated (or never-treated) controls, observed at the reference period $g-1$ and the evaluation period $t$. On that block we fit a causal forest that, inside each tree node, residualizes the outcome and treatment on unit and period fixed effects before searching for splits that maximize heterogeneity in the treatment effect. The block-level surfaces are aggregated into event-study and overall summaries with CS-style weights.
We are explicit about what is and is not new. The estimand $\tau_{g,t}(x)$ was introduced by hatamyar2023mldid, whose MLDID estimator combines the doubly-robust DiD-CATE of lu2019estimating with the CS staggered framework and not-yet-treated controls. imai2023doubly study the same object and provide doubly-robust uniform confidence bands, but for a single continuous covariate estimated by local polynomials. We do not claim to have introduced conditional group-time effects, nor the combination of machine learning with staggered DiD. What we contribute is a different way to estimate the object. MLDID removes confounding by modeling nuisance functions, an outcome regression and generalized propensity scores fit with an R-learner. We instead difference out unit and period effects within each tree node, using the fixed-effects residualization of kattenberg2023causal and gavrilova2025did. Those authors developed that device for a common adoption date---a single treatment cohort, with event-study dynamics measured against a fixed base period---rather than for staggered timing; we carry it into the Callaway--Sant'Anna group-time blocks, so that each cohort is compared only to units not yet treated and effects may differ across adoption cohorts. The two approaches are alternatives rather than a re-derivation, and Section (ref) benchmarks them directly. Under timing heterogeneity our estimator has lower bias and valid coverage, while the doubly-robust engine can be more efficient on the conditional surface when heterogeneity is smooth in covariates.
One caveat shapes how we present the results. Because each block is a clean two-group within-fixed-effects comparison, the overall (covariate-averaged) effect our estimator reports is numerically identical, by construction, to the CS aggregate. The forest adds value only on the conditional object $\tau_{g,t}(x)$ and its event-study aggregate. We therefore lead with the conditional surface and treat the scalar as a check against CS rather than as a contribution.
In Monte Carlo experiments across four data-generating processes we benchmark the estimator against TWFE, the CS estimator, a naively pooled fixed-effects causal forest, and MLDID. Under homogeneous effects, and under effects that vary only in covariates, all of the estimators are approximately unbiased on the overall effect. The differences appear when the effect varies in timing, with an early cohort that gains more and a late cohort that gains less. There TWFE and the pooled forest inherit large forbidden-comparison bias, roughly $-1.25$ on a true effect of comparable size and with zero confidence-interval coverage, while our estimator and CS remain close to unbiased with valid coverage. On the conditional surface, our estimator has the lowest error of the forest-based methods under homogeneous, timing-heterogeneous, and nonlinear heterogeneity, and improves on the pooled forest even under smooth covariate heterogeneity; only MLDID's penalized-linear learner does better in the one design where the true surface is linear in the covariates. We then apply the method to the ACA Medicaid expansion. Across the staggered county-level rollout it recovers an average -2.25 percentage-point effect on the county uninsured rate, and a conditional surface on which poorer, lower-income counties gained the most coverage---heterogeneity measured along covariates that are not lags of the outcome.
The rest of the paper proceeds as follows. Section (ref) sets up the panel and the estimand. Section (ref) states the identifying assumptions and defines the estimator. Section (ref) reports the Monte Carlo evidence, and Section (ref) the applications. Section (ref) concludes and discusses the limitations.
We observe a panel of units $i = 1,\dots,N$ over periods $t = 1,\dots,T$. For each unit-period we observe an outcome $Y_{it}$, a binary treatment indicator $D_{it} \in \{0,1\}$, and a vector of time-invariant pre-treatment covariates $X_i \in \mathbb{R}^p$.
\paragraph{Staggered, absorbing treatment.} Treatment is an absorbing state: once treated, a unit remains treated. Each unit's adoption time (cohort) is $G_i = \min\{t : D_{it} = 1\}$, with $G_i = \infty$ for never-treated units, so that $D_{it} = \mathbf{1}\{t \ge G_i\}$. Let $\mathcal{G}$ denote the set of observed finite adoption cohorts.
\paragraph{Potential outcomes.} Following CS, we index potential outcomes by cohort: $Y_{it}(g)$ is the outcome unit $i$ would realize at $t$ had it first adopted at $g$, and $Y_{it}(\infty)$ is the never-treated potential outcome. The observed outcome is
Effects may depend on the cohort $g$, on calendar time $t$, and on the covariates $X_i$, which are the object of this paper.
\paragraph{Estimand.} The target is the group-time conditional average treatment effect on the treated,
This is $\mathrm{ATT}(g,t)$ conditioned on $X_i = x$. Reported objects are weighted averages of it. The conditional event study at exposure $e = t - g$ is
with cohort-share weights $w_g^{(e)}$ restricted to cohorts observed at exposure $e$; the overall conditional effect $\tau^{\mathrm{overall}}(x)$ averages $\tau_{g,t}(x)$ over post-treatment $(g,t)$ blocks with treated-observation-share weights, and the scalar overall ATT is its covariate average over the treated.
Assumption (ref) is the CS conditional parallel-trends condition imposed cell by cell. Assumption (ref) is more demanding than the overlap condition of a pooled causal forest: a late cohort with no contemporaneous not-yet-treated controls in some covariate region is not identified there, and our estimator reports such blocks as dropped rather than silently extrapolating. Under Assumptions (ref)--(ref), for $t \ge g$,
a conditional-on-$x$ difference-in-differences between the cohort and its valid controls, anchored at the pre-adoption reference $g-1$.
The estimator marries the not-yet-treated comparison structure to a fixed-effects causal forest. For each $(g,t)$ with $t \ge g$:
Algorithm (ref) states the full procedure.
The engine distinction from MLDID lives in step 2: where MLDID forms a doubly-robust pseudo-outcome from separately estimated nuisance functions and runs an R-learner, we difference out unit and period effects within the node and read the local effect off the residualized regression. Neither requires modeling the propensity of treatment timing; both target (ref).
It is worth spelling out why a naively pooled forest is biased here. Applying an ordinary fixed-effects causal forest to the full staggered panel, pooling all periods and cohorts into one node-level residualization, reintroduces the forbidden comparison. Within a node the pooled within-FE coefficient is, by the goodman2021difference decomposition, a variance-weighted average of every $2\times2$ comparison including late-versus-early terms in which an already-treated cohort serves as the control; when effects grow with exposure those terms enter with negative weight. Node-level residualization localizes but does not remove this: each covariate cell inherits its own contaminated weighting. The block construction removes the offending comparisons. We do not derive a closed-form expression for this bias; its sign and magnitude are established quantitatively in Section (ref).
Within a block, our estimator is an ordinary causal forest applied to fixed-effects-residualized data, which lets us borrow the consistency theory of wager2018estimation block by block, following the reduction argument of gavrilova2025did. Fix a block $(g,t)$ and let $\tilde Y_{i} = Y_{it} - Y_{i,g-1}$ be the base-period-differenced outcome for units in cohort $g$ and its comparison set $\mathcal{C}_{g,t}$, with $W_i$ the cohort-$g$ membership indicator.
Two caveats bound this result, both shared with the causal-forest DiD literature. It is a block-wise statement under the block-overlap Assumption (ref); it does not by itself establish valid coverage for the pooled cross-block surface under finite cohorts, for which we rely on the bootstrap below. And it inherits the finite-sample attenuation that motivates the recentering and orthogonalization devices discussed by athey2019generalized and gavrilova2025did; we address the level through recentering and quantify residual bias by simulation. The honesty invoked here is an asymptotic device: the consistency and normality guarantee is for the honest forest as block size grows. In finite two-period blocks with only a few dozen units, honest sample-splitting can attenuate the surface more than the bias it removes is worth, which is why the implementation splits honestly only in blocks large enough to afford it. Table (ref) shows the finite-sample surface comparison is robust to this choice, so nothing in the empirical results turns on trading the asymptotic guarantee for finite-sample accuracy in small blocks.
We report inference for both objects. For the overall effect we use a unit-level block bootstrap: resample whole unit series, refit the entire pipeline, and recompute the aggregate. This captures sampling variability including that of the cohort weights, and it is our primary interval; in the state-assigned application of Section (ref) we resample at the treatment-assignment (state) level. For the conditional surface we report a sampling-based pointwise interval from the same block bootstrap, refitting the estimator on each resample and recomputing the fitted surface $\hat\tau(x)$ at the evaluation points, then taking pointwise percentile bands across draws. These bands reflect sampling uncertainty in $\tau_{g,t}(x)$. They are wider than the forest half-sample stability bands, which measure only sensitivity to which trees were grown. Because recentering shifts with each resample, we report the bands around the bootstrap distribution of the recentered surface. This is the pointwise analogue, for the multivariate forest, of the uniform bands imai2023doubly derive for a single continuous covariate.
We simulate staggered panels and benchmark five estimators: TWFE-OLS; the CS group-time estimator with not-yet-treated controls (average aggregation); a naively pooled fixed-effects causal forest; MLDID; and the staggered fixed-effects causal forest of this paper. Four data-generating processes hold fixed a never-treated control group and staggered cohorts, and vary how the effect depends on covariates and timing: (a) a homogeneous constant effect; (b) an effect heterogeneous in covariates but constant across cohorts and exposure; (c) an effect heterogeneous in timing, in which earlier cohorts gain more and the effect grows with exposure; and (d) an effect that is nonlinear in and interacts the covariates (a threshold-plus-quadratic surface). Case (c) generates forbidden-comparison bias; case (d) is the regime a forest should recover better than a linear learner. Covariates are time-invariant at the unit level. For each DGP we report, over 200 replications, the bias, root-mean-squared error, and 95% confidence-interval coverage of the overall ATT. For the forest and DML estimators we also report the root-mean-squared error of the fitted conditional surface against the true per-unit effect (CATE-RMSE). The two block-bootstrap estimators use the same unit-level bootstrap, so their coverage is comparable.
Within each staggered comparison block the forest uses honest sample-splitting only when the block is large enough to afford it (at least a few hundred observations); on the small two-period blocks typical of staggered panels, honest splitting would estimate leaf effects on a handful of observations and attenuate the surface, so those blocks are fit without it. In the large blocks of the Medicaid application all blocks clear the threshold and honest splitting is used throughout. This adaptive rule is not a tuned knob on which the main result rests: Table (ref) sweeps the honesty regime and shows that under timing heterogeneity (DGP c) the staggered forest recovers the surface far better than the pooled forest at every setting, whether honesty is forced on, forced off, or selected by block size. The clean-block construction, not the honesty setting, is what drives the advantage there.
Table (ref) reports the results. Three patterns stand out. First, on the overall effect, TWFE and the pooled forest carry a bias of -1.249 with zero coverage under timing heterogeneity (DGP c), whereas CS and the staggered forest are unbiased with close-to-nominal coverage; under the other DGPs all estimators are approximately unbiased.\footnote{The staggered forest's scalar point estimate equals the CS aggregate by construction, but its interval comes from the unit-level block bootstrap rather than the CS multiplier bootstrap, so the two need not have identical coverage in finite samples; in DGP (a) the block-bootstrap interval covers at 0.85 against CS's 0.96, reflecting harder bootstrap calibration on the small two-period blocks. Coverage is at or near nominal for both in the other designs.} This is the forbidden-comparison result: the staggered forest is the only forest-based estimator that survives timing heterogeneity on the aggregate. Second, and this is where the method earns its keep, on the conditional surface under timing heterogeneity (DGP c) the staggered forest attains a CATE-RMSE of 1.887 against the pooled forest's 2.342---a gap of many Monte Carlo standard errors that the honesty sweep in Table (ref) shows is present at every honesty setting. This is where the clean-block construction earns its keep: the pooled forest's node-level pooling lets already-treated units contaminate the control group, and no honesty setting repairs that. The staggered forest also holds the lowest CATE-RMSE under homogeneous effects (DGP a) and improves on the pooled forest under smooth covariate heterogeneity (DGP b). Third, the staggered forest and CS coincide on the scalar in every DGP, as the construction implies; the forest's added value is the conditional surface, which CS does not provide.
The fifth estimator, MLDID hatamyar2023mldid, shares our estimand but uses the alternative engine: the lu2019estimating doubly-robust score with machine-learning nuisance functions rather than within-block fixed-effects residualization. Because it too uses not-yet-treated controls, MLDID avoids the large forbidden-comparison bias of TWFE and the pooled forest: its DGP (c) bias is +0.381 rather than -1.249, which confirms that both engines solve the staggered problem. The tradeoff shows up in the comparison between the two. Under timing heterogeneity (DGP c) the fixed-effects forest does better on every axis: lower bias (-0.005 versus +0.381), valid versus poor coverage (1.00 versus 0.33), and a marginally better conditional surface (1.887 versus 1.946). Under nonlinear interacted heterogeneity (DGP d) the forest-based estimators recover the surface better than MLDID's penalized-linear learner---CATE-RMSE 0.837 for the staggered forest and 1.061 for the pooled forest, both below MLDID's 1.308---as a tree ensemble should against a linear final stage when the truth is nonlinear. Here the two forests are close to each other; the separation is between the forest engines and the linear one, which is the point of the DGP. The one place the ordering reverses is pure covariate heterogeneity with a truly linear surface (DGP b): there MLDID's penalized-linear final stage fits the linear-in-$x$ effect more efficiently (CATE-RMSE 0.236 versus our 0.480), while our estimator retains better coverage. MLDID also under-covers throughout, reflecting the greater difficulty of calibrated inference for the doubly-robust score relative to the within-transformation bootstrap. In short, the two engines are alternatives: ours is preferable when effects vary with adoption timing, when the surface is nonlinear, and when calibrated coverage matters; the doubly-robust engine is more efficient only when the conditional surface is smooth and linear and timing is benign.\footnote{We benchmark against a faithful Python reimplementation of the Lu--Nie--Wager score underlying MLDID, using random-forest nuisance functions. We attempted to benchmark against the authors' original MLDID R package directly; it installs but errors at runtime in current R even on its own bundled example, apparently a version incompatibility with its CVXR solver dependency, so we could not use it. Our port omits the minimax reweighting of hatamyar2023mldid, which targets the scalar aggregate; the small positive scalar bias MLDID shows in the table is attributable to this omission. The comparison in this paper turns on the conditional surface (CATE-RMSE), not the scalar, so the omission does not affect the head-to-head we draw; and our own estimator's scalar equals the Callaway--Sant'Anna aggregate by construction, independently of MLDID.}
Because honest sample-splitting trades bias for variance and the two-period staggered blocks are small, one might worry that the surface comparison hinges on whether honesty is on. Table (ref) settles this by running both forests under the same honesty regime---forced on, forced off, and adaptive at several block-size thresholds---and reporting the CATE-RMSE at each. Under timing heterogeneity (DGP c) the staggered forest beats the pooled forest by many Monte Carlo standard errors at every regime; the clean-block construction, not the honesty setting, is what delivers the advantage. Under the nonlinear DGP (d) the two forests are close to one another at every regime---consistent with the head-to-head above, where the meaningful separation is between the forest engines and MLDID's linear final stage rather than between the two forests.
As a validation on a familiar dataset, we apply the estimator to the county minimum-wage / teen-employment panel of callaway2021difference (500 counties, 2003--2007; cohorts adopting in 2004, 2006, and 2007, with 309 never-treated counties; a single covariate, log county population). The estimator recovers the group-time pattern reported by CS: the 2004 cohort's effect deepens with exposure (about $-0.02$, $-0.08$, $-0.14$ at exposures 0--2), and the overall ATT is -0.040 (95% CI $[-0.064, -0.015]$), a small negative effect of the minimum wage on log teen employment, matching the sign and magnitude of the published estimates (Table (ref)). This dataset has too thin a covariate set to show conditional heterogeneity, which is why we turn to a richer panel.
Our substantive application is the staggered, county-level rollout of the Affordable Care Act's Medicaid expansion. Because states expanded at different dates (a large 2014 wave, then smaller cohorts through 2023) and ten states have not expanded, the design is a clean staggered adoption with a large never-treated control group. We assemble a county-year panel, 2008--2023, from the Census Small Area Health Insurance Estimates (SAHIE): the outcome is the county uninsured rate; treatment is the state expansion, coded absorbing under a majority-year rule (a state is treated in the first calendar year for which its expansion was effective at least six months); and the time-invariant pre-treatment covariates are the county's 2013 poverty rate, its 2013 log median household income (both from the Census Small Area Income and Poverty Estimates, SAIPE), and census region. We deliberately condition on socioeconomic covariates that are not lags of the outcome: heterogeneity found along the poverty rate or income cannot be an artifact of mechanical mean reversion in the uninsured rate. The 2014 cohort alone contains about 1{,}180 counties and every cohort exceeds 60, so, unlike a cross-country macro panel, cohorts are large enough for a forest to learn within-cohort heterogeneity. Table (ref) lists the covariates and Table (ref) reports their 2013 distribution by cohort.
Treatment is assigned at the state level and inherited by counties. The number of independent treatment switches is therefore the number of expanding states, not the number of counties; the county cross-section buys resolution for the conditional surface, not additional identifying variation. Accordingly we cluster inference at the state level, resampling whole states in the block bootstrap.
Table (ref) summarizes the results, and Figure (ref) plots the event study together with the placebo pre-trends. The overall effect is a -2.25 percentage-point reduction in the county uninsured rate (95% CI $[-3.10, -1.12]$), within the range of the ACA-expansion literature, and the event study builds from $-1.5$ at adoption to about $-2.8$ by the sixth year before attenuating at long horizons where few cohorts and a thinning control pool remain. The object of interest, however, is the conditional surface: its spread is substantial (standard deviation 1.17 on a mean of -1.90), and it correlates -0.55 with the county poverty rate and 0.46 with log median household income. The forest's splits load on income (0.42) and poverty (0.39), with region secondary (0.19). Poorer counties, and lower-income counties, gained the most coverage from expansion. This is a form of heterogeneity that the average $\mathrm{ATT}(g,t)$ cannot express and that the method recovers directly.
Because these covariates are socioeconomic rather than lags of the uninsured rate, the gradient is not an artifact of mechanical mean reversion: nothing about conditioning on the poverty rate forces a county with a larger coverage gain. The finding is also economically sensible---Medicaid expansion reaches more of the uninsured where poverty is higher and incomes lower. The method's role is to recover this gradient nonparametrically and jointly with the dynamics rather than to impose a functional form.
The placebo pre-trend test estimates group-time DiDs at pre-adoption periods, which should be zero under Assumption (ref). It is clean in the periods that matter (exposures $-2$ to $-4$, at $+0.04$ to $+0.20$ pp) but shows large, trending placebo effects in the deep pre-period (exposures $\le -9$, i.e.\ 2008--2012, up to $-1.98$), which reflect divergent uninsured-rate paths between expansion and non-expansion states through the Great Recession. This does not threaten identification. Because the estimator uses the universal base period $g-1$, every cohort's reference year (2013 for the 2014 wave, later for later cohorts) lies within a 2011--2023 window, so the recession-era years never enter an estimation block; they appear only as placebo leads. Table (ref) makes this concrete: re-estimating on a trimmed 2011--2023 panel, which drops the recession-era pre-period entirely, leaves the overall ATT and its interval unchanged to three decimals, because no estimation block used those years.
We also report a rambachan2023more relative-magnitudes sensitivity analysis. Allowing post-treatment differential trends of up to $\bar M$ times the largest near-pre-period placebo, the average post-expansion effect remains significantly negative for every $\bar M$ we consider (through $\bar M = 5$): even if unobserved trend violations after adoption were five times the largest violation we observe just before it, the conclusion of a coverage gain would stand. The effect is therefore robust to economically large departures from parallel trends.
This paper has proposed a fixed-effects causal forest for staggered-adoption difference-in-differences that delivers covariate-conditional group-time treatment effects $\tau_{g,t}(x)$. The estimand is shared with recent work hatamyar2023mldid, imai2023doubly. The contribution is the estimation engine: node-level, within-block two-way fixed-effects residualization inside honest causal trees, carried into the Callaway--Sant'Anna group-time structure. We provide Monte Carlo evidence that the estimator stays unbiased and correctly covered under timing heterogeneity, where two-way fixed effects and a pooled forest fail, and an application in which it recovers interpretable heterogeneity in the effect of Medicaid expansion on insurance coverage.
Some limitations remain and point to natural next steps. The consistency result (Proposition (ref)) is block-wise, and we establish the forbidden-comparison bias of the pooled alternative by simulation rather than in closed form. Inference for the conditional surface rests on a block bootstrap (Section (ref)); sharpening it to the uniform bands imai2023doubly derive for one covariate is a natural extension to the multivariate forest. The method is a complement to existing staggered-DiD tools: when the analyst wants to know not just how large the average effect is but for whom, and is willing to maintain parallel trends rather than unconfoundedness, it delivers the conditional object within the identification the literature already uses.
The staggered estimator is implemented in the causalfe Python package. The functionality this paper introduces is not yet part of a tagged public release; the replication archive bundles the exact package source used for these results and installs it locally, so the archive reproduces the paper without depending on a future release. The archive reproduces every number, table, and figure from committed inputs: each table and macro is generated by a single script (generate_paper_artifacts.py) from the Monte Carlo output and the application result files, and a one-command driver (reproduce.sh) either runs the full pipeline (the Monte Carlo including the MLDID benchmark, the callaway2021difference minimum-wage validation, and the Medicaid application with its scripted, key-free assembly of the SAHIE county panel from public Census files) or rebuilds the paper from committed results alone. All computation is seeded, so results are identical across runs and independent of the degree of parallelism. The estimator package is installed from the bundled source (pip install -e causalfe); other dependencies are pinned in requirements.txt; see REPLICATION.md for the full mapping from scripts to paper objects.