From f4235527cb01b678d8728231aa7336e8dac0fc19 Mon Sep 17 00:00:00 2001 From: igerber Date: Sat, 5 Sep 2026 08:05:23 -0400 Subject: [PATCH] fix(sdid): shape-only pre-treatment fit RMSE, level-gap field, placebo-anchored poor-fit warning SyntheticDiD's pre_treatment_fit measured the raw level residual while the Frank-Wolfe unit weights are fit on column-centered outcomes, so a parallel treated series at a different level reported a large RMSE and a false poor-fit warning. The RMSE is now taken on the pre-period residual after removing its mean (the quantity the solver minimised), computed on the normalized arrays and rescaled; the removed mean is exposed as pre_treatment_level_gap. The poor-fit warning is anchored to an in-space placebo fit reference (Abadie, Diamond & Hainmueller 2010; Abadie 2021): up to 20 Algorithm-4 pseudo-treated draws refit at the fit-time zeta from a private RNG stream, warning when the placebo p-value <= 0.05; new fields pre_fit_placebo_rmse / pre_fit_placebo_pvalue. in_time_placebo and sensitivity_to_zeta_omega report the same shape-only statistic; DR/BR prose and the NaN (single pre-period) case are handled. Estimates, SEs and weights are unchanged. Tutorials 03/18 updated (18 re-executed); a pre-existing geo-tutorial Frank-Wolfe non-convergence surfaced by the refresh is tracked in TODO.md. --- TODO.md | 1 + .../20260904-sdid-shape-only-fit-rmse.md | 44 +++ diff_diff/business_report.py | 4 +- diff_diff/diagnostic_report.py | 25 +- diff_diff/guides/llms-full.txt | 5 +- diff_diff/practitioner.py | 3 +- diff_diff/results.py | 85 ++++- diff_diff/synthetic_did.py | 226 +++++++++++-- .../diff_diff.SyntheticDiDResults.rst | 3 + docs/methodology/REGISTRY.md | 14 +- docs/methodology/REPORTING.md | 5 +- docs/references.rst | 2 + docs/tutorials/03_synthetic_did.ipynb | 34 +- docs/tutorials/18_geo_experiments.ipynb | 223 ++++++++----- tests/test_diagnostic_report.py | 39 +++ tests/test_estimators.py | 6 + tests/test_methodology_sdid.py | 301 +++++++++++++++++- 17 files changed, 872 insertions(+), 148 deletions(-) create mode 100644 changelog.d/20260904-sdid-shape-only-fit-rmse.md diff --git a/TODO.md b/TODO.md index ef6a4b10a..96e8d23b9 100644 --- a/TODO.md +++ b/TODO.md @@ -77,6 +77,7 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m | Issue | Location | Origin | Effort | Priority | |-------|----------|--------|--------|----------| +| Geo-experiments tutorial: the SDiD fits on the `generate_factor_data(...)` panel (conversions ~1,500, sigma-hat ~310) hit the Frank-Wolfe iteration cap - the placebo fit warns once and the bootstrap cross-check warns on 100/100 draws - yet the narrative presents the bootstrap SE as a paper-faithful robustness check. Pre-existing on `main` under BOTH backends (surfaced when the notebook outputs were refreshed for the shape-only fit diagnostic); the committed April outputs predate the aggregated non-convergence warning. Either rescale / re-tune the DGP (or raise `max_iter` / loosen `min_decrease`) so the refits converge and re-execute, or label the bootstrap result unreliable and drop it from the "Validity evidence" claim. | `docs/tutorials/18_geo_experiments.ipynb` | local codex review of the SDiD shape-only fit-RMSE PR | Mid | Low | | Replicate Chang (2020) §4.2.1 ML-design RCS DGP — requires a penalized (Lasso-class) propensity learner or a maintainer fixture via the duck-typed `propensity_learner=` object route; native logit+linear verified pure noise at (N=500, p=100) across ~20 review seeds (att ~ −23..+16, SE ~4-10, EPV ~1.4-1.6 vs threshold 10, fitted out-of-fold clipping ~5-29%, mean ~14%); at (N=200, p=100) `outcome_learner="linear"` fails closed on control-fold rank deficiency even under an oracle propensity (use ridge/sieve there); §4.2.2 is replicated | `tests/test_methodology_dml_did.py`, `diff_diff/_learners.py` | DML PR-B2 | Mid | Low | | Optional scheduled end-to-end execution gate for the MMM tutorials (29/30): a cron-only workflow (or extension of `mmm-interop.yml`) that executes both notebooks in isolated exact-pin environments, so a stale/invalid committed posterior cannot stay green indefinitely - today the hybrid posture (deliberate: notebooks execute locally with committed outputs; CI smoke-tests the exporters without sampling; drift tests pin source + committed-output needles) leaves the MCMC claims un-re-executed in CI | `.github/workflows/mmm-interop.yml`, `docs/tutorials/29_mmm_calibration_pymc.ipynb`, `docs/tutorials/30_mmm_calibration_meridian.ipynb` | mmm-interop | Mid | Low | | Committed `fixest::feols` event-study golden for TWFE `event_study=True` (within + pooled specs, unbalanced + covariate panels, matched CR1 cluster convention, per-period effects + vcov block) - the in-suite gates are shared-core cross-checks (TWFE-within == MPD-absorb, pooled == MPD bit-exact), so a defect common to the shared core would pass; the live-R harness (`benchmarks/R/benchmark_multiperiod.R`, `feols(y ~ treated * time_f \| unit)`) validated the within design in `docs/benchmarks.rst` but is not a committed regression test - follow the `fixest_did_twfe_golden.json` committed-golden pattern (pytest.skip when absent) | `tests/test_fixest_did_twfe_parity.py`, `benchmarks/R/` | 3(a) R2 | Mid | Medium | diff --git a/changelog.d/20260904-sdid-shape-only-fit-rmse.md b/changelog.d/20260904-sdid-shape-only-fit-rmse.md new file mode 100644 index 000000000..1f2e452d2 --- /dev/null +++ b/changelog.d/20260904-sdid-shape-only-fit-rmse.md @@ -0,0 +1,44 @@ +### Fixed +- **SyntheticDiD `pre_treatment_fit` is now shape-only**: the reported + pre-treatment RMSE (and the "Pre-treatment fit is poor" warning) previously + measured the raw level residual between the treated mean and the synthetic + control, while the Frank-Wolfe unit weights are fit on column-centered + outcomes (`intercept=True`, matching R `synthdid`) and deliberately leave a + constant level gap to the DiD step. A parallel treated series sitting at a + different level therefore reported a large RMSE and a false poor-fit + warning even when the ATT was recovered exactly. The RMSE is now taken on + the pre-period residual after removing its mean, which is the data-fit + component of the centered Frank-Wolfe objective, computed on the + normalized outcome scale and rescaled (so a large common outcome level + cannot perturb it); `in_time_placebo()` and `sensitivity_to_zeta_omega()` + report the same shape-only `pre_fit_rmse`. Results pickled before this + release are migrated on load: the shape RMSE and the level gap are + recomputed from the stored trajectories and replace the stale level RMSE + (which is cleared when no trajectories were stored), never relabeled. Estimates, standard errors + and weights are unchanged. + +### Behavioral Changes +- **SyntheticDiD pre-fit diagnostic redefinition**: `pre_treatment_fit` and + the `pre_fit_rmse` diagnostic columns drop the level gap, so their values + fall for any design with a treated-vs-synthetic level offset, and are NaN + with a single pre-period. The poor-fit warning is now anchored to an + in-space placebo fit reference (Abadie, Diamond & Hainmueller 2010; Abadie + 2021): the treated fit is compared with the same statistic over up to 20 + placebo fits of control units treated as if treated (Algorithm 4 draws at + the fit-time zeta), and the warning fires when the placebo p-value is at + or below 0.05 (at the 20-draw default: worse than every placebo draw), + replacing the unreachable `1 x std(treated pre-outcomes)` rule. The reference is computed for every variance method (one extra + Frank-Wolfe solve per draw) from a private RNG stream, so SE draws are + unchanged; it needs at least 19 successful draws (`n_bootstrap >= 19`) to + fire and is absent when no pseudo-control remains. A treated unit far + noisier than every control still fits worse than every placebo and warns; + the warning text says so. `summary()` labels the values + `Pre-fit RMSE (shape)`, `Pre-fit level gap` and `Pre-fit placebo p-value`. + +### Added +- **`SyntheticDiDResults.pre_treatment_level_gap`**: signed mean pre-period + gap (treated minus synthetic), the constant offset absorbed by the DiD + step, reported in `summary()` and `to_dict()` for inspection. +- **`SyntheticDiDResults.pre_fit_placebo_rmse` / `pre_fit_placebo_pvalue`**: + the placebo pre-fit reference distribution and the treated fit's placebo + p-value behind the poor-fit warning (p-value also in `to_dict()`). diff --git a/diff_diff/business_report.py b/diff_diff/business_report.py index 42eed4b9b..c08e4f6c8 100644 --- a/diff_diff/business_report.py +++ b/diff_diff/business_report.py @@ -2348,8 +2348,8 @@ def _render_summary(schema: Dict[str, Any]) -> str: else: sentences.append( "The synthetic control is designed to match the treated " - "group's pre-period trajectory (SDiD's weighted-parallel-" - "trends analogue)." + "group's pre-period trend (SDiD's weighted-parallel-" + "trends analogue; a constant level gap is differenced out)." ) elif verdict == "inconclusive": # Round-35 P1 CI review on PR #318: a ``verdict=="inconclusive"`` diff --git a/diff_diff/diagnostic_report.py b/diff_diff/diagnostic_report.py index 9f818cdd0..2be30b326 100644 --- a/diff_diff/diagnostic_report.py +++ b/diff_diff/diagnostic_report.py @@ -3647,8 +3647,16 @@ def _pt_synthetic_fit(self) -> Dict[str, Any]: SDiD's design-enforced fit quality substitutes for a standard PT test: the synthetic control is explicitly constructed to match the treated - group's pre-period trajectory, so small ``pre_treatment_fit`` RMSE + group's pre-period trend, so small ``pre_treatment_fit`` RMSE means the weighted-PT analogue is satisfied. + + ``pre_treatment_fit`` is SHAPE-ONLY (the constant pre-period level + gap is removed; SDID differences it out), whereas classic SCM's + ``pre_rmspe`` in ``_pt_scm_fit`` is level-inclusive by design. Both + are emitted under the shared ``pre_treatment_fit_rmse`` key as a + fit-quality RMSE in each estimator's own metric. NaN (a single + pre-period) is reported as ``skipped`` so the narrative never renders + ``RMSE = nan``. """ r = self._results fit = _to_python_float(getattr(r, "pre_treatment_fit", None)) @@ -3658,6 +3666,15 @@ def _pt_synthetic_fit(self) -> Dict[str, Any]: "method": "synthetic_fit", "reason": "SyntheticDiDResults.pre_treatment_fit is not populated " "on this fit.", } + if not np.isfinite(fit): + return { + "status": "skipped", + "method": "synthetic_fit", + "reason": ( + "SyntheticDiDResults.pre_treatment_fit is not defined on this " + "fit (fewer than 2 pre-periods)." + ), + } # Proxy verdict: unlike a classical PT p-value, this is a fit-quality # metric. Classify conservatively — phrasing in BR will explain that # this is SDiD's design-enforced analogue, not a PT hypothesis test. @@ -4576,10 +4593,10 @@ def _render_overall_interpretation(schema: Dict[str, Any], labels: Dict[str, str ) else: sentences.append( - f"The synthetic control matches the treated group's " - f"pre-period trajectory with RMSE = " + f"The synthetic control tracks the treated group's " + f"pre-period trend with shape-only RMSE = " f"{rmse:.3g} (SDiD's design-enforced analogue of parallel " - f"trends)." + f"trends; a constant level gap is differenced out by SDID)." if isinstance(rmse, (int, float)) else "SDiD's synthetic control is designed to satisfy the " "weighted parallel-trends analogue." diff --git a/diff_diff/guides/llms-full.txt b/diff_diff/guides/llms-full.txt index 34370378b..39e3cb9f5 100644 --- a/diff_diff/guides/llms-full.txt +++ b/diff_diff/guides/llms-full.txt @@ -1880,7 +1880,10 @@ Returned by `SyntheticDiD.fit()`. | `noise_level` | `float` | Estimated noise level | | `zeta_omega` | `float` | Unit weight regularization | | `zeta_lambda` | `float` | Time weight regularization | -| `pre_treatment_fit` | `float` | Pre-treatment RMSE | +| `pre_treatment_fit` | `float` | Shape-only pre-treatment RMSE (constant level gap removed: the data-fit component of the centered Frank-Wolfe objective; NaN with one pre-period) | +| `pre_treatment_level_gap` | `float` | Signed mean pre-period gap, treated minus synthetic (absorbed by the DiD step; diagnostic only) | +| `pre_fit_placebo_rmse` | `np.ndarray` | Reference distribution: shape-only pre-fit RMSE of up to 20 placebo (control-as-treated) fits; `None` when no placebo fit is possible | +| `pre_fit_placebo_pvalue` | `float` | Plus-one-adjusted rank of the treated fit among the placebo fits, `(1 + #{placebo >= treated}) / (1 + n_draws)`; a fit diagnostic, not a treatment-effect p-value; the poor-fit warning fires at <= 0.05 (Abadie, Diamond & Hainmueller 2010 in-space placebo fit assessment) | **Methods:** `summary()`, `print_summary()`, `to_dict()`, `to_dataframe()`, `get_unit_weights_df()`, `get_time_weights_df()` diff --git a/diff_diff/practitioner.py b/diff_diff/practitioner.py index 78b406986..15b4ca119 100644 --- a/diff_diff/practitioner.py +++ b/diff_diff/practitioner.py @@ -798,7 +798,8 @@ def _handle_synthetic(results: Any): "counterfactual well." ), code=( - "print(f'Pre-treatment fit (RMSE): {results.pre_treatment_fit:.4f}')\n" + "print(f'Pre-fit RMSE (shape): {results.pre_treatment_fit:.4f}')\n" + "print(f'Pre-fit level gap: {results.pre_treatment_level_gap:.4f}')\n" "concentration = results.get_weight_concentration()\n" "print(f\"Effective N: {concentration['effective_n']:.1f}\")\n" "print(f\"Top-5 weight share: {concentration['top_k_share']:.2%}\")" diff --git a/diff_diff/results.py b/diff_diff/results.py index 688bf688d..fa739f535 100644 --- a/diff_diff/results.py +++ b/diff_diff/results.py @@ -1149,6 +1149,31 @@ class SyntheticDiDResults(BaseResults): (for ``"jackknife"``). The ``variance_method`` field disambiguates the contents. (The deprecated read-only alias ``placebo_effects`` returns this array and is removed in v4.0.0.) + pre_treatment_fit : float, optional + Shape-only pre-treatment fit: the RMSE of the pre-period residual + (treated mean minus synthetic control) after removing its mean. The + Frank-Wolfe unit weights are fit on column-centered outcomes, so a + constant level gap is not a fit failure (SDID differences it out) + and is excluded here. NaN with a single pre-period. + pre_treatment_level_gap : float, optional + Signed mean pre-period gap, treated minus synthetic control. This is + the constant offset absorbed by the DiD step; it is reported for + inspection and never enters the poor-fit warning. + pre_fit_placebo_rmse : np.ndarray, optional + Reference distribution for ``pre_treatment_fit``: the shape-only + pre-fit RMSE of each placebo fit in which a random set of + ``n_treated`` control units is treated as if treated and the unit + weights are re-estimated on the remaining controls (Algorithm 4 + draws, fit-time zeta). At most 20 draws (``min(n_bootstrap, 20)``); + ``None`` when no placebo fit is possible (no pseudo-control left, or + a single pre-period). + pre_fit_placebo_pvalue : float, optional + ``(1 + #{placebo RMSE >= treated RMSE}) / (1 + n_draws)``: the + plus-one-adjusted rank of the treated units' fit among the placebo + fits — a fit diagnostic, not a treatment-effect p-value. The + poor-fit warning fires at ``<= 0.05`` (in-space placebo fit + assessment, Abadie, Diamond & Hainmueller 2010; Abadie 2021), which + needs at least 19 successful draws to be reachable. synthetic_pre_trajectory : np.ndarray, optional Synthetic control trajectory in pre-treatment periods, shape ``(n_pre,)``. Equal to ``Y_pre_control @ omega_eff`` where @@ -1186,6 +1211,9 @@ class SyntheticDiDResults(BaseResults): zeta_omega: Optional[float] = field(default=None) zeta_lambda: Optional[float] = field(default=None) pre_treatment_fit: Optional[float] = field(default=None) + pre_treatment_level_gap: Optional[float] = field(default=None) + pre_fit_placebo_rmse: Optional[np.ndarray] = field(default=None) + pre_fit_placebo_pvalue: Optional[float] = field(default=None) variance_effects: Optional[np.ndarray] = field(default=None) n_bootstrap: Optional[int] = field(default=None) # Survey design metadata (SurveyMetadata instance from diff_diff.survey) @@ -1256,6 +1284,33 @@ def __setstate__(self, state: Dict[str, Any]) -> None: if "placebo_effects" in state and "variance_effects" not in state: state = dict(state) state["variance_effects"] = state.pop("placebo_effects") + # Pre-v3.11.2 pickles: ``pre_treatment_fit`` was the LEVEL-inclusive + # RMSE and the shape-only / placebo-reference fields did not exist. + # Never relabel the stale value as shape-only: recompute both + # statistics from the stored trajectories when they are present + # (retained on results since v3.8), otherwise clear it. The placebo + # reference cannot be rebuilt without the panel, so it stays None + # (the poor-fit warning is a fit-time event and is not replayed). + if "pre_treatment_level_gap" not in state: + state = dict(state) + treated_pre = state.get("treated_pre_trajectory") + synthetic_pre = state.get("synthetic_pre_trajectory") + if treated_pre is not None and synthetic_pre is not None: + resid = np.asarray(treated_pre, dtype=np.float64) - np.asarray( + synthetic_pre, dtype=np.float64 + ) + level_gap = float(np.mean(resid)) + state["pre_treatment_level_gap"] = level_gap + state["pre_treatment_fit"] = ( + float(np.sqrt(np.mean((resid - level_gap) ** 2))) + if resid.shape[0] >= 2 + else float("nan") + ) + else: + state["pre_treatment_level_gap"] = None + state["pre_treatment_fit"] = None + state.setdefault("pre_fit_placebo_rmse", None) + state.setdefault("pre_fit_placebo_pvalue", None) self.__dict__.update(state) @property @@ -1326,7 +1381,11 @@ def summary(self, alpha: Optional[float] = None) -> str: lines.append(f"{'Noise level:':<25} {self.noise_level:>10.4f}") if self.pre_treatment_fit is not None: - lines.append(f"{'Pre-treatment fit (RMSE):':<25} {self.pre_treatment_fit:>10.4f}") + lines.append(f"{'Pre-fit RMSE (shape):':<25} {self.pre_treatment_fit:>10.4f}") + if self.pre_treatment_level_gap is not None: + lines.append(f"{'Pre-fit level gap:':<25} {self.pre_treatment_level_gap:>10.4f}") + if self.pre_fit_placebo_pvalue is not None: + lines.append(f"{'Pre-fit placebo p-value:':<25} {self.pre_fit_placebo_pvalue:>10.3f}") # Variance method info lines.append(f"{'Variance method:':<25} {self.variance_method:>10}") @@ -1416,6 +1475,8 @@ def to_dict(self) -> Dict[str, Any]: "zeta_omega": self.zeta_omega, "zeta_lambda": self.zeta_lambda, "pre_treatment_fit": self.pre_treatment_fit, + "pre_treatment_level_gap": self.pre_treatment_level_gap, + "pre_fit_placebo_pvalue": self.pre_fit_placebo_pvalue, } if self.n_bootstrap is not None: result["n_bootstrap"] = self.n_bootstrap @@ -1613,7 +1674,9 @@ def in_time_placebo( Columns: - ``fake_treatment_period`` — the shifted date - ``att`` — placebo ATT (ideally near 0) - - ``pre_fit_rmse`` — RMSE on the fake pre-window + - ``pre_fit_rmse`` — shape-only RMSE on the fake pre-window + (pre-period mean gap removed, matching + ``pre_treatment_fit``) - ``n_pre_fake`` — periods before the fake date - ``n_post_fake`` — periods from the fake date onward @@ -1749,7 +1812,10 @@ def in_time_placebo( lambda_fake, ) synthetic_pre_fake_n = Y_pre_c_n @ omega_eff_fake - pre_fit_n = float(np.sqrt(np.mean((y_pre_t_mean_n - synthetic_pre_fake_n) ** 2))) + # Shape-only RMSE (pre-window mean gap removed), matching the + # fit-time ``pre_treatment_fit`` definition. n_pre_fake >= 2 here. + resid_fake_n = y_pre_t_mean_n - synthetic_pre_fake_n + pre_fit_n = float(np.sqrt(np.mean((resid_fake_n - resid_fake_n.mean()) ** 2))) # ATT is scale-equivariant and shift-invariant in Y; RMSE is # scale-equivariant. Rescale back to original-Y units. row["att"] = float(att_fake_n * Y_scale) @@ -1789,7 +1855,9 @@ def sensitivity_to_zeta_omega( Columns: - ``zeta_omega`` — the regularization value evaluated - ``att`` — resulting ATT - - ``pre_fit_rmse`` — RMSE on the original pre-period + - ``pre_fit_rmse`` — shape-only RMSE on the original + pre-period (pre-period mean gap removed, matching + ``pre_treatment_fit``; NaN with a single pre-period) - ``max_unit_weight`` — max element of the composed ``omega_eff`` (sensitivity indicator: close to 1 means near-one-hot solutions; close to ``1/n_control`` means @@ -1895,7 +1963,14 @@ def sensitivity_to_zeta_omega( time_weights, ) synthetic_pre_n = Y_pre_control_n @ omega_eff - pre_fit_n = float(np.sqrt(np.mean((y_pre_t_mean_n - synthetic_pre_n) ** 2))) + # Shape-only RMSE (pre-period mean gap removed), matching the + # fit-time ``pre_treatment_fit`` definition; undefined (NaN) with + # a single pre-period, which is a legal fit. + resid_n = y_pre_t_mean_n - synthetic_pre_n + if resid_n.shape[0] >= 2: + pre_fit_n = float(np.sqrt(np.mean((resid_n - resid_n.mean()) ** 2))) + else: + pre_fit_n = float("nan") herf = float(np.sum(omega_eff**2)) rows.append( { diff --git a/diff_diff/synthetic_did.py b/diff_diff/synthetic_did.py index 56146c7e4..a04c01f56 100644 --- a/diff_diff/synthetic_did.py +++ b/diff_diff/synthetic_did.py @@ -26,6 +26,24 @@ validate_n_bootstrap, ) +# Poor-fit warning anchor: the treated units' shape-only pre-treatment RMSE +# is compared against the distribution of the same statistic over placebo +# (control-as-treated) fits — the in-space placebo fit assessment of Abadie, +# Diamond & Hainmueller (2010) / Abadie (2021), transposed onto SDID's +# Algorithm 4 pseudo-treated draws. ``_PRE_FIT_REFERENCE_DRAWS`` caps the +# number of placebo pre-fits (each is one Frank-Wolfe solve; under +# ``variance_method="placebo"`` they are the first draws of the SE loop's +# permutation stream); the warning fires when the placebo p-value +# ``(1 + #{placebo_rmse >= treated_rmse}) / (1 + n_draws)`` is at or below +# ``_PRE_FIT_PLACEBO_ALPHA``, so it needs at least 19 successful draws +# (n_bootstrap >= 19) to be reachable at all. The cap of 20 keeps the cost +# at ~20 Frank-Wolfe solves per fit (the pure-Python test suite doubled in +# wall-clock at 50) while leaving the rule reachable: at 20 draws it fires +# exactly when the treated fit is worse than every placebo draw +# (p = 1/21 ~ 0.048). +_PRE_FIT_REFERENCE_DRAWS = 20 +_PRE_FIT_PLACEBO_ALPHA = 0.05 + class SyntheticDiD(DifferenceInDifferences): """ @@ -82,7 +100,11 @@ class SyntheticDiD(DifferenceInDifferences): Number of replications for variance estimation. Used for: - Bootstrap: Number of bootstrap samples - Placebo: Number of random permutations (matches R's `replications` argument) - Ignored when ``variance_method="jackknife"``. + Ignored by jackknife variance estimation. For every variance method + it also caps the pre-treatment fit reference distribution + (``min(n_bootstrap, 20)`` placebo pre-fit draws behind + ``pre_fit_placebo_rmse`` / ``pre_fit_placebo_pvalue`` and the + poor-fit warning, which needs at least 19 draws to fire). seed : int, optional Random seed for reproducibility. If None (default), results will vary between runs. @@ -827,31 +849,38 @@ def fit( # type: ignore[override] Y_pre_treated_mean = Y_pre_treated_mean_n * Y_scale + Y_shift Y_post_treated_mean = Y_post_treated_mean_n * Y_scale + Y_shift - # Compute pre-treatment fit (RMSE) using composed weights on the + # Pre-treatment fit diagnostics using composed weights on the # original Y (user-visible scale). omega_eff is a simplex — applies # cleanly to any linear rescale of Y — so trajectories live on the # original outcome scale for plotting and the poor-fit warning. + # + # The fit statistic is SHAPE-ONLY: the Frank-Wolfe unit-weight + # objective column-centers the pre-period matrix (intercept=True, + # matching R synthdid), so omega is chosen to match the treated + # pre-period *movements* and a constant level gap is deliberately + # left to the DiD step. The RMSE is therefore taken on the residual + # after removing its pre-period mean — on the non-survey path this + # is exactly the residual the solver minimised; on the survey path + # omega_eff is the post-hoc composed/renormalised vector, so it is the + # fit of the weights actually used. The removed mean is reported + # separately as the signed level gap (treated minus synthetic). synthetic_pre_trajectory = Y_pre_control @ omega_eff synthetic_post_trajectory = Y_post_control @ omega_eff - pre_fit_rmse = np.sqrt(np.mean((Y_pre_treated_mean - synthetic_pre_trajectory) ** 2)) - - # Warn if pre-treatment fit is poor (Registry requirement). - # Threshold: 1× SD of treated pre-treatment outcomes — a natural baseline - # since RMSE exceeding natural variation indicates the synthetic control - # fails to reproduce the treated series' level or trend. - pre_treatment_sd = ( - np.std(Y_pre_treated_mean, ddof=1) if len(Y_pre_treated_mean) > 1 else 0.0 - ) - if pre_treatment_sd > 0 and pre_fit_rmse > pre_treatment_sd: - warnings.warn( - f"Pre-treatment fit is poor: RMSE ({pre_fit_rmse:.4f}) exceeds " - f"the standard deviation of treated pre-treatment outcomes " - f"({pre_treatment_sd:.4f}). The synthetic control may not " - f"adequately reproduce treated unit trends. Consider adding " - f"more control units or adjusting regularization.", - UserWarning, - stacklevel=2, - ) + # The two scalar diagnostics are computed on the NORMALIZED arrays + # (same contract as the estimator and the post-fit diagnostics in + # results.py) and rescaled by Y_scale, so a large common outcome + # level cannot perturb them through floating-point cancellation. + pre_resid_n = Y_pre_treated_mean_n - Y_pre_control_n @ omega_eff + pre_level_gap_n = float(np.mean(pre_resid_n)) + if len(pre_resid_n) >= 2: + pre_fit_rmse_n = float(np.sqrt(np.mean((pre_resid_n - pre_level_gap_n) ** 2))) + else: + # Shape is undefined with a single pre-period. + pre_fit_rmse_n = float("nan") + pre_fit_rmse = pre_fit_rmse_n * Y_scale + pre_level_gap = pre_level_gap_n * Y_scale + # The poor-fit warning is emitted after the variance block: its + # reference distribution (placebo pre-fits) is computed there. # Treated-unit trajectories (the pre/post means already computed above). treated_pre_trajectory = Y_pre_treated_mean @@ -1204,6 +1233,55 @@ def fit( # type: ignore[override] variance_effects = np.asarray(variance_effects_n) * Y_scale inference_method = "placebo" + # Pre-treatment fit reference distribution (Registry requirement): + # the treated units' shape-only pre-fit RMSE is placed within the + # distribution of the same statistic over placebo (control-as- + # treated) fits — Algorithm 4's pseudo-treated draws, refit with the + # fit-time zeta — the in-space placebo fit assessment of Abadie, + # Diamond & Hainmueller (2010) / Abadie (2021). Independent of the + # variance method (its own generator seeded from ``self.seed``, so + # the SE loops' RNG streams are untouched; with a fixed seed on the + # non-survey / pweight-only placebo path the draws coincide with the + # SE loop's first permutations — not on full-design surveys, whose SE + # loop permutes within strata, and not when ``seed`` is None). + pre_fit_placebo_rmse: Optional[np.ndarray] = None + pre_fit_placebo_pvalue: Optional[float] = None + if np.isfinite(pre_fit_rmse_n): + pre_fit_placebo_rmse_n = self._placebo_pre_fit_reference( + Y_pre_control_n, + n_treated=len(treated_units), + zeta_omega=zeta_omega_n, + min_decrease=min_decrease, + w_control=w_control, + init_omega=unit_weights, + n_draws=min(self.n_bootstrap, _PRE_FIT_REFERENCE_DRAWS), + ) + if len(pre_fit_placebo_rmse_n) > 0: + # Placebo p-value (compared on the normalized scale; Y_scale + # is a positive constant so the ordering is unchanged). + n_ge = int(np.sum(pre_fit_placebo_rmse_n >= pre_fit_rmse_n)) + pre_fit_placebo_pvalue = (1.0 + n_ge) / (1.0 + len(pre_fit_placebo_rmse_n)) + pre_fit_placebo_rmse = pre_fit_placebo_rmse_n * Y_scale + pre_fit_placebo_rmse.flags.writeable = False + if pre_fit_placebo_pvalue <= _PRE_FIT_PLACEBO_ALPHA: + warnings.warn( + f"Pre-treatment fit is poor: the treated units' shape-only " + f"pre-fit RMSE ({pre_fit_rmse:.4f}) is worse than " + f"{100 * (1 - pre_fit_placebo_pvalue):.0f}% of " + f"{len(pre_fit_placebo_rmse_n)} placebo fits of control units " + f"treated as if treated (placebo p-value " + f"{pre_fit_placebo_pvalue:.3f}; median placebo RMSE " + f"{np.median(pre_fit_placebo_rmse):.4f}). The synthetic control " + f"may not reproduce the treated units' pre-treatment trend (a " + f"constant level gap is excluded - SDID differences it out; see " + f"pre_treatment_level_gap). This can also occur when treated " + f"units are much noisier than the controls: inspect " + f"treated_pre_trajectory against synthetic_pre_trajectory. " + f"Consider adding more control units or adjusting regularization.", + UserWarning, + stacklevel=2, + ) + # Compute test statistics t_stat, p_value_analytical, conf_int = safe_inference(att, se, alpha=self.alpha) # Empirical p-value is valid only for placebo (Algorithm 4): control @@ -1288,6 +1366,9 @@ def fit( # type: ignore[override] zeta_omega=zeta_omega, zeta_lambda=zeta_lambda, pre_treatment_fit=pre_fit_rmse, + pre_treatment_level_gap=pre_level_gap, + pre_fit_placebo_rmse=pre_fit_placebo_rmse, + pre_fit_placebo_pvalue=pre_fit_placebo_pvalue, variance_effects=variance_effects if len(variance_effects) > 0 else None, n_bootstrap=self.n_bootstrap if inference_method == "bootstrap" else None, survey_metadata=survey_metadata, @@ -2006,6 +2087,109 @@ def _placebo_variance_se( return se, placebo_estimates + def _placebo_pre_fit_reference( + self, + Y_pre_control: np.ndarray, + n_treated: int, + zeta_omega: float, + min_decrease: float, + w_control: Optional[np.ndarray], + init_omega: Optional[np.ndarray], + n_draws: int, + ) -> np.ndarray: + """Shape-only pre-fit RMSEs of placebo (control-as-treated) fits. + + Runs steps 1-3 of Algorithm 4 (Arkhangelsky et al. 2021): permute the + control indices, designate the last ``n_treated`` as pseudo-treated, + re-estimate the unit weights on the pseudo-controls with the fit-time + ``zeta_omega`` (warm-started from ``init_omega`` exactly as the + placebo SE loop is), and return the shape-only RMSE of each + pseudo-treated pre-period residual (mean removed) — the same + statistic ``pre_treatment_fit`` reports for the treated units. The + treated units' RMSE is then placed within this distribution, the + in-space placebo fit assessment of Abadie, Diamond & Hainmueller + (2010) / Abadie (2021) transposed onto SDID. + + Diagnostic only: per-draw Frank-Wolfe non-convergence warnings are + suppressed, failed draws are skipped, and the permutation stream + (``np.random.default_rng(self.seed)``) is private to this loop. Uses + the unstratified permutation for every survey design (pweight-only + composition ``omega * w_control`` mirrors the placebo SE loop); + returns an empty array when there is no pseudo-control left or fewer + than 2 pre-periods. + + Parameters + ---------- + Y_pre_control : np.ndarray + Normalized control pre-outcomes, shape ``(n_pre, n_control)``. + n_treated : int + Number of treated units (size of each pseudo-treated set). + zeta_omega : float + Fit-time unit-weight regularization on the normalized scale. + min_decrease : float + Frank-Wolfe convergence threshold used at fit time. + w_control : np.ndarray, optional + Per-control survey weights (pweight-only composition). + init_omega : np.ndarray, optional + Fit-time raw Frank-Wolfe unit weights used as warm-start. + n_draws : int + Number of placebo draws. + + Returns + ------- + np.ndarray + Shape-only RMSEs (normalized scale) of the successful draws. + """ + n_pre, n_control = Y_pre_control.shape + n_pseudo_control = n_control - n_treated + if n_pseudo_control < 1 or n_pre < 2 or n_draws < 1: + return np.array([]) + rng = np.random.default_rng(self.seed) + out: List[float] = [] + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + for _rep in range(n_draws): + try: + perm = rng.permutation(n_control) + pseudo_control_idx = perm[:n_pseudo_control] + pseudo_treated_idx = perm[n_pseudo_control:] + Y_pre_pseudo_control = Y_pre_control[:, pseudo_control_idx] + if w_control is not None: + y_pre_pseudo_treated = np.average( + Y_pre_control[:, pseudo_treated_idx], + axis=1, + weights=w_control[pseudo_treated_idx], + ) + else: + y_pre_pseudo_treated = np.mean(Y_pre_control[:, pseudo_treated_idx], axis=1) + pseudo_omega_init = ( + _sum_normalize(init_omega[pseudo_control_idx]) + if init_omega is not None + else None + ) + pseudo_omega = compute_sdid_unit_weights( + Y_pre_pseudo_control, + y_pre_pseudo_treated, + zeta_omega=zeta_omega, + min_decrease=min_decrease, + init_weights=pseudo_omega_init, + ) + if w_control is not None: + pseudo_omega_eff = pseudo_omega * w_control[pseudo_control_idx] + denom = float(pseudo_omega_eff.sum()) + if denom <= 0: + continue + pseudo_omega_eff = pseudo_omega_eff / denom + else: + pseudo_omega_eff = pseudo_omega + resid = y_pre_pseudo_treated - Y_pre_pseudo_control @ pseudo_omega_eff + rmse = float(np.sqrt(np.mean((resid - resid.mean()) ** 2))) + if np.isfinite(rmse): + out.append(rmse) + except (ValueError, LinAlgError, ZeroDivisionError): + continue + return np.asarray(out, dtype=np.float64) + def _placebo_variance_se_survey( self, Y_pre_control: np.ndarray, diff --git a/docs/api/_autosummary/diff_diff.SyntheticDiDResults.rst b/docs/api/_autosummary/diff_diff.SyntheticDiDResults.rst index 16b1d35d4..1e6aa81b4 100644 --- a/docs/api/_autosummary/diff_diff.SyntheticDiDResults.rst +++ b/docs/api/_autosummary/diff_diff.SyntheticDiDResults.rst @@ -36,7 +36,10 @@ ~SyntheticDiDResults.n_bootstrap ~SyntheticDiDResults.noise_level ~SyntheticDiDResults.placebo_effects + ~SyntheticDiDResults.pre_fit_placebo_pvalue + ~SyntheticDiDResults.pre_fit_placebo_rmse ~SyntheticDiDResults.pre_treatment_fit + ~SyntheticDiDResults.pre_treatment_level_gap ~SyntheticDiDResults.significance_stars ~SyntheticDiDResults.survey_metadata ~SyntheticDiDResults.synthetic_post_trajectory diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 7e700e3c1..e19c8601c 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -3235,7 +3235,7 @@ the finite-dimensional `p_0` is handled by the variance correction below. *Assumption checks / warnings:* - Requires balanced panel (same units observed in all periods) -- Warns if pre-treatment fit is poor (high RMSE) +- Warns if pre-treatment fit is poor (shape-only RMSE worse than ≥95% of placebo control-as-treated fits; see the fit-diagnostic Note below) - Treatment must be "block" structure: all treated units treated at same time *Estimator equation (as implemented):* @@ -3263,6 +3263,8 @@ where A = Y_unit[:, :N_co], b = Y_unit[:, N_co], and centering is column-wise (i The sparsification step concentrates weights on the most important control units, improving interpretability and stability. +- **Note:** Pre-treatment fit diagnostic (`pre_treatment_fit`, `pre_treatment_level_gap`, the poor-fit warning; v3.11.2). The fit statistic is measured on the same centered residual the `intercept=True` Frank-Wolfe objective minimises: `pre_treatment_fit = RMSE(resid - mean(resid))` with `resid = treated_pre_trajectory - synthetic_pre_trajectory`, and `pre_treatment_level_gap = mean(resid)` (signed, treated minus synthetic). A constant treated-vs-synthetic level gap is permitted by SDID — it is absorbed by the DiD step, which is exactly why the unit weights are fit on column-centered outcomes — so it is reported separately and never counted as poor fit. On the non-survey path this is exactly the residual the solver minimised; on the survey path `ω_eff` is the post-hoc composed/renormalised vector (see the survey Notes below), so it is the fit of the weights actually used. R has no user-facing fit statistic (`synthdid_rmse_plot` plots `sqrt(vals)`, the centered Frank-Wolfe objective, and `synthdid_plot` shifts the synthetic trajectory by an intercept offset), so this is consistent with R. HISTORY: before v3.11.2 the RMSE was taken on raw levels, so a parallel treated series at a different level reported a large RMSE and a false "poor fit" warning even when the ATT was recovered exactly. Both scalars are computed on the normalized arrays (the Y-normalization contract below) and rescaled by `Y_scale`, so a large common outcome level cannot perturb them through floating-point cancellation, and the fit-time value equals the multiplier-1.0 row of `sensitivity_to_zeta_omega()` exactly. **Warning threshold — in-space placebo fit reference** (Abadie, Diamond & Hainmueller 2010, Section on placebo tests; Abadie 2021, JEL): the treated units' `pre_treatment_fit` is placed within the distribution of the same statistic over placebo fits in which a random set of `n_treated` control units is treated as if treated and ω is re-estimated on the remaining controls at the fit-time ζ (steps 1–3 of Algorithm 4, warm-started like the placebo SE loop, `n_draws = min(n_bootstrap, 20)`; `results.pre_fit_placebo_rmse` holds the draws on the outcome scale). The warning fires when `pre_fit_placebo_pvalue = (1 + #{placebo_rmse ≥ treated_rmse}) / (1 + n_draws) ≤ 0.05`, i.e. the treated fit is worse than at least 95% of the placebo fits — the SDID transposition of the SCM practice of judging a treated unit's pre-RMSPE against the in-space placebo pre-RMSPE distribution. The 20-draw cap is a cost decision (each draw is one Frank-Wolfe solve; the pure-Python test suite doubled in wall-clock at 50) chosen just above the minimum of 19 draws at which α = 0.05 is reachable (`1/(19+1) = 0.05`), so at the default the rule reduces to "worse than every placebo draw" (p = 1/21 ≈ 0.048); `n_bootstrap < 19` silences it. The reference loop owns a private permutation stream seeded from `seed` (the SE loops' RNG streams are untouched; with a fixed `seed` on the non-survey / pweight-only `variance_method="placebo"` path the draws coincide with the SE loop's first permutations — not on full-design surveys, whose SE loop permutes within strata, and not when `seed=None`), runs for every variance method, suppresses per-draw Frank-Wolfe non-convergence warnings, skips failed draws, and uses the unstratified permutation for every survey design (pweight-only composition `ω · w_control` as in the placebo SE loop). Calibration (2 treated, 8 controls, T=8, iid noise 0.5, 20 draws): adequately fit designs (parallel trends with a level offset, flat series at different levels, a 1e9 common offset, T_pre = 4) give p ≈ 0.48–0.67; misfit designs (treated trends vs flat controls, quadratic or sinusoidal treated vs flat controls, opposite trends, mild trend vs flat, noiseless flat controls vs a trending treated series) all give p = 1/21 ≈ 0.048 (worse than every draw). **Known limitation** (inherent to any control-side reference, accepted and documented): treated units much noisier than the controls (e.g. 50x) fit worse than every placebo and trigger the warning on an adequate fit; the warning text says so and points to the trajectories. Dispersed control trends inflate the placebo fits too, so a treated misfit of comparable size is borderline (a quadratic treated series against controls with slopes ~U(-3, 3): p = 0.048 at 20 draws, 0.06 at 50). HISTORY of rejected anchors: the pre-v3.11.2 `1 x std(treated_pre, ddof=1)` rule is unreachable for a shape-only RMSE (a constant synthetic against a trending treated series gives exactly the population SD, below the ddof=1 SD) and fires at random when the treated series is flat noise; a `2 x noise_level` (σ̂) rule was calibrated (good 0.2–0.5, misfit 2.9–14.9) but σ̂ is control-side and pools first differences across units, so it needs the same caveats without the paper anchor; a residual-whiteness anchor (`2 * sd(diff(resid)) / sqrt(2)`) is robust to noisy treated units but blind to oscillatory misfit in short windows. `pre_treatment_fit` is NaN with a single pre-period (shape undefined; `pre_treatment_level_gap` is still defined; no reference distribution), and there is no reference distribution when no pseudo-control remains (`n_control ≤ n_treated`). **Legacy pickles**: `SyntheticDiDResults.__setstate__` detects a pre-v3.11.2 state (no `pre_treatment_level_gap` key), recomputes the shape-only RMSE and the level gap from the stored trajectories (or clears the stale level RMSE to `None` when no trajectories were stored), and defaults the placebo-reference fields to `None` — the stale level RMSE is never relabeled as shape-only. + *Time weights λ (Frank-Wolfe on collapsed form):* Build collapsed-form matrix Y_time of shape (N_co, T_pre + 1), where the last column is the per-control post-period mean (averaged across post-periods for each control unit). Solve: @@ -3356,7 +3358,7 @@ Convergence criterion: stop when objective decrease < min_decrease² (default mi - **Placebo p-value floor**: `p_value = max(empirical_p, 1/(n_replications + 1))` to avoid reporting exactly zero. - **Varying treatment within unit**: Raises `ValueError`. SDID requires block treatment (constant within each unit). Suggests CallawaySantAnna or ImputationDiD for staggered adoption. - **Unbalanced panel**: Raises `ValueError`. SDID requires all units observed in all periods. Suggests `balance_panel()`. -- **Poor pre-treatment fit**: Warns (`UserWarning`) when `pre_fit_rmse > std(treated_pre_outcomes, ddof=1)`. Diagnostic only; estimation proceeds. +- **Poor pre-treatment fit**: Warns (`UserWarning`) when the treated units' shape-only `pre_treatment_fit` (RMSE of the pre-period residual after removing its mean; NaN when `n_pre < 2`) is at least as poor as ≥95% of placebo (control-as-treated) fits: `pre_fit_placebo_pvalue = (1 + #{placebo_rmse ≥ treated_rmse}) / (1 + n_draws) ≤ 0.05`, with `n_draws = min(n_bootstrap, 20)` Algorithm-4 pseudo-treated draws refit at the fit-time ζ (at the default 20 draws this fires exactly when the treated fit is worse than every draw, p = 1/21 ≈ 0.048). `pre_treatment_level_gap = mean(resid)` is reported alongside and never enters the rule. No reference distribution (and no warning) when no pseudo-control remains (`n_control ≤ n_treated`) or `n_pre < 2`; the rule needs ≥19 successful draws to be reachable (`n_bootstrap < 19` silences it). Diagnostic only; estimation proceeds. See the fit-diagnostic Note under *Unit weights ω*. - **Jackknife with single treated unit**: Returns NaN SE. Cannot leave-one-out with N_tr=1; R returns NA for the same condition. - **Jackknife with single nonzero-weight control**: Returns NaN SE. Leaving out the only effective control is not meaningful. - **Jackknife with non-finite LOO estimate**: Returns NaN SE. Unlike bootstrap/placebo, jackknife is deterministic and cannot skip failed iterations; NaN propagates through `var()` (matches R behavior). @@ -3490,17 +3492,17 @@ Convergence criterion: stop when objective decrease < min_decrease² (default mi The schema smoke test is `TestCoverageMCArtifact::test_coverage_artifacts_present`; regenerate the JSON via `python benchmarks/python/coverage_sdid.py --n-seeds 500 --n-bootstrap 200 --output benchmarks/data/sdid_coverage.json` (~15–40 min on M-series Mac, Rust backend — warm-start convergence makes newer runs faster than the original cold-start one). **Artifact cadence (documentation-grade, not pinned-regression-grade):** this file is a documentation substrate — the table above transcribes from it, and the schema test just checks structure. It is **not** numerically pinned by any regression test, because MC noise at 500 seeds × B=200 (2σ ≈ 0.02–0.05 per cell) makes tight bounds fragile and loose bounds uninformative. The runtime characterization test that guards the one regression class worth catching (dispatch breakage → rejection rate ≈0 or ≈0.5, or SE-collapse) is `TestPValueSemantics::test_bootstrap_p_value_null_dispersion` (slow; calibration-agnostic dispersion + loose rejection-rate band). Regenerate the JSON when a methodology change materially shifts per-draw numerics — SE formula, new variance method, FW solver swap, paper-procedure correction. **Do not** regenerate on refactors that preserve the FW global optimum (warm-start vs cold-start, backend migration, pure renames, docstring fixes) — those stay within MC noise of the committed numbers. Per-seed bit-identity on the captured fixture is the cheaper, stricter check: see `TestScaleEquivariance::test_baseline_parity_small_scale[bootstrap]` at `rel=1e-14`. -- **Note:** Internal Y normalization. Before weight optimization, the estimator, and variance procedures, `fit()` centers Y by `mean(Y_pre_control)` and scales by `std(Y_pre_control)`; `Y_scale` falls back to `1.0` when std is non-finite or below `1e-12 * max(|mean|, 1)`. Auto-regularization and `noise_level` are computed on normalized Y; user-supplied `zeta_omega` / `zeta_lambda` are divided by `Y_scale` internally for Frank-Wolfe. τ, SE, CI, the placebo/bootstrap/jackknife effect vectors, `results_.noise_level`, and `results_.zeta_omega` / `results_.zeta_lambda` are all reported on the user's original outcome scale (user-supplied zetas are echoed back exactly to avoid float roundoff). Mathematically a no-op — τ is location-invariant and scale-equivariant, and FW weights are invariant under `(Y, ζ) → (Y/s, ζ/s)` — but prevents catastrophic cancellation in the SDID double-difference when outcomes span millions-to-billions (see synth-inference/synthdid#71 for the R-package version of this issue). Normalization constants are derived from controls' pre-period only so the reference is unaffected by treatment. `in_time_placebo()` and `sensitivity_to_zeta_omega()` reuse the exact same `Y_shift` / `Y_scale` captured on the fit snapshot: they normalize the re-sliced arrays before re-running Frank-Wolfe, pass `zeta / Y_scale` to the weight solvers, and rescale the returned `att` and `pre_fit_rmse` by `Y_scale` before reporting; unit-weight diagnostics (`max_unit_weight`, `effective_n`) are scale-invariant and reported directly. +- **Note:** Internal Y normalization. Before weight optimization, the estimator, and variance procedures, `fit()` centers Y by `mean(Y_pre_control)` and scales by `std(Y_pre_control)`; `Y_scale` falls back to `1.0` when std is non-finite or below `1e-12 * max(|mean|, 1)`. Auto-regularization and `noise_level` are computed on normalized Y; user-supplied `zeta_omega` / `zeta_lambda` are divided by `Y_scale` internally for Frank-Wolfe. τ, SE, CI, the placebo/bootstrap/jackknife effect vectors, `results_.noise_level`, and `results_.zeta_omega` / `results_.zeta_lambda` are all reported on the user's original outcome scale (user-supplied zetas are echoed back exactly to avoid float roundoff). Mathematically a no-op — τ is location-invariant and scale-equivariant, and FW weights are invariant under `(Y, ζ) → (Y/s, ζ/s)` — but prevents catastrophic cancellation in the SDID double-difference when outcomes span millions-to-billions (see synth-inference/synthdid#71 for the R-package version of this issue). Normalization constants are derived from controls' pre-period only so the reference is unaffected by treatment. `in_time_placebo()` and `sensitivity_to_zeta_omega()` reuse the exact same `Y_shift` / `Y_scale` captured on the fit snapshot: they normalize the re-sliced arrays before re-running Frank-Wolfe, pass `zeta / Y_scale` to the weight solvers, and rescale the returned `att` and (shape-only) `pre_fit_rmse` by `Y_scale` before reporting; unit-weight diagnostics (`max_unit_weight`, `effective_n`) are scale-invariant and reported directly. *Validation diagnostics (post-fit methods on `SyntheticDiDResults`):* -- **Trajectories** (`synthetic_pre_trajectory`, `synthetic_post_trajectory`, `treated_pre_trajectory`, `treated_post_trajectory`): retained on results to support plotting and custom fit metrics. `synthetic_pre_trajectory = Y_pre_control @ ω_eff`; `treated_pre_trajectory` is the survey-weighted treated mean (matches the Frank-Wolfe target). `pre_treatment_fit` is recoverable as `RMSE(treated_pre_trajectory, synthetic_pre_trajectory)`. +- **Trajectories** (`synthetic_pre_trajectory`, `synthetic_post_trajectory`, `treated_pre_trajectory`, `treated_post_trajectory`): retained on results to support plotting and custom fit metrics. `synthetic_pre_trajectory = Y_pre_control @ ω_eff`; `treated_pre_trajectory` is the survey-weighted treated mean (matches the Frank-Wolfe target). With `resid = treated_pre_trajectory - synthetic_pre_trajectory`, `pre_treatment_fit` is recoverable as `RMSE(resid - mean(resid))` (shape-only) and `pre_treatment_level_gap` as `mean(resid)`. - **`get_loo_effects_df()`**: user-facing join of the jackknife leave-one-out pseudo-values (stored in `variance_effects`) to the underlying unit identities. **Unit-level LOO only** — available on the non-survey and pweight-only jackknife paths (classical Algorithm 3: one LOO per unit, first `n_control` positions map to `control_unit_ids`, next `n_treated` to `treated_unit_ids`; `att_loo` is NaN when the zero-sum composed-weight guard fired for that unit; `delta_from_full = att_loo - att`). Under the full-design survey jackknife path (PSU-level LOO with stratum aggregation, Rust & Rao 1996), the underlying replicates are PSU-level rather than unit-level — the accessor raises `NotImplementedError` pointing to `result.variance_effects` for the raw PSU-level replicate array. Dispatch is gated by an explicit `_loo_granularity` flag set at fit-time (`"unit"` vs `"psu"`). Requires `variance_method='jackknife'`; raises `ValueError` otherwise. - **`get_weight_concentration(top_k=5)`**: returns `effective_n = 1/Σω²` (inverse Herfindahl), `herfindahl`, `top_k_share`, `top_k`. Operates on `self.unit_weights` which stores the composed `ω_eff`; for survey-weighted fits the metrics reflect the population-weighted concentration, not the raw Frank-Wolfe solution. - **`in_time_placebo(fake_treatment_periods=None, zeta_omega_override=None, zeta_lambda_override=None)`**: re-slices the pre-window at each fake treatment period and re-fits both ω and λ via Frank-Wolfe. Default sweeps every feasible pre-period (position index `i ≥ 2` so ≥2 pre-fake periods remain for weight estimation, `i ≤ n_pre - 1` so ≥1 post-fake period exists). Credible designs produce near-zero placebo ATTs; departures indicate pre-treatment dynamics the estimator is picking up. - **Note:** Regularization reuses `self.zeta_omega` / `self.zeta_lambda` from the original fit (matches R `synthdid` convention of treating regularization as a property of the fit). `*_override` re-fits with new values. - - **Note:** Infeasibility-only NaN — the method emits NaN for dimensional infeasibility (e.g., survey composition producing zero weight sum on the fake window); Frank-Wolfe non-convergence is not detectable mid-solver, so `pre_fit_rmse` is the user-facing signal for poor refit quality. Passing a `fake_treatment_period` in `post_periods` raises `ValueError` (not a placebo). -- **`sensitivity_to_zeta_omega(zeta_grid=None, multipliers=(0.25, 0.5, 1.0, 2.0, 4.0))`**: re-fits ω at each zeta value on the original pre-window. Default grid is `multipliers * self.zeta_omega` — a 5-point grid spanning 16x from smallest to largest multiplier, symmetric in log space around 1.0. Returns `att`, `pre_fit_rmse`, `max_unit_weight`, `effective_n` per row. + - **Note:** Infeasibility-only NaN — the method emits NaN for dimensional infeasibility (e.g., survey composition producing zero weight sum on the fake window); Frank-Wolfe non-convergence is not detectable mid-solver, so the shape-only `pre_fit_rmse` (pre-window mean gap removed, same definition as `pre_treatment_fit`) is the user-facing signal for poor refit quality. Passing a `fake_treatment_period` in `post_periods` raises `ValueError` (not a placebo). +- **`sensitivity_to_zeta_omega(zeta_grid=None, multipliers=(0.25, 0.5, 1.0, 2.0, 4.0))`**: re-fits ω at each zeta value on the original pre-window. Default grid is `multipliers * self.zeta_omega` — a 5-point grid spanning 16x from smallest to largest multiplier, symmetric in log space around 1.0. Returns `att`, `pre_fit_rmse` (shape-only, same definition as `pre_treatment_fit`; NaN with a single pre-period), `max_unit_weight`, `effective_n` per row. - **Note:** Time weights are held fixed at the original Frank-Wolfe output (`self.time_weights_array`), not re-fit. This isolates sensitivity to `zeta_omega` specifically; sensitivity to `zeta_lambda` is not currently exposed. - **Note:** At `multiplier=1.0` (or `zeta_grid` containing `self.zeta_omega`), the ATT reproduces `self.att` to machine precision with the same seeded draw. diff --git a/docs/methodology/REPORTING.md b/docs/methodology/REPORTING.md index 62faca605..71bcf94f1 100644 --- a/docs/methodology/REPORTING.md +++ b/docs/methodology/REPORTING.md @@ -318,8 +318,9 @@ a library setting. - **Note:** Estimator-native validation surfaces are surfaced rather than duplicated. `SyntheticDiDResults` routes parallel-trends to - `pre_treatment_fit` (the RMSE of the synthetic-control fit on the - pre-period), and routes sensitivity to `in_time_placebo()` + + `pre_treatment_fit` (the shape-only RMSE of the synthetic-control fit on + the pre-period; the constant level gap is reported separately as + `pre_treatment_level_gap`), and routes sensitivity to `in_time_placebo()` + `sensitivity_to_zeta_omega()`. `TROPResults` surfaces factor-model diagnostics (`effective_rank`, `loocv_score`, selected `lambda_*`) under `estimator_native_diagnostics`. `SyntheticControlResults` diff --git a/docs/references.rst b/docs/references.rst index f57aada20..5f7f7347d 100644 --- a/docs/references.rst +++ b/docs/references.rst @@ -164,6 +164,8 @@ Synthetic Control Method - **Abadie, A., Diamond, A., & Hainmueller, J. (2015).** "Comparative Politics and the Synthetic Control Method." *American Journal of Political Science*, 59(2), 495-510. https://doi.org/10.1111/ajps.12116 +- **Abadie, A. (2021).** "Using Synthetic Controls: Feasibility, Data Requirements, and Methodological Aspects." *Journal of Economic Literature*, 59(2), 391-425. https://doi.org/10.1257/jel.20191450 + - **Chernozhukov, V., Wüthrich, K., & Zhu, Y. (2021).** "An Exact and Robust Conformal Inference Method for Counterfactual and Synthetic Controls." *Journal of the American Statistical Association*, 116(536), 1849-1864. https://doi.org/10.1080/01621459.2021.1920957 - **Firpo, S., & Possebom, V. (2018).** "Synthetic Control Method: Inference, Sensitivity Analysis and Confidence Sets." *Journal of Causal Inference*, 6(2), 20160026. https://doi.org/10.1515/jci-2016-0026 diff --git a/docs/tutorials/03_synthetic_did.ipynb b/docs/tutorials/03_synthetic_did.ipynb index 9c456c9f1..4cbb23ebe 100644 --- a/docs/tutorials/03_synthetic_did.ipynb +++ b/docs/tutorials/03_synthetic_did.ipynb @@ -337,7 +337,10 @@ "source": [ "## Pre-treatment Fit\n", "\n", - "A key diagnostic is how well the synthetic control matches the treated units in the pre-treatment period." + "A key diagnostic is how well the synthetic control tracks the treated units' pre-treatment *trend*. SDiD fits the unit weights on centered outcomes, so a constant level gap between the treated units and the synthetic control is not a fit failure - the difference-in-differences step removes it. The library therefore reports two numbers:\n", + "\n", + "- `pre_treatment_fit`: the **shape-only** RMSE of the pre-period residual (treated mean minus synthetic) after removing its mean. Lower is better. The estimator judges it against the same statistic for placebo fits in which control units are treated as if treated (the in-space placebo check of Abadie, Diamond and Hainmueller 2010) and warns when the treated fit is worse than 95% of them (`pre_fit_placebo_pvalue` at or below 0.05).\n", + "- `pre_treatment_level_gap`: the signed mean pre-period gap (treated minus synthetic). Reported for inspection; it never counts as poor fit." ] }, { @@ -346,8 +349,10 @@ "metadata": {}, "outputs": [], "source": [ - "print(f\"Pre-treatment fit (RMSE): {results.pre_treatment_fit:.4f}\")\n", - "print(f\"\\nLower values indicate better fit.\")" + "print(f\"Pre-fit RMSE (shape): {results.pre_treatment_fit:.4f}\")\n", + "print(f\"Pre-fit level gap: {results.pre_treatment_level_gap:.4f}\")\n", + "print(f\"Placebo fit p-value: {results.pre_fit_placebo_pvalue:.3f} (rank of the treated fit among placebo fits, plus-one adjusted)\")\n", + "print(\"\\nLower shape RMSE indicates the synthetic tracks the treated pre-trend more closely.\")" ] }, { @@ -500,7 +505,7 @@ " 'ATT': res.att,\n", " 'SE': res.se,\n", " 'Eff. N controls': eff_n,\n", - " 'Pre-fit RMSE': res.pre_treatment_fit\n", + " 'Pre-fit RMSE (shape)': res.pre_treatment_fit\n", " })\n", "\n", "reg_df = pd.DataFrame(results_list)\n", @@ -599,7 +604,24 @@ { "cell_type": "markdown", "metadata": {}, - "source": "## Summary\n\nKey takeaways for Synthetic DiD:\n\n1. **Best use cases**: Few treated units, many controls, long pre-period\n2. **Unit weights**: Identify which controls are most similar to treated (Frank-Wolfe with sparsification)\n3. **Time weights**: Determine which pre-periods are most informative (Frank-Wolfe on collapsed form)\n4. **Pre-treatment fit**: Lower RMSE indicates better synthetic match\n5. **Inference options**:\n - Placebo (`variance_method=\"placebo\"`, default): Placebo-based variance from controls. Library default (R's default is bootstrap; we deviate for survey availability + perf).\n - Bootstrap (`variance_method=\"bootstrap\"`): Paper-faithful pairs bootstrap re-estimating ω and λ via Frank-Wolfe per draw (Algorithm 2 step 2; matches R's default `vcov`). ~5–30× slower than placebo.\n - Jackknife (`variance_method=\"jackknife\"`): Algorithm 3 — fixed-weight leave-one-out.\n6. **Regularization**: Auto-computed from data noise level by default. Override with `zeta_omega`/`zeta_lambda`.\n\nReference:\n- Arkhangelsky, D., Athey, S., Hirshberg, D. A., Imbens, G. W., & Wager, S. (2021). Synthetic difference-in-differences. American Economic Review, 111(12), 4088-4118." + "source": [ + "## Summary\n", + "\n", + "Key takeaways for Synthetic DiD:\n", + "\n", + "1. **Best use cases**: Few treated units, many controls, long pre-period\n", + "2. **Unit weights**: Identify which controls are most similar to treated (Frank-Wolfe with sparsification)\n", + "3. **Time weights**: Determine which pre-periods are most informative (Frank-Wolfe on collapsed form)\n", + "4. **Pre-treatment fit**: Lower shape-only RMSE (`pre_treatment_fit`) indicates the synthetic tracks the treated pre-trend; a constant level gap (`pre_treatment_level_gap`) is expected and is differenced out\n", + "5. **Inference options**:\n", + " - Placebo (`variance_method=\"placebo\"`, default): Placebo-based variance from controls. Library default (R's default is bootstrap; we deviate for survey availability + perf).\n", + " - Bootstrap (`variance_method=\"bootstrap\"`): Paper-faithful pairs bootstrap re-estimating ω and λ via Frank-Wolfe per draw (Algorithm 2 step 2; matches R's default `vcov`). ~5–30× slower than placebo.\n", + " - Jackknife (`variance_method=\"jackknife\"`): Algorithm 3 — fixed-weight leave-one-out.\n", + "6. **Regularization**: Auto-computed from data noise level by default. Override with `zeta_omega`/`zeta_lambda`.\n", + "\n", + "Reference:\n", + "- Arkhangelsky, D., Athey, S., Hirshberg, D. A., Imbens, G. W., & Wager, S. (2021). Synthetic difference-in-differences. American Economic Review, 111(12), 4088-4118." + ] } ], "metadata": { @@ -609,4 +631,4 @@ }, "nbformat": 4, "nbformat_minor": 4 -} \ No newline at end of file +} diff --git a/docs/tutorials/18_geo_experiments.ipynb b/docs/tutorials/18_geo_experiments.ipynb index 4f5fd5265..68696f7e1 100644 --- a/docs/tutorials/18_geo_experiments.ipynb +++ b/docs/tutorials/18_geo_experiments.ipynb @@ -40,7 +40,13 @@ "cell_type": "markdown", "id": "t18-cell-005", "metadata": {}, - "source": "**Why diff-diff.** The diff-diff library implements Synthetic Difference-in-Differences following [Arkhangelsky et al. (2021)](https://www.aeaweb.org/articles?id=10.1257/aer.20190159) - both the unit weights and the time weights, the placebo standard error procedure from the paper, and full panel-data interpretability. Implementation details and any documented deviations from the R `synthdid` reference are tracked in [`docs/methodology/REGISTRY.md`](https://github.com/igerber/diff-diff/blob/main/docs/methodology/REGISTRY.md).\n\nThis tutorial sits in SDiD's documented sweet spot: a small number of treated markets in a larger pool of donor controls, where basic DiD's averaging doesn't help and you need a counterfactual built specifically for your treated markets. The library's [practitioner decision tree](../practitioner_decision_tree.rst#few-test-markets) puts SyntheticDiD on the \"Few Test Markets\" branch for exactly this reason.\n\nIf you've been using GeoLift, CausalImpact, or rolling your own synthetic control in pandas, this tutorial gives you the canonical SDiD implementation in Python with the diagnostics, inference, and stakeholder packaging in one place." + "source": [ + "**Why diff-diff.** The diff-diff library implements Synthetic Difference-in-Differences following [Arkhangelsky et al. (2021)](https://www.aeaweb.org/articles?id=10.1257/aer.20190159) - both the unit weights and the time weights, the placebo standard error procedure from the paper, and full panel-data interpretability. Implementation details and any documented deviations from the R `synthdid` reference are tracked in [`docs/methodology/REGISTRY.md`](https://github.com/igerber/diff-diff/blob/main/docs/methodology/REGISTRY.md).\n", + "\n", + "This tutorial sits in SDiD's documented sweet spot: a small number of treated markets in a larger pool of donor controls, where basic DiD's averaging doesn't help and you need a counterfactual built specifically for your treated markets. The library's [practitioner decision tree](../practitioner_decision_tree.rst#few-test-markets) puts SyntheticDiD on the \"Few Test Markets\" branch for exactly this reason.\n", + "\n", + "If you've been using GeoLift, CausalImpact, or rolling your own synthetic control in pandas, this tutorial gives you the canonical SDiD implementation in Python with the diagnostics, inference, and stakeholder packaging in one place." + ] }, { "cell_type": "code", @@ -48,10 +54,10 @@ "id": "t18-cell-006", "metadata": { "execution": { - "iopub.execute_input": "2026-04-22T21:17:39.827552Z", - "iopub.status.busy": "2026-04-22T21:17:39.827479Z", - "iopub.status.idle": "2026-04-22T21:17:40.857211Z", - "shell.execute_reply": "2026-04-22T21:17:40.856867Z" + "iopub.execute_input": "2026-09-05T11:55:57.628830Z", + "iopub.status.busy": "2026-09-05T11:55:57.628460Z", + "iopub.status.idle": "2026-09-05T11:55:58.473335Z", + "shell.execute_reply": "2026-09-05T11:55:58.472995Z" } }, "outputs": [], @@ -100,10 +106,10 @@ "id": "t18-cell-009", "metadata": { "execution": { - "iopub.execute_input": "2026-04-22T21:17:40.858677Z", - "iopub.status.busy": "2026-04-22T21:17:40.858527Z", - "iopub.status.idle": "2026-04-22T21:17:40.868516Z", - "shell.execute_reply": "2026-04-22T21:17:40.868282Z" + "iopub.execute_input": "2026-09-05T11:55:58.474808Z", + "iopub.status.busy": "2026-09-05T11:55:58.474708Z", + "iopub.status.idle": "2026-09-05T11:55:58.480790Z", + "shell.execute_reply": "2026-09-05T11:55:58.480556Z" } }, "outputs": [ @@ -155,10 +161,10 @@ "id": "t18-cell-010", "metadata": { "execution": { - "iopub.execute_input": "2026-04-22T21:17:40.869577Z", - "iopub.status.busy": "2026-04-22T21:17:40.869495Z", - "iopub.status.idle": "2026-04-22T21:17:40.874320Z", - "shell.execute_reply": "2026-04-22T21:17:40.874100Z" + "iopub.execute_input": "2026-09-05T11:55:58.481850Z", + "iopub.status.busy": "2026-09-05T11:55:58.481777Z", + "iopub.status.idle": "2026-09-05T11:55:58.486281Z", + "shell.execute_reply": "2026-09-05T11:55:58.486095Z" } }, "outputs": [ @@ -253,10 +259,10 @@ "id": "t18-cell-011", "metadata": { "execution": { - "iopub.execute_input": "2026-04-22T21:17:40.875296Z", - "iopub.status.busy": "2026-04-22T21:17:40.875214Z", - "iopub.status.idle": "2026-04-22T21:17:40.882826Z", - "shell.execute_reply": "2026-04-22T21:17:40.882608Z" + "iopub.execute_input": "2026-09-05T11:55:58.487269Z", + "iopub.status.busy": "2026-09-05T11:55:58.487202Z", + "iopub.status.idle": "2026-09-05T11:55:58.492660Z", + "shell.execute_reply": "2026-09-05T11:55:58.492461Z" } }, "outputs": [ @@ -373,10 +379,10 @@ "id": "t18-cell-012", "metadata": { "execution": { - "iopub.execute_input": "2026-04-22T21:17:40.883777Z", - "iopub.status.busy": "2026-04-22T21:17:40.883709Z", - "iopub.status.idle": "2026-04-22T21:17:40.988502Z", - "shell.execute_reply": "2026-04-22T21:17:40.988246Z" + "iopub.execute_input": "2026-09-05T11:55:58.493652Z", + "iopub.status.busy": "2026-09-05T11:55:58.493586Z", + "iopub.status.idle": "2026-09-05T11:55:58.570567Z", + "shell.execute_reply": "2026-09-05T11:55:58.570329Z" } }, "outputs": [ @@ -432,7 +438,16 @@ "cell_type": "markdown", "id": "t18-cell-015", "metadata": {}, - "source": "Synthetic Difference-in-Differences finds two sets of weights:\n\n1. **Unit weights** ($\\omega_j$): a weighted blend of control markets whose pre-period trajectory matches the treated markets' pre-period trajectory.\n2. **Time weights** ($\\lambda_t$): a weighting of pre-treatment periods that emphasizes the baseline weeks most informative for the comparison.\n\nThe ATT estimator combines both: take the time-weighted average of (treated mean minus unit-weighted control mean), then subtract the same quantity computed in the pre-period. The unit weights make the synthetic control match the treated group; the time weights make the comparison robust to pre-treatment level differences.\n\nThis is the method introduced in [Arkhangelsky, Athey, Hirshberg, Imbens, & Wager (2021)](https://www.aeaweb.org/articles?id=10.1257/aer.20190159). Algorithmic details and any documented deviations from the R `synthdid` reference live in [`docs/methodology/REGISTRY.md`](https://github.com/igerber/diff-diff/blob/main/docs/methodology/REGISTRY.md)." + "source": [ + "Synthetic Difference-in-Differences finds two sets of weights:\n", + "\n", + "1. **Unit weights** ($\\omega_j$): a weighted blend of control markets whose pre-period trajectory matches the treated markets' pre-period trajectory.\n", + "2. **Time weights** ($\\lambda_t$): a weighting of pre-treatment periods that emphasizes the baseline weeks most informative for the comparison.\n", + "\n", + "The ATT estimator combines both: take the time-weighted average of (treated mean minus unit-weighted control mean), then subtract the same quantity computed in the pre-period. The unit weights make the synthetic control match the treated group; the time weights make the comparison robust to pre-treatment level differences.\n", + "\n", + "This is the method introduced in [Arkhangelsky, Athey, Hirshberg, Imbens, & Wager (2021)](https://www.aeaweb.org/articles?id=10.1257/aer.20190159). Algorithmic details and any documented deviations from the R `synthdid` reference live in [`docs/methodology/REGISTRY.md`](https://github.com/igerber/diff-diff/blob/main/docs/methodology/REGISTRY.md)." + ] }, { "cell_type": "code", @@ -440,10 +455,10 @@ "id": "t18-cell-016", "metadata": { "execution": { - "iopub.execute_input": "2026-04-22T21:17:40.989623Z", - "iopub.status.busy": "2026-04-22T21:17:40.989544Z", - "iopub.status.idle": "2026-04-22T21:17:41.228821Z", - "shell.execute_reply": "2026-04-22T21:17:41.228493Z" + "iopub.execute_input": "2026-09-05T11:55:58.571847Z", + "iopub.status.busy": "2026-09-05T11:55:58.571774Z", + "iopub.status.idle": "2026-09-05T11:55:58.807846Z", + "shell.execute_reply": "2026-09-05T11:55:58.807535Z" } }, "outputs": [], @@ -467,10 +482,10 @@ "id": "t18-cell-017", "metadata": { "execution": { - "iopub.execute_input": "2026-04-22T21:17:41.230093Z", - "iopub.status.busy": "2026-04-22T21:17:41.230011Z", - "iopub.status.idle": "2026-04-22T21:17:41.231757Z", - "shell.execute_reply": "2026-04-22T21:17:41.231511Z" + "iopub.execute_input": "2026-09-05T11:55:58.809069Z", + "iopub.status.busy": "2026-09-05T11:55:58.808992Z", + "iopub.status.idle": "2026-09-05T11:55:58.810865Z", + "shell.execute_reply": "2026-09-05T11:55:58.810652Z" } }, "outputs": [ @@ -490,17 +505,19 @@ "Zeta (unit weights): 724.4726\n", "Zeta (time weights): 0.000310\n", "Noise level: 309.5577\n", - "Pre-treatment fit (RMSE): 33.7142\n", + "Pre-fit RMSE (shape): 33.4243\n", + "Pre-fit level gap: 4.4120\n", + "Pre-fit placebo p-value: 0.381\n", "Variance method: placebo\n", "\n", "---------------------------------------------------------------------------\n", "Parameter Estimate Std. Err. t-stat P>|t| \n", "---------------------------------------------------------------------------\n", - "ATT 311.8536 7.3125 42.647 0.0099 **\n", + "ATT 311.8536 7.2659 42.920 0.0099 **\n", "---------------------------------------------------------------------------\n", "\n", - "95% Confidence Interval: [297.5213, 326.1858]\n", - "CV (SE/abs(ATT)): 0.0234\n", + "95% Confidence Interval: [297.6127, 326.0944]\n", + "CV (SE/abs(ATT)): 0.0233\n", "\n", "---------------------------------------------------------------------------\n", " Top Unit Weights (Synthetic Control) \n", @@ -544,7 +561,7 @@ "id": "t18-cell-020", "metadata": {}, "source": [ - "SDiD's interpretability advantage: you can see exactly which control markets are doing the work. The unit weights tell you the weighted blend; the time weights tell you which baseline weeks the method emphasized; the pre-treatment fit RMSE tells you how well the synthetic match worked.\n", + "SDiD's interpretability advantage: you can see exactly which control markets are doing the work. The unit weights tell you the weighted blend; the time weights tell you which baseline weeks the method emphasized; the shape-only pre-treatment fit RMSE tells you how well the synthetic blend tracks the treated pre-*trend* (a constant level gap is differenced out by SDiD and reported separately as `pre_treatment_level_gap`).\n", "\n", "This is not a black box - you can show stakeholders the receipts." ] @@ -555,10 +572,10 @@ "id": "t18-cell-021", "metadata": { "execution": { - "iopub.execute_input": "2026-04-22T21:17:41.232798Z", - "iopub.status.busy": "2026-04-22T21:17:41.232728Z", - "iopub.status.idle": "2026-04-22T21:17:41.235765Z", - "shell.execute_reply": "2026-04-22T21:17:41.235563Z" + "iopub.execute_input": "2026-09-05T11:55:58.812003Z", + "iopub.status.busy": "2026-09-05T11:55:58.811939Z", + "iopub.status.idle": "2026-09-05T11:55:58.814742Z", + "shell.execute_reply": "2026-09-05T11:55:58.814512Z" } }, "outputs": [ @@ -598,10 +615,10 @@ "id": "t18-cell-022", "metadata": { "execution": { - "iopub.execute_input": "2026-04-22T21:17:41.236762Z", - "iopub.status.busy": "2026-04-22T21:17:41.236687Z", - "iopub.status.idle": "2026-04-22T21:17:41.294144Z", - "shell.execute_reply": "2026-04-22T21:17:41.293879Z" + "iopub.execute_input": "2026-09-05T11:55:58.815723Z", + "iopub.status.busy": "2026-09-05T11:55:58.815652Z", + "iopub.status.idle": "2026-09-05T11:55:58.870296Z", + "shell.execute_reply": "2026-09-05T11:55:58.870072Z" } }, "outputs": [ @@ -637,10 +654,10 @@ "id": "t18-cell-023", "metadata": { "execution": { - "iopub.execute_input": "2026-04-22T21:17:41.295247Z", - "iopub.status.busy": "2026-04-22T21:17:41.295156Z", - "iopub.status.idle": "2026-04-22T21:17:41.340650Z", - "shell.execute_reply": "2026-04-22T21:17:41.340382Z" + "iopub.execute_input": "2026-09-05T11:55:58.871418Z", + "iopub.status.busy": "2026-09-05T11:55:58.871323Z", + "iopub.status.idle": "2026-09-05T11:55:58.913722Z", + "shell.execute_reply": "2026-09-05T11:55:58.913499Z" } }, "outputs": [ @@ -695,10 +712,10 @@ "id": "t18-cell-024", "metadata": { "execution": { - "iopub.execute_input": "2026-04-22T21:17:41.341768Z", - "iopub.status.busy": "2026-04-22T21:17:41.341675Z", - "iopub.status.idle": "2026-04-22T21:17:41.343483Z", - "shell.execute_reply": "2026-04-22T21:17:41.343271Z" + "iopub.execute_input": "2026-09-05T11:55:58.914792Z", + "iopub.status.busy": "2026-09-05T11:55:58.914709Z", + "iopub.status.idle": "2026-09-05T11:55:58.916864Z", + "shell.execute_reply": "2026-09-05T11:55:58.916634Z" } }, "outputs": [ @@ -706,22 +723,32 @@ "name": "stdout", "output_type": "stream", "text": [ - "Pre-treatment fit RMSE: 33.71\n", + "Pre-fit RMSE (shape): 33.42\n", + "Pre-fit level gap: 4.41\n", + "Placebo fit p-value: 0.381\n", "\n", - "The library auto-computes a noise level from the data and warns if the\n", - "pre-fit RMSE exceeds the standard deviation of treated pre-period outcomes.\n", - "Our fit is well within that envelope, so the synthetic control is tracking\n", - "the treated pre-trend closely - good fit means a trustworthy post-period effect.\n" + "The shape-only RMSE measures how well the synthetic blend tracks the treated\n", + "pre-trend after removing the constant level gap (SDiD differences that gap out).\n", + "The library compares it with the same statistic for placebo fits of control\n", + "markets treated as if treated, and warns when the treated fit is worse than 95%\n", + "of them. A good pre-fit supports the plausibility of the design; it does not by\n", + "itself validate the counterfactual (weighted parallel trends, no anticipation and\n", + "the absence of post-launch shocks remain assumptions).\n" ] } ], "source": [ - "print(f\"Pre-treatment fit RMSE: {results.pre_treatment_fit:.2f}\")\n", + "print(f\"Pre-fit RMSE (shape): {results.pre_treatment_fit:.2f}\")\n", + "print(f\"Pre-fit level gap: {results.pre_treatment_level_gap:.2f}\")\n", + "print(f\"Placebo fit p-value: {results.pre_fit_placebo_pvalue:.3f}\")\n", "print()\n", - "print(\"The library auto-computes a noise level from the data and warns if the\")\n", - "print(\"pre-fit RMSE exceeds the standard deviation of treated pre-period outcomes.\")\n", - "print(\"Our fit is well within that envelope, so the synthetic control is tracking\")\n", - "print(\"the treated pre-trend closely - good fit means a trustworthy post-period effect.\")" + "print(\"The shape-only RMSE measures how well the synthetic blend tracks the treated\")\n", + "print(\"pre-trend after removing the constant level gap (SDiD differences that gap out).\")\n", + "print(\"The library compares it with the same statistic for placebo fits of control\")\n", + "print(\"markets treated as if treated, and warns when the treated fit is worse than 95%\")\n", + "print(\"of them. A good pre-fit supports the plausibility of the design; it does not by\")\n", + "print(\"itself validate the counterfactual (weighted parallel trends, no anticipation and\")\n", + "print(\"the absence of post-launch shocks remain assumptions).\")" ] }, { @@ -730,10 +757,10 @@ "id": "t18-cell-025", "metadata": { "execution": { - "iopub.execute_input": "2026-04-22T21:17:41.344422Z", - "iopub.status.busy": "2026-04-22T21:17:41.344352Z", - "iopub.status.idle": "2026-04-22T21:17:41.398064Z", - "shell.execute_reply": "2026-04-22T21:17:41.397777Z" + "iopub.execute_input": "2026-09-05T11:55:58.917861Z", + "iopub.status.busy": "2026-09-05T11:55:58.917783Z", + "iopub.status.idle": "2026-09-05T11:55:58.968236Z", + "shell.execute_reply": "2026-09-05T11:55:58.967977Z" } }, "outputs": [ @@ -796,10 +823,10 @@ "id": "t18-cell-026", "metadata": { "execution": { - "iopub.execute_input": "2026-04-22T21:17:41.399142Z", - "iopub.status.busy": "2026-04-22T21:17:41.399060Z", - "iopub.status.idle": "2026-04-22T21:17:41.401344Z", - "shell.execute_reply": "2026-04-22T21:17:41.401110Z" + "iopub.execute_input": "2026-09-05T11:55:58.969290Z", + "iopub.status.busy": "2026-09-05T11:55:58.969225Z", + "iopub.status.idle": "2026-09-05T11:55:58.971631Z", + "shell.execute_reply": "2026-09-05T11:55:58.971412Z" } }, "outputs": [ @@ -823,7 +850,12 @@ " f\"CI ({results.conf_int[0]:.2f}, {results.conf_int[1]:.2f}) does not cover \"\n", " f\"the true effect of 300 - DGP or estimator drift\"\n", ")\n", - "assert results.pre_treatment_fit < 60, f\"Pre-fit RMSE drifted to {results.pre_treatment_fit:.2f}\"\n", + "assert results.pre_treatment_fit < 50, (\n", + " f\"Pre-fit shape RMSE drifted to {results.pre_treatment_fit:.2f}\"\n", + ")\n", + "assert abs(results.pre_treatment_level_gap) < 10, (\n", + " f\"Pre-fit level gap drifted to {results.pre_treatment_level_gap:.2f}\"\n", + ")\n", "print(\"All drift guards passed.\")" ] }, @@ -839,7 +871,17 @@ "cell_type": "markdown", "id": "t18-cell-028", "metadata": {}, - "source": "diff-diff's `SyntheticDiD` supports three standard error methods, and the difference between the two paper-based ones is *what gets resampled* per replication:\n\n- **Placebo SE** (default): permutes which control units are pretended to be \"treated\", then **re-estimates both the unit weights and the time weights** (Frank-Wolfe) on each permutation and recomputes SDiD. The standard deviation of those placebo effects is the SE. This is Algorithm 4 in Arkhangelsky et al. (2021) and matches R's `synthdid::vcov(method=\"placebo\")`.\n- **Bootstrap SE**: pairs-bootstrap resampling of all units with replacement, then **re-estimates both the unit weights and the time weights** via Frank-Wolfe on each resampled panel and recomputes SDiD. This is Algorithm 2 step 2 in Arkhangelsky et al. (2021) and matches R's default `synthdid::vcov(method=\"bootstrap\")` behavior (which rebinds `attr(estimate, \"opts\")` so the renormalized ω is only Frank-Wolfe initialization). Expect ~5–30× slower per fit than placebo (panel-size dependent).\n- **Jackknife SE**: deterministic Algorithm 3 — fixed-weight leave-one-out across all units. Faster than bootstrap; mildly anti-conservative on smaller panels.\n\nBoth bootstrap and placebo re-estimate the weights per replication, so each reflects the full uncertainty in the weighting procedure. They differ in *how* they resample: placebo permutes the control-vs-treated assignment, bootstrap draws with replacement. On exchangeable DGPs the two SEs typically track each other; on small panels with non-exchangeable factor structure (like the marketing geo-experiment here), they can differ in magnitude while still agreeing on significance and CI direction.\n\nAll three methods are configured on the `SyntheticDiD` *constructor*, not on `.fit()`. Use placebo by default (it's the library default; R's default is bootstrap); switch to bootstrap if you want a cross-check from a different resampling protocol; switch to jackknife if you need a deterministic, fast alternative." + "source": [ + "diff-diff's `SyntheticDiD` supports three standard error methods, and the difference between the two paper-based ones is *what gets resampled* per replication:\n", + "\n", + "- **Placebo SE** (default): permutes which control units are pretended to be \"treated\", then **re-estimates both the unit weights and the time weights** (Frank-Wolfe) on each permutation and recomputes SDiD. The standard deviation of those placebo effects is the SE. This is Algorithm 4 in Arkhangelsky et al. (2021) and matches R's `synthdid::vcov(method=\"placebo\")`.\n", + "- **Bootstrap SE**: pairs-bootstrap resampling of all units with replacement, then **re-estimates both the unit weights and the time weights** via Frank-Wolfe on each resampled panel and recomputes SDiD. This is Algorithm 2 step 2 in Arkhangelsky et al. (2021) and matches R's default `synthdid::vcov(method=\"bootstrap\")` behavior (which rebinds `attr(estimate, \"opts\")` so the renormalized ω is only Frank-Wolfe initialization). Expect ~5–30× slower per fit than placebo (panel-size dependent).\n", + "- **Jackknife SE**: deterministic Algorithm 3 — fixed-weight leave-one-out across all units. Faster than bootstrap; mildly anti-conservative on smaller panels.\n", + "\n", + "Both bootstrap and placebo re-estimate the weights per replication, so each reflects the full uncertainty in the weighting procedure. They differ in *how* they resample: placebo permutes the control-vs-treated assignment, bootstrap draws with replacement. On exchangeable DGPs the two SEs typically track each other; on small panels with non-exchangeable factor structure (like the marketing geo-experiment here), they can differ in magnitude while still agreeing on significance and CI direction.\n", + "\n", + "All three methods are configured on the `SyntheticDiD` *constructor*, not on `.fit()`. Use placebo by default (it's the library default; R's default is bootstrap); switch to bootstrap if you want a cross-check from a different resampling protocol; switch to jackknife if you need a deterministic, fast alternative." + ] }, { "cell_type": "code", @@ -847,13 +889,21 @@ "id": "t18-cell-029", "metadata": { "execution": { - "iopub.execute_input": "2026-04-22T21:17:41.402325Z", - "iopub.status.busy": "2026-04-22T21:17:41.402255Z", - "iopub.status.idle": "2026-04-22T21:17:41.662155Z", - "shell.execute_reply": "2026-04-22T21:17:41.661887Z" + "iopub.execute_input": "2026-09-05T11:55:58.972780Z", + "iopub.status.busy": "2026-09-05T11:55:58.972709Z", + "iopub.status.idle": "2026-09-05T11:55:59.259481Z", + "shell.execute_reply": "2026-09-05T11:55:59.259227Z" } }, "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/igerber/diff-diff/diff_diff/synthetic_did.py:1110: UserWarning: Frank-Wolfe did not converge on 100 of 100 valid bootstrap draws (variance_method='bootstrap'). SE is still reported from the final iterate of each draw, but non-convergent draws may be noisier; consider relaxing min_decrease or increasing pre-period length if regularization is already moderate.\n", + " se_n, bootstrap_estimates_n = self._bootstrap_se(\n" + ] + }, { "data": { "text/html": [ @@ -887,9 +937,9 @@ " 0\n", " Placebo (default)\n", " 311.85\n", - " 7.31\n", - " 297.52\n", - " 326.19\n", + " 7.27\n", + " 297.61\n", + " 326.09\n", " \n", " \n", " 1\n", @@ -905,7 +955,7 @@ ], "text/plain": [ " Method ATT SE CI_low CI_high\n", - "0 Placebo (default) 311.85 7.31 297.52 326.19\n", + "0 Placebo (default) 311.85 7.27 297.61 326.09\n", "1 Bootstrap 311.85 4.44 303.15 320.56" ] }, @@ -942,10 +992,10 @@ "id": "t18-cell-030", "metadata": { "execution": { - "iopub.execute_input": "2026-04-22T21:17:41.663265Z", - "iopub.status.busy": "2026-04-22T21:17:41.663187Z", - "iopub.status.idle": "2026-04-22T21:17:41.664937Z", - "shell.execute_reply": "2026-04-22T21:17:41.664722Z" + "iopub.execute_input": "2026-09-05T11:55:59.260593Z", + "iopub.status.busy": "2026-09-05T11:55:59.260521Z", + "iopub.status.idle": "2026-09-05T11:55:59.262361Z", + "shell.execute_reply": "2026-09-05T11:55:59.262154Z" } }, "outputs": [ @@ -971,7 +1021,8 @@ "\n", " * [HIGH] Step 6: Check pre-treatment fit and weight concentration\n", " Why: Synthetic DiD relies on pre-treatment fit to construct weights. Poor fit or highly concentrated unit weights suggest the synthetic control may not approximate the counterfactual well.\n", - " >>> print(f'Pre-treatment fit (RMSE): {results.pre_treatment_fit:.4f}')\n", + " >>> print(f'Pre-fit RMSE (shape): {results.pre_treatment_fit:.4f}')\n", + " >>> print(f'Pre-fit level gap: {results.pre_treatment_level_gap:.4f}')\n", " >>> concentration = results.get_weight_concentration()\n", " >>> print(f\"Effective N: {concentration['effective_n']:.1f}\")\n", " >>> print(f\"Top-5 weight share: {concentration['top_k_share']:.2%}\")\n", @@ -998,7 +1049,7 @@ " >>> sens_df = results.sensitivity_to_zeta_omega()\n", " >>> print(sens_df)\n", "\n", - " * [HIGH] Step 8: Compare with staggered estimators (CS, SA)\n", + " * [HIGH] Step 8: Compare with staggered estimators (CallawaySantAnna, SunAbraham)\n", " Why: SyntheticDiD is for few treated units; compare with staggered estimators if applicable. Use TROP only if factor confounding is suspected (different use case).\n", " >>> from diff_diff import CallawaySantAnna\n", " >>> cs = CallawaySantAnna()\n", @@ -1044,7 +1095,7 @@ ">\n", "> **Sample size and design.** 5 pilot markets, 75 control markets. 12 weeks of weekly data: 6 weeks pre-launch, 6 weeks post-launch. Outcome: weekly conversions per market. Method: Synthetic Difference-in-Differences (Arkhangelsky et al. 2021), the canonical generalization of synthetic control to multi-treated panel settings. The 75 control markets serve as the donor pool that SDiD reweights to construct a counterfactual specific to the 5 pilot markets.\n", ">\n", - "> **Validity evidence.** The synthetic control's pre-treatment fit RMSE is well below the standard deviation of treated pre-period outcomes (the library would warn otherwise), which means the weighted blend of donor markets tracks the treated pre-trend closely. The placebo standard error matches the published Arkhangelsky et al. (2021) method, and we cross-checked with paper-faithful refit bootstrap inference (see Inference and Trustworthiness) — both methods agree on the point estimate (311.85) and on the result's significance, with bootstrap producing a narrower CI than placebo on this small panel (5 treated × 6 pre-periods, factor-model heterogeneity). The estimate is statistically significant under both inference methods, and the placebo 95% CI cleanly covers the true treatment effect on the synthetic data we used to demonstrate the workflow.\n", + "> **Validity evidence.** The synthetic control's shape-only pre-treatment fit RMSE sits inside the range of placebo fits of control markets treated as if treated (the library would warn if it were worse than 95% of them), which means the weighted blend of donor markets tracks the treated pre-trend as well as a typical donor market is tracked; the constant pre-period level gap between the pilot markets and the synthetic blend is reported separately and is differenced out by the method. Good pre-fit supports the plausibility of the design but does not by itself validate the counterfactual. The placebo standard error matches the published Arkhangelsky et al. (2021) method, and we cross-checked with paper-faithful refit bootstrap inference (see Inference and Trustworthiness) — both methods agree on the point estimate (311.85) and on the result's significance, with bootstrap producing a narrower CI than placebo on this small panel (5 treated × 6 pre-periods, factor-model heterogeneity). The estimate is statistically significant under both inference methods, and the placebo 95% CI cleanly covers the true treatment effect on the synthetic data we used to demonstrate the workflow.\n", ">\n", "> **What \"312 conversions per market per week\" means in business terms.** Across 5 pilot markets and 6 weeks, that's roughly 9,400 incremental conversions attributable to the campaign in this small pilot. Translate to your own revenue-per-conversion to compare against the pilot's campaign spend, then use the per-market lift estimate to project what a broader rollout would deliver.\n", ">\n", diff --git a/tests/test_diagnostic_report.py b/tests/test_diagnostic_report.py index 7a4c386b2..c5206519a 100644 --- a/tests/test_diagnostic_report.py +++ b/tests/test_diagnostic_report.py @@ -2026,6 +2026,45 @@ def test_sdid_pt_uses_synthetic_fit_method(self, sdid_fit): assert pt["method"] == "synthetic_fit" assert pt["verdict"] == "design_enforced_pt" assert isinstance(pt.get("pre_treatment_fit_rmse"), float) + assert np.isfinite(pt["pre_treatment_fit_rmse"]) + + def test_sdid_single_pre_period_fit_is_skipped_not_nan(self): + """A 1-pre-period SDiD fit is legal but its shape-only + ``pre_treatment_fit`` is NaN; the PT analogue must report ``skipped`` + (never a narrative "RMSE = nan"), while the raw NaN still flows into + the native-diagnostics data surface by the library's NaN convention.""" + rng = np.random.default_rng(3) + rows = [] + for u in range(10): + is_treated = 1 if u < 2 else 0 + for t in range(3): + rows.append( + { + "unit": u, + "time": t, + "outcome": 2.0 * t + (5.0 if is_treated else 0.0) + rng.normal(0, 0.5), + "treated": is_treated, + } + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + fit = dd.SyntheticDiD().fit( + pd.DataFrame(rows), + outcome="outcome", + treatment="treated", + unit="unit", + time="time", + post_periods=[1, 2], + ) + assert np.isnan(fit.pre_treatment_fit) + report = DiagnosticReport(fit).to_dict() + pt = report["parallel_trends"] + assert pt["status"] == "skipped" + assert "pre_treatment_fit_rmse" not in pt + assert "fewer than 2 pre-periods" in pt["reason"] + native = report["estimator_native_diagnostics"] + assert native["status"] == "ran" + assert np.isnan(native["pre_treatment_fit"]) def test_sdid_native_section_populated(self, sdid_fit): fit, _ = sdid_fit diff --git a/tests/test_estimators.py b/tests/test_estimators.py index 62d30afa3..e0e506721 100644 --- a/tests/test_estimators.py +++ b/tests/test_estimators.py @@ -2666,6 +2666,8 @@ def test_pre_treatment_fit(self, sdid_panel_data): assert results.pre_treatment_fit is not None assert results.pre_treatment_fit >= 0 + assert isinstance(results.pre_treatment_level_gap, float) + assert np.isfinite(results.pre_treatment_level_gap) def test_summary_output(self, sdid_panel_data, ci_params): """Test that summary produces string output.""" @@ -2705,6 +2707,10 @@ def test_to_dict(self, sdid_panel_data, ci_params): assert "n_pre_periods" in result_dict assert "n_post_periods" in result_dict assert "pre_treatment_fit" in result_dict + assert "pre_treatment_level_gap" in result_dict + assert np.isfinite(result_dict["pre_treatment_level_gap"]) + assert "pre_fit_placebo_pvalue" in result_dict + assert 0.0 < result_dict["pre_fit_placebo_pvalue"] <= 1.0 def test_to_dataframe(self, sdid_panel_data, ci_params): """Test conversion to DataFrame.""" diff --git a/tests/test_methodology_sdid.py b/tests/test_methodology_sdid.py index ff5089a10..39beddd60 100644 --- a/tests/test_methodology_sdid.py +++ b/tests/test_methodology_sdid.py @@ -14,6 +14,7 @@ import pandas as pd import pytest +from diff_diff.results import SyntheticDiDResults from diff_diff.synthetic_did import SyntheticDiD from diff_diff.utils import ( _compute_noise_level, @@ -1400,10 +1401,25 @@ def test_jackknife_same_att_as_placebo(self): assert abs(res_jk.att - res_pl.att) < 1e-10 def test_jackknife_n_bootstrap_ignored(self): - """n_bootstrap=1 should not raise for jackknife (it's ignored).""" + """n_bootstrap is ignored by jackknife VARIANCE estimation (n_bootstrap=1 + must not raise; ATT and SE are identical), but it still caps the + pre-fit placebo reference (min(n_bootstrap, 20) draws), so the + diagnostic fields differ and the warning is unreachable at 1 draw.""" sdid = SyntheticDiD(variance_method="jackknife", n_bootstrap=1) assert sdid.n_bootstrap == 1 assert sdid.variance_method == "jackknife" + df = _make_panel(n_control=15, n_treated=3, seed=42) + kw = dict(outcome="outcome", treatment="treated", unit="unit", time="period") + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res1 = SyntheticDiD(variance_method="jackknife", n_bootstrap=1, seed=42).fit(df, **kw) + res20 = SyntheticDiD(variance_method="jackknife", n_bootstrap=20, seed=42).fit(df, **kw) + assert res1.att == res20.att + assert res1.se == res20.se + assert res1.n_bootstrap is None and res20.n_bootstrap is None + assert len(res1.pre_fit_placebo_rmse) == 1 + assert len(res20.pre_fit_placebo_rmse) == 20 + assert res1.pre_fit_placebo_pvalue >= 0.5 # 1 draw: p >= 1/2, never warns def test_jackknife_n_bootstrap_none_in_results(self): """Results should have n_bootstrap=None for jackknife.""" @@ -2506,13 +2522,78 @@ def test_balanced_panel_passes(self): class TestPreTreatmentFitWarning: """Test that poor pre-treatment fit emits a warning.""" + @staticmethod + def _panel(treated_fn, control_fn, *, n_treated=2, n_control=8, T=8, noise=0.5, seed=42): + """Panel where treated/control outcomes follow ``f(t) + N(0, noise)``.""" + rng = np.random.default_rng(seed) + rows = [] + for u in range(n_treated + n_control): + is_treated = 1 if u < n_treated else 0 + f = treated_fn if is_treated else control_fn + for t in range(T): + rows.append( + { + "unit": u, + "time": t, + "outcome": f(t) + rng.normal(0, noise), + "treated": is_treated, + } + ) + return pd.DataFrame(rows) + + @staticmethod + def _fit(data, post_periods): + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + res = SyntheticDiD().fit( + data, + outcome="outcome", + treatment="treated", + unit="unit", + time="time", + post_periods=post_periods, + ) + fit_warnings = [x for x in w if "Pre-treatment fit is poor" in str(x.message)] + return res, fit_warnings + def test_poor_fit_emits_warning(self): - """Treated units at very different level from controls should warn.""" + """Treated series trends while controls are flat: no simplex weighting + of flat controls can follow the trend, so the shape-only residual is + many multiples of the control noise level and the warning fires. + (Under the retired 1x treated-SD rule this design could NOT fire for + a shape-only RMSE: a flat synthetic gives exactly the population SD.) + """ + data = self._panel(lambda t: 2.0 * t, lambda t: 0.0) + res, fit_warnings = self._fit(data, post_periods=[6, 7]) + assert len(fit_warnings) >= 1, "Expected poor pre-treatment fit warning" + assert res.pre_fit_placebo_pvalue <= 0.05 + # Worse than every one of the 20 placebo fits: p = 1 / 21. + assert res.pre_fit_placebo_pvalue == pytest.approx(1.0 / 21.0) + assert res.pre_treatment_fit > np.max(res.pre_fit_placebo_rmse) + msg = str(fit_warnings[0].message) + assert "shape-only" in msg and "pre_treatment_level_gap" in msg and "placebo" in msg + + def test_level_offset_does_not_warn(self): + """The reported scenario: treated = control trend + constant. The FW + objective is column-centered, so the level gap is not a fit failure; + the RMSE must be small, the gap reported, and the ATT recovered.""" + data = self._panel(lambda t: 2.0 * t + 50.0, lambda t: 2.0 * t) + data.loc[(data["treated"] == 1) & (data["time"] >= 6), "outcome"] += 5.0 + res, fit_warnings = self._fit(data, post_periods=[6, 7]) + assert len(fit_warnings) == 0, f"Unexpected fit warning: {fit_warnings[0].message}" + assert res.pre_treatment_level_gap == pytest.approx(50.0, abs=2.0) + assert res.pre_fit_placebo_pvalue > 0.05 + assert res.pre_treatment_fit < 0.05 * abs(res.pre_treatment_level_gap) + assert res.att == pytest.approx(5.0, abs=1.5) + + def test_flat_noise_treated_does_not_warn(self): + """The pre-v3.11.2 'poor fit' fixture (treated ~100, controls ~10, both + flat) is a textbook GOOD SDID design: a pure level offset. It must no + longer warn, and the offset must surface as the level gap.""" np.random.seed(42) rows = [] for u in range(10): is_treated = 1 if u < 2 else 0 - # Large level difference: treated ~100, control ~10 level = 100.0 if is_treated else 10.0 for t in range(8): rows.append( @@ -2523,22 +2604,213 @@ def test_poor_fit_emits_warning(self): "treated": is_treated, } ) - data = pd.DataFrame(rows) - sdid = SyntheticDiD() + res, fit_warnings = self._fit(pd.DataFrame(rows), post_periods=[6, 7]) + assert len(fit_warnings) == 0, f"Unexpected fit warning: {fit_warnings[0].message}" + assert res.pre_treatment_level_gap == pytest.approx(90.0, abs=1.0) + + def test_single_pre_period_fit_is_nan(self): + """Shape is undefined with one pre-period: NaN RMSE, finite level gap, + no warning; sensitivity_to_zeta_omega reaches the same NaN branch.""" + data = self._panel(lambda t: 2.0 * t + 5.0, lambda t: 2.0 * t, T=3) + res, fit_warnings = self._fit(data, post_periods=[1, 2]) + assert len(fit_warnings) == 0 + assert np.isnan(res.pre_treatment_fit) + assert np.isfinite(res.pre_treatment_level_gap) + assert res.pre_fit_placebo_rmse is None and res.pre_fit_placebo_pvalue is None + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + sens = res.sensitivity_to_zeta_omega() + assert sens["pre_fit_rmse"].isna().all() + assert np.isfinite(sens["att"]).all() + + def test_two_pre_periods_single_control_no_warning(self): + """With a single control there is no pseudo-control left once it is + treated as if treated, so no placebo reference exists and the fit + warning must stay silent (the fit itself is legal).""" + data = self._panel( + lambda t: 2.0 * t + 3.0, lambda t: 2.0 * t, n_treated=1, n_control=1, T=4, noise=0.05 + ) with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always") - sdid.fit( + res = SyntheticDiD().fit( data, outcome="outcome", treatment="treated", unit="unit", time="time", - post_periods=[6, 7], + post_periods=[2, 3], ) - fit_warnings = [x for x in w if "Pre-treatment fit is poor" in str(x.message)] - assert ( - len(fit_warnings) >= 1 - ), "Expected warning about poor pre-treatment fit but none was raised" + # The single-control placebo variance path emits its own + # "Not enough control units" warning; only the fit warning is asserted. + fit_warnings = [x for x in w if "Pre-treatment fit" in str(x.message)] + assert res.noise_level == 0.0 + assert np.isfinite(res.pre_treatment_fit) + # No pseudo-control remains once the single control is pseudo-treated. + assert res.pre_fit_placebo_rmse is None and res.pre_fit_placebo_pvalue is None + assert len(fit_warnings) == 0, f"Unexpected fit warning: {fit_warnings[0].message}" + + def test_noiseless_controls_still_warn(self): + """Noiseless, exactly parallel controls are fit exactly by every + placebo draw (placebo RMSE 0), so a trending treated series is worse + than all of them and warns (this input warned under the pre-v3.11.2 + level rule too).""" + rows = [] + for u in range(10): + is_treated = 1 if u < 2 else 0 + for t in range(8): + y = 2.0 * t if is_treated else 10.0 + 3.0 * u + rows.append({"unit": u, "time": t, "outcome": y, "treated": is_treated}) + res, fit_warnings = self._fit(pd.DataFrame(rows), post_periods=[6, 7]) + assert res.noise_level == 0.0 + assert res.pre_treatment_fit > 1e-6 + # Flat placebo sets are fit exactly by flat controls: every placebo + # RMSE is 0, so the trending treated series is worse than all of them. + assert np.max(res.pre_fit_placebo_rmse) == pytest.approx(0.0, abs=1e-10) + assert res.pre_fit_placebo_pvalue == pytest.approx(1.0 / 21.0) + assert len(fit_warnings) >= 1, "Expected poor pre-treatment fit warning" + + def test_treated_level_shift_invariance(self): + """Adding a constant to the TREATED units only leaves the centered FW + objective and the control-derived normalisation unchanged, so the + shape-only RMSE must be invariant on all three surfaces (fit, + sensitivity_to_zeta_omega, in_time_placebo) while the level gap + shifts by exactly the constant. Under the retired level formula every + one of these moved by an amount of order the offset.""" + df0 = _make_panel(seed=17) + df1 = df0.copy() + df1.loc[df1["treated"] == 1, "outcome"] += 25.0 + kw = dict(outcome="outcome", treatment="treated", unit="unit", time="period") + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res0 = SyntheticDiD(variance_method="jackknife", seed=17).fit(df0, **kw) + res1 = SyntheticDiD(variance_method="jackknife", seed=17).fit(df1, **kw) + sens0, sens1 = res0.sensitivity_to_zeta_omega(), res1.sensitivity_to_zeta_omega() + itp0, itp1 = res0.in_time_placebo(), res1.in_time_placebo() + assert res1.pre_treatment_fit == pytest.approx(res0.pre_treatment_fit, rel=1e-8) + assert res1.pre_treatment_level_gap - res0.pre_treatment_level_gap == pytest.approx( + 25.0, abs=1e-8 + ) + np.testing.assert_allclose( + sens1["pre_fit_rmse"].to_numpy(), sens0["pre_fit_rmse"].to_numpy(), rtol=1e-8 + ) + np.testing.assert_allclose( + itp1["pre_fit_rmse"].to_numpy(), itp0["pre_fit_rmse"].to_numpy(), rtol=1e-8 + ) + # The multiplier-1.0 grid point re-fits at the fit-time zeta on the + # same window, so it must reproduce the fit-time statistic exactly. + row = sens0[np.isclose(sens0["zeta_omega"], res0.zeta_omega)] + assert len(row) == 1 + assert float(row["pre_fit_rmse"].iloc[0]) == pytest.approx(res0.pre_treatment_fit, rel=1e-8) + + def test_pre_fit_placebo_reference_fields(self): + """The reference distribution is capped at min(n_bootstrap, 20) draws, + finite, read-only, identical across variance methods for the same + seed (private RNG stream), and silenced below 19 draws.""" + data = self._panel(lambda t: 2.0 * t + 50.0, lambda t: 2.0 * t) + kw = dict(outcome="outcome", treatment="treated", unit="unit", time="time") + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res_p = SyntheticDiD(variance_method="placebo", seed=7).fit( + data, post_periods=[6, 7], **kw + ) + res_j = SyntheticDiD(variance_method="jackknife", seed=7).fit( + data, post_periods=[6, 7], **kw + ) + res_small = SyntheticDiD(variance_method="placebo", n_bootstrap=10, seed=7).fit( + data, post_periods=[6, 7], **kw + ) + assert len(res_p.pre_fit_placebo_rmse) == 20 + assert np.all(np.isfinite(res_p.pre_fit_placebo_rmse)) + assert np.all(res_p.pre_fit_placebo_rmse >= 0) + assert not res_p.pre_fit_placebo_rmse.flags.writeable + np.testing.assert_array_equal(res_j.pre_fit_placebo_rmse, res_p.pre_fit_placebo_rmse) + assert res_j.pre_fit_placebo_pvalue == res_p.pre_fit_placebo_pvalue + assert 0.0 < res_p.pre_fit_placebo_pvalue <= 1.0 + assert res_p.to_dict()["pre_fit_placebo_pvalue"] == res_p.pre_fit_placebo_pvalue + assert "Pre-fit placebo p-value" in res_p.summary() + # n_bootstrap=10 -> 10 draws -> p >= 1/11 > 0.05: the rule cannot fire. + assert len(res_small.pre_fit_placebo_rmse) == 10 + assert res_small.pre_fit_placebo_pvalue >= 1.0 / 11.0 + + def test_large_common_offset_location_invariance(self): + """CI review P1 on PR #818: the fit-time diagnostics are computed on + the normalized arrays, so a large common outcome level (1e9) leaves + the shape RMSE, level gap, placebo p-value and warning decision + unchanged, and the fit-time RMSE equals the multiplier-1.0 + sensitivity row exactly.""" + base = self._panel(lambda t: 2.0 * t + 5.0, lambda t: 2.0 * t, noise=1.0) + shifted = base.copy() + shifted["outcome"] += 1.0e9 + kw = dict(outcome="outcome", treatment="treated", unit="unit", time="time") + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res0 = SyntheticDiD(seed=3).fit(base, post_periods=[6, 7], **kw) + res1 = SyntheticDiD(seed=3).fit(shifted, post_periods=[6, 7], **kw) + sens1 = res1.sensitivity_to_zeta_omega() + assert res1.pre_treatment_fit == pytest.approx(res0.pre_treatment_fit, rel=1e-6) + assert res1.pre_treatment_level_gap == pytest.approx(res0.pre_treatment_level_gap, rel=1e-6) + assert res1.pre_fit_placebo_pvalue == res0.pre_fit_placebo_pvalue + np.testing.assert_allclose(res1.pre_fit_placebo_rmse, res0.pre_fit_placebo_rmse, rtol=1e-6) + row = sens1[np.isclose(sens1["zeta_omega"], res1.zeta_omega)] + assert len(row) == 1 + assert float(row["pre_fit_rmse"].iloc[0]) == pytest.approx( + res1.pre_treatment_fit, rel=1e-10 + ) + + def test_legacy_pickle_state_migrates_level_rmse(self): + """CI review P1 on PR #818: a results object pickled before v3.11.2 + carries the LEVEL-inclusive RMSE in ``pre_treatment_fit`` and none of + the new fields. ``__setstate__`` must recompute the shape-only RMSE + and the level gap from the stored trajectories (never relabel the + stale value), default the placebo fields to None, and leave summary() + callable; with no trajectories the stale value is cleared.""" + data = self._panel(lambda t: 2.0 * t + 50.0, lambda t: 2.0 * t) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = SyntheticDiD(seed=5).fit( + data, outcome="outcome", treatment="treated", unit="unit", time="time" + ) + resid = res.treated_pre_trajectory - res.synthetic_pre_trajectory + legacy_level_rmse = float(np.sqrt(np.mean(resid**2))) + state = res.__getstate__() + for key in ("pre_treatment_level_gap", "pre_fit_placebo_rmse", "pre_fit_placebo_pvalue"): + state.pop(key) + state["pre_treatment_fit"] = legacy_level_rmse # what an old pickle stores + assert legacy_level_rmse > 10 * res.pre_treatment_fit # level-dominated + + revived = SyntheticDiDResults.__new__(SyntheticDiDResults) + revived.__setstate__(state) + assert revived.pre_treatment_fit == pytest.approx(res.pre_treatment_fit, rel=1e-12) + assert revived.pre_treatment_level_gap == pytest.approx( + res.pre_treatment_level_gap, rel=1e-12 + ) + assert revived.pre_fit_placebo_rmse is None and revived.pre_fit_placebo_pvalue is None + summary = revived.summary() + assert "Pre-fit RMSE (shape)" in summary and "placebo p-value" not in summary + assert revived.to_dict()["pre_fit_placebo_pvalue"] is None + + # No trajectories stored (older still): clear rather than relabel. + bare = dict(state) + for key in ( + "treated_pre_trajectory", + "synthetic_pre_trajectory", + "treated_post_trajectory", + "synthetic_post_trajectory", + ): + bare.pop(key) + revived_bare = SyntheticDiDResults.__new__(SyntheticDiDResults) + revived_bare.__setstate__(bare) + assert revived_bare.pre_treatment_fit is None + assert revived_bare.pre_treatment_level_gap is None + assert "Pre-fit RMSE" not in revived_bare.summary() + + # A current pickle round-trips unchanged. + import pickle + + rt = pickle.loads(pickle.dumps(res)) + assert rt.pre_treatment_fit == res.pre_treatment_fit + assert rt.pre_fit_placebo_pvalue == res.pre_fit_placebo_pvalue + np.testing.assert_array_equal(rt.pre_fit_placebo_rmse, res.pre_fit_placebo_rmse) def test_good_fit_no_warning(self): """Parallel trends data with similar levels should not warn.""" @@ -2813,10 +3085,11 @@ def test_pre_fit_rmse_recoverable(self): df = _make_panel(seed=17) sdid = SyntheticDiD(variance_method="jackknife", seed=17) res = sdid.fit(df, outcome="outcome", treatment="treated", unit="unit", time="period") - rmse = float( - np.sqrt(np.mean((res.treated_pre_trajectory - res.synthetic_pre_trajectory) ** 2)) - ) + resid = res.treated_pre_trajectory - res.synthetic_pre_trajectory + level_gap = float(np.mean(resid)) + rmse = float(np.sqrt(np.mean((resid - level_gap) ** 2))) assert abs(rmse - res.pre_treatment_fit) < 1e-10 + assert abs(level_gap - res.pre_treatment_level_gap) < 1e-10 class TestLooEffectsDf: